/* 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 */
/*ParseElement parse options on element command */
#include "cddefines.h"
#include "optimize.h"
#include "abund.h"
#include "dense.h"
#include "input.h"
#include "iso.h"
#include "hmi.h"
#include "mole.h"
#include "elementnames.h"
#include "parser.h"

void ParseElement( Parser &p)
{
	/* this will be used to remember how many levels are active in any element we turn off,
	 * so that we can retain state if turned back on  */
	static bool lgFirst = true;
	static long levels[NISO][LIMELM];
	bool lgEnd, 
	  lgHIT;

	bool lgForceLog=false, lgForceLinear=false;

	long int nelem, 
	  j, 
	  k;
	double param=0.0;

	DEBUG_ENTRY( "ParseElement()" );

	/* zero out this array if first call */
	if( lgFirst )
	{
		lgFirst = false;
		for( long i=0; i<NISO; ++i )
		{
			for( nelem=0; nelem<LIMELM; ++nelem )
			{
				levels[i][nelem] = 0;
			}
		}
	}
	/* say that abundances have been changed */
	abund.lgAbnSolar = false;

	/* read in list of elements for abundances command */
	if( p.nMatch("READ") )
	{
		abund.npSolar = 0;
		for( long i=1; i < LIMELM; i++ )
		{
			p.getline();

			if( p.m_lgEOF )
			{
				fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" );
				cdEXIT(EXIT_FAILURE);
			}

			/* this would be a line starting with END to say end on list */
			if( p.strcmp("END") == 0 )
			{
				return;
			}

			j = 1;
			lgHIT = false;
			while( j < LIMELM && !lgHIT )
			{
				j += 1;

				if( p.strcmp( elementnames.chElementNameShort[j-1] ) == 0 )
				{
					abund.npSolar += 1;
					abund.ipSolar[abund.npSolar-1] = j;
					lgHIT = true;
				}
			}

			if( !lgHIT )
			{
				p.PrintLine( ioQQQ);
				fprintf( ioQQQ, " Sorry, but I did not recognize element name on this line.\n" );
				fprintf( ioQQQ, " Here is the list of names I recognize.\n" );
				fprintf( ioQQQ, " " );

				for( k=2; k <= LIMELM; k++ )
				{
					fprintf( ioQQQ, "%4.4s\n", elementnames.chElementNameShort[k-1] );
				}

				cdEXIT(EXIT_FAILURE);
			}
		}

		/* fell through, make sure one more line, the end line */
		p.getline();

		if( p.strcmp("END") == 0 )
		{
			return;
		}

		else
		{
			fprintf( ioQQQ, " Too many elements were entered.\n" );
			fprintf( ioQQQ, " I only know about%3d elements.\n", 
			  LIMELM );
			fprintf( ioQQQ, " Sorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* find which element - will be used in remainder of routine to 
	 * adjust aspects of this element */
	nelem = p.GetElem( );

	bool lgElementSet = false;
	if( nelem >= ipHYDROGEN )
	{
		lgElementSet = true;
		/* nelem is now the atomic number of the element, and must be correct */
		ASSERT( nelem>=0 && nelem < LIMELM );
	}

	/* look for log or linear keywords */
	if( p.nMatch(" LOG") )
		lgForceLog = true;
	else if( p.nMatch("LINE") )
		lgForceLinear = true;

	if( p.nMatch("SCAL") )
	{
		param = p.FFmtRead();

		/* enter abundance as scale factor, relative to what is in now */
		if( p.lgEOL() )
		{
			fprintf( ioQQQ, " There must be a number on this line.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}

		else
		{
			if( !lgElementSet )
			{
				fprintf( ioQQQ, 
					" ParseElement did not find an element on the following line:\n" );
				p.PrintLine( ioQQQ );
				cdEXIT(EXIT_FAILURE);
			}
			/* interpret as log unless forced linear */
			if( (lgForceLog || param <= 0.) && !lgForceLinear )
			{
				/* option for log of scale factor */
				param = pow(10.,param);
			}
			abund.ScaleElement[nelem] = (realnum)param;
		}
	}

	else if( p.nMatch("ABUN") )
	{
		param = p.FFmtRead();

		/* log of actual abundance */
		if( p.lgEOL() )
		{
			fprintf( ioQQQ, " There must be a number on this line.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}

		else
		{
			if( !lgElementSet )
			{
				fprintf( ioQQQ, 
					" ParseElement did not find an element on the following line:\n" );
				p.PrintLine( ioQQQ );
				cdEXIT(EXIT_FAILURE);
			}
			/* interpret as log unless forced linear */
			if( !lgForceLinear )
			{
				/* option for log of scale factor */
				param = pow(10.,param);
			}
			abund.solar[nelem] = (realnum)param;

			if( abund.solar[nelem] > 1. )
			{
				fprintf( ioQQQ, 
					" Please check the abundance of this element.  It seems high to me.\n" );
			}
		}
	}

	else if( p.nMatch(" OFF") )
	{
		param = p.FFmtRead();
	
		/* keyword LIMIT sets lower limit to abundance of element
		 * that will be included */
		/* option to turn off this element, set abundance to zero */
		if( p.nMatch( "LIMI") )
		{
			if( p.lgEOL() )
			{
				fprintf( ioQQQ, " There must be an abundances on the ELEMENT OFF LIMIT command.\n" );
				fprintf( ioQQQ, " Sorry.\n" );
				cdEXIT(EXIT_FAILURE);
			}
			dense.AbundanceLimit = (realnum)pow(10., param);
		}
		else
		{
			if( !lgElementSet )
			{
				fprintf( ioQQQ, 
					" ParseElement did not find an element on the following line:\n" );
				p.PrintLine ( ioQQQ );
				cdEXIT(EXIT_FAILURE);
			}
			/* specify element name */
			dense.lgElmtOn[nelem] = false;
			/* another flag that element is off */
			dense.IonHigh[nelem] = -1;
			dense.lgSetIoniz[nelem] = false;
			/* set no levels for all elements that are turned off, except for helium itself, which always exists */
			if( nelem > ipHELIUM )
			{
				levels[ipHYDROGEN][nelem] = iso_sp[ipH_LIKE][nelem].numLevels_max;
				levels[ipHELIUM][nelem] = iso_sp[ipHE_LIKE][nelem].numLevels_max;
				iso_sp[ipH_LIKE][nelem].numLevels_max = 0;
				iso_sp[ipHE_LIKE][nelem].numLevels_max = 0;
			}

			if( nelem == ipHYDROGEN )
			{
				fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" );
				fprintf( ioQQQ, " Sorry.\n" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	/* specify an ionization distribution */
	else if( p.nMatch("IONI") )
	{
		bool lgLogSet = false;
		int ion;
		int ihi , low;


		if( !lgElementSet )
		{
			fprintf( ioQQQ, 
				" ParseElement did not find an element on the following line:\n" );
			p.PrintLine( ioQQQ );
			cdEXIT(EXIT_FAILURE);
		}
		if( !dense.lgElmtOn[nelem] )
		{
			fprintf(ioQQQ,"Sorry, you cannot set the ionization of %s since it has been turned off.\n" ,
				elementnames.chElementName[nelem] );
			cdEXIT(EXIT_FAILURE);
		}

		/* option to specify ionization distribution for this element */
		dense.lgSetIoniz[nelem] = true;
		ion = 0;
		while( ion<nelem+2 )
		{
			/* the ionization fractions that are set when above set true,
			 * gas phase abundance is this times total abundance
			 * Ionization fraction for [nelem][ion] */
			dense.SetIoniz[nelem][ion] = (realnum)p.FFmtRead();

			if( p.lgEOL() ) 
				break;

			/* all are log if any are negative */
			if( dense.SetIoniz[nelem][ion] < 0. )
				lgLogSet = true;
			++ion;
		}
		param = dense.SetIoniz[nelem][0];

		/* zero out ones not specified */
		for( long i=ion; i<nelem+2; ++i )
			dense.SetIoniz[nelem][i] = 0.;

		/* convert rest to linear if any were negative */
		if( lgLogSet )
		{
			for( long i=0; i<ion; ++i )
				dense.SetIoniz[nelem][i] = (realnum)pow((realnum)10.f, dense.SetIoniz[nelem][i]);
		}

		/* now check that no zero abundances exist between lowest and highest non-zero values */
		low = 0;
		while( dense.SetIoniz[nelem][low]==0 && low < nelem+1 )
			++low;
		ihi = nelem+1;
		while( dense.SetIoniz[nelem][ihi]==0 && ihi > low )
			--ihi;

		if( ihi==low && dense.SetIoniz[nelem][ihi]==0 )
		{
			fprintf(ioQQQ," element ionization command has all zero ionization fractions.  This is not possible.\n Sorry\n");
			cdEXIT(EXIT_FAILURE);
		}
		realnum AbundSum=0.;
		for(long i=low; i<=ihi; ++i )
		{
			if( dense.SetIoniz[nelem][i]==0 )
			{
				fprintf(ioQQQ," element abundance command has zero abundance between positive values.  This is not possible.\n Sorry\n");
				cdEXIT(EXIT_FAILURE);
			}
			AbundSum += dense.SetIoniz[nelem][i];
		}

		// renormalize so that sum is unity - needed for various sum rules
		for(long i=low; i<=ihi; ++i )
			dense.SetIoniz[nelem][i] /= AbundSum;

	}

	/* option to turn element on - some ini files turn off elements,
	 * this can turn them back on */
	else if( p.nMatch(" ON ") )
	{

		if( !lgElementSet )
		{
			fprintf( ioQQQ, 
				" ParseElement did not find an element on the following line:\n" );
			p.PrintLine( ioQQQ );
			cdEXIT(EXIT_FAILURE);
		}
		/* option to turn off this element, set abundance to zero */
		dense.lgElmtOn[nelem] = true;
		/* reset levels to default if they were ever turned off with element off command */
		if( levels[ipHYDROGEN][nelem] )
		{
			iso_sp[ipH_LIKE][nelem].numLevels_max = levels[ipHYDROGEN][nelem];
			iso_sp[ipHE_LIKE][nelem].numLevels_max = levels[ipHELIUM][nelem];
		}
	}

	else if( p.nMatch("TABL") )
	{

		if( !lgElementSet )
		{
			fprintf( ioQQQ, 
				" ParseElement did not find an element on the following line:\n" );
			p.PrintLine( ioQQQ );
			cdEXIT(EXIT_FAILURE);
		}
		/* read in table of position-dependent abundances for a particular element. */
		abund.lgAbunTabl[nelem] = true;

		/* general flag saying this option turned on */
		abund.lgAbTaON = true;

		/* does the table give depth or radius?  keyword DEPTH
		 * says it is depth, default is radius */
		if( p.nMatch("DEPT") )
			abund.lgAbTaDepth[nelem] = true;
		else
			abund.lgAbTaDepth[nelem] = false;

		/* make sure not trying to change hydrogen */
		if( nelem == ipHYDROGEN )
		{
			fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}

		/* read pair giving radius and abundance */
		p.getline();	  
		abund.AbTabRad[0][nelem] = (realnum)p.FFmtRead();
		abund.AbTabFac[0][nelem] = (realnum)p.FFmtRead();

		if( p.lgEOL() )
		{
			fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" );
			cdEXIT(EXIT_FAILURE);
		}

		/* number of points in the table */
		abund.nAbunTabl = 2;

		lgEnd = false;
		/* LIMTAB is limit to number of points we can store and is
		 * set to 500 in abundances */
		while( !lgEnd && abund.nAbunTabl < LIMTABD )
		{
			p.getline();
			lgEnd = p.m_lgEOF;
			if( !lgEnd )
			{
				/* convert first 4 char to caps, into chCap */
				if( p.strcmp("END") == 0 )
					lgEnd = true;
			}

			/* lgEnd may have been set within above if, if end line encountered*/
			if( !lgEnd )
			{
				abund.AbTabRad[abund.nAbunTabl-1][nelem] = 
					(realnum)p.FFmtRead();

				abund.AbTabFac[abund.nAbunTabl-1][nelem] = 
					(realnum)p.FFmtRead();
				abund.nAbunTabl += 1;
			}
		}

		abund.nAbunTabl -= 1;

		/* now check that radii are in increasing order */
		for( long i=1; i < abund.nAbunTabl; i++ )
		{
			/* the radius values are assumed to be strictly increasing */
			if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] )
			{
				fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii "
					"must be in increasing order\n" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	else
	{
		fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be "
			"a keyword on this line.\n" );
		fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, "
			"_ON_, IONIZATION, and ABUNDANCE.\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		cdEXIT(EXIT_FAILURE);
	}

	/* vary option */
	if( optimize.lgVarOn )
	{
		if( p.nMatch("SCAL") )
		{
			/* vary scale factor */
			sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG", 
			  elementnames.chElementNameShort[nelem] );

		}
		else if( p.nMatch("ABUN") )
		{
			/* vary absolute abundance */
			sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUNDANCE %%f LOG", 
			  elementnames.chElementNameShort[nelem] );
		}
		/* pointer to where to write */
		optimize.nvfpnt[optimize.nparm] = input.nRead;
		if( param >=0. )
			optimize.vparm[0][optimize.nparm] = (realnum)log10(param);
		else
			optimize.vparm[0][optimize.nparm] = (realnum)param;
		optimize.vincr[optimize.nparm] = 0.2f;
		optimize.nvarxt[optimize.nparm] = 1;
		++optimize.nparm;
	}
	return;
}