/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ParseAbundances parse and read in composition as set by abundances command */
#include "cddefines.h"
#include "grains.h"
#include "grainvar.h"
#include "abund.h"
#include "phycon.h"
#include "called.h"
#include "elementnames.h"
#include "input.h"
#include "parser.h"

void ParseAbundances(Parser &p)
			/* following set true by grains command,
			* so this will not set more grains if grains already on. */
{
	bool lgLog;
	long int i;
	double absav[LIMELM], 
	  chk;
	GrainPar gp;

	DEBUG_ENTRY( "ParseAbundances()" );

	if( p.nMatch("STAR") )
	{
		/* Fred Hamann's star burst galaxy mixture -- includes a number which isn't an abundance */
		abund_starburst(p);
		return;
	}

	absav[0] = p.FFmtRead();
	/* this branch at least one number on the line - an abundance */
	if( !p.lgEOL() )
	{
		absav[1] = p.FFmtRead();
		if( p.lgEOL() )
		{
			/* this branch, we found one, but not two, numbers -
			 * must be one of a few special cases */
			if( p.nMatch(" ALL") )
			{
				/* special option, all abundances will be this number */
				if( absav[0] <= 0. )
				{
					absav[0] = pow(10.,absav[0]);
				}
				for( i=1; i < LIMELM; i++ )
				{
					abund.solar[i] = (realnum)absav[0];
				}

			}
			else if( p.nMatch("OLD ") && p.nMatch("SOLA") )
			{
				i = (int)absav[0];
				/* old solar - better be the number "84" on the line */
				if( i!=84 )
				{
					fprintf( ioQQQ, 
						" The only old abundance set I have is for version 84 - %3ld was requested.  Sorry.\n", 
						i );
					cdEXIT(EXIT_FAILURE);
				}
				for( i=1; i < LIMELM; i++ )
				{
					/* these are the old abundances */
					abund.SolarSave[i] = abund.OldSolar84[i];
					abund.solar[i] = abund.OldSolar84[i];
				}
			}
			else if( p.nMatch("GASS10"))
			{
				/* Grevesse, Asplund, Sauval, and Scott solar */
				/* >>refer	solar	abund	Grevesse, N., Asplund, M., Sauval, A. J., & Scott, P. 2010, Ap&SS, 48 */
				for( i=1; i < LIMELM; i++ )
				{
					/* these are the GASS10 solar abundances */
					abund.SolarSave[i] = abund.GASS10[i];
					abund.solar[i] = abund.GASS10[i];
				}
			}
			else
			{
				fprintf( ioQQQ, 
					" I did not recognize a sub-keyword - options are ALL, OLD SOLAR 84, and GASS10.  Sorry.\n");
				cdEXIT(EXIT_FAILURE);
			}

			/* normal return */
			return;
		}

		/* we get here if there is a second number - read in all abundances */
		for( i=2; i < abund.npSolar; i++ )
		{
			absav[i] = p.FFmtRead();
			if( p.lgEOL() )
			{
				/*   read CONTINUE line if J scanned off end before reading all abund */
				do 
				{
					p.getline();
					if( p.m_lgEOF )
					{
						fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
									abund.npSolar, i );
						cdEXIT(EXIT_FAILURE);
					}
				} while( p.isComment() ); 

				p.echo();

				if( p.strcmp("CONT") != 0 )
				{
					fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
						abund.npSolar, i );
					cdEXIT(EXIT_FAILURE);
				}
				else
				{
					absav[i] = p.FFmtRead();
					if( p.lgEOL() )
					{
						fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
							abund.npSolar, i);
						cdEXIT(EXIT_FAILURE);
					}
				}
			}
		}

		/* fell through to here after reading in N abundances for N elements
		 * check that there are no more abundances on the line - that would
		 * be an error - a typo */
		chk = p.FFmtRead();
		if( !p.lgEOL() || (chk!=0.) )
		{
			/* got another number, not lgEOL, so too many numbers */
			fprintf( ioQQQ, " There were more than %ld abundances entered\n", 
				abund.npSolar );
			fprintf( ioQQQ, " Could there have been a typo somewhere?\n" );
		}

		/* are numbers scale factors, or log of abund rel to H?? */
		lgLog = false;
		for( i=0; i < abund.npSolar; i++ )
			if( absav[i] < 0. )
				lgLog = true;

		if( lgLog )
		{
			/* entered as log of number rel to hydrogen */
			for( i=0; i < abund.npSolar; i++ )
				abund.solar[abund.ipSolar[i]-1] = (realnum)pow(10.,absav[i]);
		}
		else
		{
			/* scale factors relative to solar */
			for( i=0; i < abund.npSolar; i++ )
				abund.solar[abund.ipSolar[i]-1] *= (realnum)absav[i];
		}

		/* check whether the abundances are reasonable */
		for( i=1; i < LIMELM; i++ )
		{
			if( abund.solar[i] > 0.2 )
			{
				fprintf( ioQQQ, " Is an abundance of %.3e relative to H reasonable for %2.2s?\n", 
					abund.solar[i], elementnames.chElementSym[i] );
			}
		}
		return;
	}

	/* following set of clauses -  no numbers on the line - 
	 * use one of several stored abundance sets */
	if( p.nMatch(" AGB") || p.nMatch("AGB ") || p.nMatch("PLAN") )
	{
		/* AGB/planetary nebula abundances */
		/* only turn on grains if "no grains" is not present */
		if( !p.nMatch("NO GR") )
		{
			/* turn on grains if not already done */
			if( !p.m_lgDSet )
			{
				/* return bins allocated by previous abundances ... commands */
				gv.clear();
				/* now allocate new grain bins */
				gp.dep = 1.;
				// standard dust to gas ratio
				gp.nDustFunc = DF_STANDARD;
				gp.lgRequestQHeating = true;
				gp.lgForbidQHeating = false;
				gp.lgGreyGrain = false;

				/* NO QHEAT option to turn off quantum heating for grains */
				if( p.nMatch("NO QH") ) 
				{
					gp.lgForbidQHeating = true;
					gp.lgRequestQHeating = false;
					phycon.lgPhysOK = false;
				}

				/* actually set up the grains */
				mie_read_opc("graphite_ism_10.opc",gp);
				mie_read_opc("silicate_ism_10.opc",gp);
			}
		}

		for( i=0; i < LIMELM; i++ )
		{
			abund.solar[i] = abund.apn[i];
			if( !abund.lgElmONapn[i] )
			{
				/* turn off elements - do this way to make sure,
				* that iso sequence limits are handled properly */
				char chDUMMY[INPUT_LINE_LENGTH];
				sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
				p.setline(chDUMMY);				
				ParseElement( p );
			}
		}
	}

	else if( p.nMatch("CAME") )
	{
		/* mix from Cameron 1982, "Essays on Nuclear Astrophysics" */
		/* now turn off the heavy elements */
		for( i=0; i < LIMELM; i++ )
			abund.solar[i] = abund.camern[i];
	}

	else if( p.nMatch("CRAB") )
	{
		// Crab Nebula "typical" abundances
		for( i=0; i < LIMELM; i++ )
		{
			abund.solar[i] = abund.aCrab[i];
			if( !abund.lgElmONaCrab[i] )
			{
				/* turn off elements - do this way to make sure,
				* that iso sequence limits are handled properly */
				char chDUMMY[INPUT_LINE_LENGTH];
				sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
				p.setline(chDUMMY);
				ParseElement( p );
			}
		}
	}
	else if( p.nMatch("HII ") || p.nMatch("H II") || p.nMatch("ORIO") )
	{
		/* H II region abundances - turn on grains by default
		 * "no grains" to not do this */
		if( !p.nMatch("NO GR") )
		{
			/* option to turn on grains */
			if( !p.m_lgDSet )
			{
				/* return bins allocated by previous abundances ... commands */
				gv.clear();
				/* now allocate new grain bins */
				gp.dep = 1.;
				// standard dust to gas ratio
				gp.nDustFunc = DF_STANDARD;
				gp.lgRequestQHeating = true;
				gp.lgForbidQHeating = false;
				gp.lgGreyGrain = false;

				/* NO QHEAT option to turn off quantum heating for grains */
				if( p.nMatch("NO QH") ) 
				{
					gp.lgForbidQHeating = true;
					gp.lgRequestQHeating = false;
					phycon.lgPhysOK = false;
				}
				/* This scales the Orion grain abundance so that the observed 
				 * dust to gas ratio that Cloudy predicts is in agreement with
				 * that observed in the Veil,
				 *>>refer	grain	Abel, N., Brogan, C., Ferland, G., O'Dell, C.R., 
				 *>>refercon	Shaw, G., Troland, T., 2004, ApJ, submitted */
				gp.dep *= 0.85;

				mie_read_opc("graphite_orion_10.opc",gp);
				mie_read_opc("silicate_orion_10.opc",gp);
			}
		}

		for( i=0; i < LIMELM; i++ )
		{
			abund.solar[i] = abund.ahii[i];
			if( !abund.lgElmONahii[i] )
			{
				/* turn off elements - do this way to make sure,
				* that iso sequence limits are handled properly */
				char chDUMMY[INPUT_LINE_LENGTH];
				sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
				p.setline(chDUMMY);
				ParseElement( p );
			}
		}
	}

	else if( p.nMatch("ISM ") || p.nMatch(" ISM") )
	{
		/* ISM abundances from Cowie and Songaila Ann Rev '86 */
		/* only turn on grains if "no grains" is not present */
		if( !p.nMatch("NO GR") )
		{
			if( !p.m_lgDSet )
			{
				/* return bins allocated by previous abundances ... commands */
				gv.clear();
				/* now allocate new grain bins */
				gp.dep = 1.;
				// standard dust to gas ratio
				gp.nDustFunc = DF_STANDARD;
				gp.lgRequestQHeating = true;
				gp.lgForbidQHeating = false;
				gp.lgGreyGrain = false;

				/* NO QHEAT option to turn off quantum heating */
				if( p.nMatch("NO QH") ) 
				{
					gp.lgForbidQHeating = true;
					gp.lgRequestQHeating = false;
					phycon.lgPhysOK = false;
				}

				/* actually set up the grains */
				mie_read_opc("graphite_ism_10.opc",gp);
				mie_read_opc("silicate_ism_10.opc",gp);
			}
		}

		for( i=0; i < LIMELM; i++ )
		{
			abund.solar[i] = abund.aism[i];
			if( !abund.lgElmONaism[i] )
			{
				/* turn off elements - do this way to make sure,
				 * that iso sequence limits are handled properly */
				char chDUMMY[INPUT_LINE_LENGTH];
				sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
				p.setline(chDUMMY);
				ParseElement( p );
			}
		}
	}

	else if( p.nMatch("NOVA") )
	{
		/* Nova Cyg abundances */
		for( i=0; i < LIMELM; i++ )
			abund.solar[i] = abund.anova[i];
	}

	else if( p.nMatch("PRIM") )
	{
		/* roughly primordial abundances: He/H=.072 */
		for( i=0; i < 4; i++ )
		{
			abund.solar[i] = abund.aprim[i];
		}

		/* now turn off the heavy elements */
		for( i=4; i < LIMELM; i++ )
		{
			/* turn off heavy elements - do this way to make sure,
			 * that H-like and He-like level limits are handled properly */
			char chDUMMY[INPUT_LINE_LENGTH];
			sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
			p.setline(chDUMMY);
			ParseElement( p );
		}
	}

	else
	{
		fprintf( ioQQQ, 
			" ABUNDances must have GASS10, PLAN, H II, CAMERON, CRAB, NOVA, ALL, STARBURST, OLD SOLAR 84 or PRIMORDIAL.  Sorry.\n" );
		cdEXIT(EXIT_FAILURE);
	}
	return;
}