/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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 "parse.h"
/*lint -e661 creation of  out of bound pointer */
/*lint -e817 Conceivably negative subscript */

void ParseElement(char *chCard )
{
	/* 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 int lgFirst = TRUE;
	static long levels[NISO][LIMELM];
	char chCap[INPUT_LINE_LENGTH];
	int lgEOL, 
	  lgEnd, 
	  lgHIT;

	int lgForceLog=FALSE, lgForceLinear=FALSE;

	long int i, 
	  nelem, 
	  j, 
	  k;
	double param;

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

	/* zero out this array if first call */
	if( lgFirst )
	{
		lgFirst = FALSE;
		for( 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( lgMatch("READ",chCard) )
	{
		abund.npSolar = 0;
		for( i=1; i < LIMELM; i++ )
		{
			input_readarray(chCard,&lgEnd);

			if( lgEnd )
			{
				fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" );
				puts( "[Stop in ParseElement]" );
				cdEXIT(EXIT_FAILURE);
			}

			strcpy( chCap, chCard );
			caps(chCap);

			/* this would be a line starting with END to say end on list */
			if( strncmp(chCap,"END",3) == 0 )
			{
				
#				ifdef DEBUG_FUN
				fputs( " <->ParseElement()\n", debug_fp );
#				endif
				return;
			}

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

				if( strncmp( chCap,elementnames.chElementNameShort[j-1] , 4) == 0 )
				{
					abund.npSolar += 1;
					abund.ipSolar[abund.npSolar-1] = j;
					lgHIT = TRUE;
				}
			}

			if( !lgHIT )
			{
				fprintf( ioQQQ, "%80.80s\n", chCard );
				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] );
				}

				puts( "[Stop in ParseElement]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* fell through, make sure one more line, the end line */
		input_readarray(chCard,&lgEnd);
		strcpy( chCap, chCard );
		caps(chCap);

		if( strncmp(chCap,"END",3) == 0 )
		{
			
#			ifdef DEBUG_FUN
			fputs( " <->ParseElement()\n", debug_fp );
#			endif
			return;
		}

		else
		{
			fprintf( ioQQQ, " Too many elements were entered.\n" );
			fprintf( ioQQQ, " I only know about%3ld elements.\n", 
			  LIMELM );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseElement]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* find which element */
	i = 0;

	/* will be used in remainder of routine to point to particular element */
	nelem = GetElem(chCard );

	if( nelem < 0 )
	{
		fprintf( ioQQQ, 
			" ParseElement did not find an element on the following line:\n" );
		fprintf( ioQQQ, 
			"%80.80s\n", chCard );
		puts( "[Stop in ParseElement]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* nelem is now the atomic number of the element, and must be correct */
	ASSERT( nelem>=0 && nelem < LIMELM );

	/* look for log or linear keywords */
	/* >>chng 01 jul 23, add this option to force linear (log had been tehre) */
	if( lgMatch(" LOG",chCard) )
	{
		lgForceLog = TRUE;
	}
	else if( lgMatch("LINE",chCard) )
	{
		lgForceLinear = TRUE;
	}

	i = 4;
	param = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

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

		else
		{
			/* interpret as log unless forced linear */
			if( (lgForceLog || param <= 0.) && !lgForceLinear )
			{
				/* option for log of scale factor */
				param = pow(10.,param);
			}
			abund.ScaleElement[nelem] = (float)param;
		}
	}

	else if( lgMatch("ABUN",chCard) )
	{
		/* log of actual abundance */
		if( lgEOL )
		{
			fprintf( ioQQQ, " There must be a number on this line.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseElement]" );
			cdEXIT(EXIT_FAILURE);
		}

		else
		{
			if( lgForceLinear )
			{
				abund.solar[nelem] = (float)param;
			}
			else
			{
				abund.solar[nelem] = (float)pow(10.,param);
			}

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

	else if( lgMatch(" OFF",chCard) )
	{
		/* option to turn off this element, set abundance to zero */
		dense.lgElmtOn[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.numLevels_max[ipH_LIKE][nelem];
			levels[ipHELIUM][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];
			iso.numLevels_max[ipH_LIKE][nelem] = 0;
			iso.numLevels_max[ipHE_LIKE][nelem] = 0;
		}

		if( nelem == 0 )
		{
			fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseElement]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* >>chng 01 may 08, add ionization option to specify an ionization distribution */
	else if( lgMatch("IONI",chCard) )
	{
		int lgLogSet = FALSE;
		int ion;
		int ihi , low;

		i = 4;
		/* 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 aabundance
			 * Ionization fraction for [nelem][ion] */
			dense.SetIoniz[nelem][ion] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

			if( lgEOL ) 
				break;

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

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

		/* convert rest to linear if any were negative */
		if( lgLogSet )
		{
			for( i=0; i<ion; ++i )
			{
				dense.SetIoniz[nelem][i] = (float)pow(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 abundance command has all zero abundances.  This is not possible.\n Sorry\n");
			puts( "[Stop in ParseElement]" );
			cdEXIT(EXIT_FAILURE);
		}
		for(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");
				puts( "[Stop in ParseElement]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
		/* >>chng 05 mar 13, add this
		 * finally turn off any part of the chemical networks that might use this element,
		 * since chemistry must come into equilibrium with the ionization, and vice vese -
		 * ion distribution cannot come into equilbrium with chemistry when ions are set,
		 * so much turn off chemistry */
		if( nelem==ipHYDROGEN )
		{
			fprintf(ioQQQ," Hydrogen ionization has been set, H chemistry is disabled.\n");
			/* turn off only H2 */
			hmi.lgNoH2Mole = TRUE;
		}
		/* loop over all atoms in co chem net - these occur as the last NUM_ELEMENTS 
		 * elements in the full set of NUM_HEAVY_MOLEC */
		for( i=NUM_HEAVY_MOLEC+NUM_ELEMENTS; i<NUM_COMOLE_CALC; ++i )
		{
			if( nelem==co.nelem_hevmol[i] )
			{
				fprintf(ioQQQ," The ionization of an element in the CO chemistry network has been set, CO chemistry is disabled.\n");
				/* turn off CO network */
				co.lgNoCOMole = TRUE;
			}
		}
	}

	/* >>chng 01 mar 29, add on option */
	else if( lgMatch(" ON ",chCard) )
	{
		/* 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.numLevels_max[ipH_LIKE][nelem] = levels[ipHYDROGEN][nelem];
			iso.numLevels_max[ipHE_LIKE][nelem] = levels[ipHELIUM][nelem];
		}
	}

	else if( lgMatch("TABL",chCard) )
	{
		/* >>chng 97 jun 02, add table option
		 * when called, read in densities from input stream */
		abund.lgAbunTabl[nelem] = TRUE;

		/* general flag saying this option turned on */
		abund.lgAbTaON = TRUE;
		if( lgMatch("DEPT",chCard) )
		{
			abund.lgAbTaDepth[nelem] = TRUE;
		}
		else
		{
			abund.lgAbTaDepth[nelem] = FALSE;
		}

		/* make sure not trying to change hydrogen */
		if( nelem == 0 )
		{
			fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseElement]" );
			cdEXIT(EXIT_FAILURE);
		}

		input_readarray(chCard,&lgEnd);
		i = 1;
		abund.AbTabRad[0][nelem] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		abund.AbTabFac[0][nelem] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		if( lgEOL )
		{
			fprintf( ioQQQ, " no pairs entered - cant interpolate\n" );
			puts( "[Stop in ParseElement]" );
			cdEXIT(EXIT_FAILURE);
		}

		abund.nAbunTabl = 2;
		lgEnd = FALSE;
		/* LIMTAB is set to 500 in abundances */
		while( !lgEnd && abund.nAbunTabl < LIMTABD )
		{
			input_readarray(chCard,&lgEnd);
			if( !lgEnd )
			{
				/* convert first 4 char to caps, into chCap */
				cap4(chCap , chCard);
				if( strncmp(chCap,"END",3) == 0 )
					lgEnd = TRUE;
			}

			/* lgEnd may have been set within above if, if end line encountered*/
			if( !lgEnd )
			{
				i = 1;
				abund.AbTabRad[abund.nAbunTabl-1][nelem] = 
					(float)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);

				abund.AbTabFac[abund.nAbunTabl-1][nelem] = 
					(float)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
				abund.nAbunTabl += 1;
			}
		}

		abund.nAbunTabl -= 1;

		/* now chec that abundances are in increasing order */
		for( 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, " abun radii must be in increasing order\n" );
				puts( "[Stop in ParseElement]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	else
	{
		fprintf( ioQQQ, " 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" );
		puts( "[Stop in ParseElement]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* vary option */
	if( optimize.lgVarOn )
	{
		optimize.nvarxt[optimize.nparm] = 1;
		/* pointer to where to write */
		optimize.nvfpnt[optimize.nparm] = input.nRead;

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

			/* param is linear scale factor */
			optimize.vparm[0][optimize.nparm] = (float)log10(param);
			optimize.vincr[optimize.nparm] = 0.2f;
		}

		else if( lgMatch("ABUN",chCard) )
		{
			/* vary absolute abundance */
			sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUND %%f LOG ", 
			  elementnames.chElementNameShort[nelem] );

			/* param is log of abundance by number relative to hydrogen */
			optimize.vparm[0][optimize.nparm] = (float)param;
			optimize.vincr[optimize.nparm] = 0.2f;
		}
		++optimize.nparm;
	}

#	ifdef DEBUG_FUN
	fputs( " <->ParseElement()\n", debug_fp );
#	endif
	return;
}
/*lint +e661 creation of  out of bound pointer */
/*lint +e817 Conceivably negative subscript */

