/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * 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 "abund.h"
#include "phycon.h"
#include "called.h"
#include "elementnames.h"
#include "input.h"
#include "parse.h"

void ParseAbundances(char *chCard,
			/* following set TRUE by grains command,
			* so this will not set more grains if grains already on. */
		     int lgDSet)
{
	int lgEOF, 
	  lgEOL, 
	  lgLog;
	long int i, 
	  j, 
	  k;
	double absav[LIMELM], 
	  chk;
	GrainPar gp;


#	ifdef DEBUG_FUN
	fputs( "<+>ParseAbundances()\n", debug_fp );
#	endif

	j = 5;
	absav[0] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
	/* check whether it scanned off end of line on first try
	 * >>chng 97 jun 8, get rid of go to
	 * if( lgEOL ) go to10 */
	/* at least one number on the line if it passes this test */
	if( !lgEOL )
	{
		absav[1] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			/* this branch, we found one, but not two, numbers -
			 * must be one of a few special cases */
			if( lgMatch(" ALL",chCard) )
			{
				/* 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] = (float)absav[0];
				}
				
			}
			else if( lgMatch("OLD ",chCard) && lgMatch("SOLA",chCard) )
			{
				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 );
					puts( "[Stop in ParseAbundances]" );
					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( lgMatch("STAR",chCard) )
			{
				/* Fred Hamonn's starburst galaxy mixture */
				abund_starburst(chCard);
			}
			else
			{
				fprintf( ioQQQ, 
					" I did not recognize a sub-keyword - options are ALL and OLD SOLAR 84.  Sorry.\n");
				puts( "[Stop in ParseAbundances]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* normal return */
#			ifdef DEBUG_FUN
			fputs( " <->ParseAbundances()\n", debug_fp );
#			endif
			return;
		}

		/* do this is there was a second number */
		for( i=2; i < abund.npSolar; i++ )
		{
			absav[i] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				/*   read CONTINUE line if J scanned off end before reading all abund */
				input_readarray(chCard,&lgEOF);
				if( lgEOF )
				{
					fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
						abund.npSolar, i );
					puts( "[Stop in ParseAbundances]" );
					cdEXIT(EXIT_FAILURE);
				}

				/*   option to ignore all lines starting with #, *, or %. */
				while( (chCard[0] == '#' || chCard[0] == '*') || chCard[0] == '%' )
				{
					if( lgEOF )
					{
						fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
							abund.npSolar, i );
						puts( "[Stop in ParseAbundances]" );
						cdEXIT(EXIT_FAILURE);
					}
					/* read in another line */
					input_readarray(chCard,&lgEOF);
				}

				/* we have the line image, print it */
				if( called.lgTalk )
				{
					fprintf( ioQQQ, "                       * ");
					k=0;
					while( chCard[k]!='\0' )
					{
						fprintf(ioQQQ,"%c",chCard[k]);
						++k;
					}
					while( k<80 )
					{
						fprintf(ioQQQ,"%c",' ');
						++k;
					}
					fprintf( ioQQQ,"*\n"); 
				}

				/* now convert to caps */
				caps(chCard);
				if( strncmp(chCard,"CONT",4) != 0 )
				{
					fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
						abund.npSolar, i );
					puts( "[Stop in ParseAbundances]" );
					cdEXIT(EXIT_FAILURE);
				}
				else
				{
					j = 1;
					absav[i] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
					if( lgEOF )
					{
						fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
							abund.npSolar, i);
						puts( "[Stop in ParseAbundances]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
			}
		}

		/* fell through to here after reading in abundances
			* check that there are no more - could mean typo further upstream */
		chk = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);

		if( !lgEOL || (chk!=0.) )
		{
			/* got another number, not lgEOL, so too many numbers */
			fprintf( ioQQQ, " There were more than %3ld 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] = (float)pow(10.,absav[i]);
			}
		}
		else
		{
			/* scale factors relative to solar */
			for( i=0; i < abund.npSolar; i++ )
			{
				abund.solar[abund.ipSolar[i]-1] *= (float)absav[i];
			}
		}

		/* check whether the abundances are reasonable */
		if( abund.solar[1] > 0.2 )
		{
			fprintf( ioQQQ, " Is an abundance of%10.3e for %2.2s reasonable?\n", 
				abund.solar[1], elementnames.chElementSym[1] );
		}

		for( i=2; i < LIMELM; i++ )
		{
			if( abund.solar[i] > 0.1 )
			{
				fprintf( ioQQQ, " Is an abundance of%10.3e for %2.2s reasonable?\n", 
					abund.solar[i], elementnames.chElementSym[i] );
			}
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->ParseAbundances()\n", debug_fp );
#		endif
		return;
	}
	/* no numbers on the line if it got this far */
	/* enter std abundances for various objects */
	if( lgMatch(" AGB",chCard) || lgMatch("AGB ",chCard) || lgMatch("PLAN",chCard) )
	{
		/* AGB/planetary nebula abundances */
		/* only turn on grains if "no grains" is not present */
		if( !lgMatch("O GR",chCard) )
		{
			/* turn on grains if not already done */
			if( !lgDSet )
			{
				/* return bins allocated by previous abundances ... commands */
				ReturnGrainBins();
				/* now allocate new grain bins */
				gp.dep = 1.;
				gp.lgAbunVsDepth = FALSE;
				gp.lgRequestQHeating = TRUE;
				gp.lgForbidQHeating = FALSE;

				/* >>chng 03 may 21, add this option to abund command */
				/* the keyword "no qheat" always takes precedence */
				if( lgMatch("O QH",chCard) ) 
				{
					gp.lgForbidQHeating = TRUE;
					gp.lgRequestQHeating = FALSE;
					phycon.lgPhysOK = FALSE;
				}

				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];
		}
	}

	else if( lgMatch("HII ",chCard) || lgMatch("H II",chCard) || lgMatch("ORIO",chCard) )
	{
		/* H II region abundances */
		/* only turn on grains if "no grains" is not present */
		if( !lgMatch("O GR",chCard) )
		{
			/* option to turn on grains */
			if( !lgDSet )
			{
				/* return bins allocated by previous abundances ... commands */
				ReturnGrainBins();
				/* now allocate new grain bins */
				gp.dep = 1.;
				gp.lgAbunVsDepth = FALSE;
				gp.lgRequestQHeating = TRUE;
				gp.lgForbidQHeating = FALSE;

				/* >>chng 03 may 21, add this option to abund command */
				/* the keyword "no qheat" always takes precedence */
				if( lgMatch("O QH",chCard) ) 
				{
					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];
		}
	}

	else if( lgMatch("ISM ",chCard) || lgMatch(" ISM",chCard) )
	{
		/* ISM abundances from Cowie and Songaila Ann Rev '86 */
		/* only turn on grains if "no grains" is not present */
		if( !lgMatch("O GR",chCard) )
		{
			if( !lgDSet )
			{
				/* return bins allocated by previous abundances ... commands */
				ReturnGrainBins();
				/* now allocate new grain bins */
				gp.dep = 1.;
				gp.lgAbunVsDepth = FALSE;
				gp.lgRequestQHeating = TRUE;
				gp.lgForbidQHeating = FALSE;

				/* >>chng 03 may 21, add this option to abund command */
				/* the keyword "no qheat" always takes precedence */
				if( lgMatch("O QH",chCard) ) 
				{
					gp.lgForbidQHeating = TRUE;
					gp.lgRequestQHeating = FALSE;
					phycon.lgPhysOK = FALSE;
				}

				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];
		}
	}

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

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

		/* now turn off the heavy elements */
		for( i=4; i < LIMELM; i++ )
		{
			/* >>chng 01 dec 17, from setting off to call to ParseElement,
			 * so that H-like and He-like level limits are handled properly */
			char chDUMMY[LIMELM];
			sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
			/* now convert to caps */
			caps(chDUMMY);
			ParseElement( chDUMMY );
			/*dense.lgElmtOn[i] = FALSE;*/
		}
	}

	else if( lgMatch("CAME",chCard) )
	{
		/* mix from Cameron 1982, "Essays on Nuclear Astrophysics" */
		for( i=0; i < LIMELM; i++ )
		{
			abund.solar[i] = abund.camern[i];
		}
	}

	else
	{
		fprintf( ioQQQ, 
			" ABUND must have PLAN, H II, CAMERON, NOVA, NLR, ALL, STARBURST, OLD SOLAR 84 or PRIMORDIAL.  Sorry.\n" );
		puts( "[Stop in ParseAbundances]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->ParseAbundances()\n", debug_fp );
#	endif
	return;
}

