/* 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 */
/*AbundancesPrt print all abundances, both gas phase and grains */
/*AbundancesSet sets initial abundances after parameters are entered by reading input */
/*AbundancesTable interpolate on table of points to do 'element table' command, */
/*PrtElem print chemical composition at start of calculation */
#include "cddefines.h"
#include "physconst.h"
#include "called.h"
#include "trace.h"
#include "wind.h"
#include "elementnames.h"
#include "dense.h"
#include "radius.h"
#include "grainvar.h"
#include "abund.h"

/*PrtElem print chemical composition at start of calculation */
static void PrtElem(
  /* the job to do, the options are "init", "fill", "flus" */
  char *chJob, 
  /* label for the element */
  char *chLabl, 
  /* its abundance */
  double abund_prt);

/*AbundancesPrt print all abundances, both gas phase and grains */
void AbundancesPrt( void )
{
	long int i,
		nd;
	double grain_num_rel_hydr[3];

#	ifdef DEBUG_FUN
	fputs( "<+>AbundancesSet()\n", debug_fp );
#	endif
	
	/* this is main loop to print abundances of each element */
	if( called.lgTalk )
	{
		PrtElem("initG","  ",0.);/* initialize print routine for gas*/
		for( i=0; i < LIMELM; i++ )
		{
			if( dense.lgElmtOn[i] )
			{
				/* fill in print buffer with abundances */
				PrtElem("fill",(char*)elementnames.chElementSym[i],
				  abund.solar[i]);
			}
		}

		/* flush the print buffer */
		PrtElem("flus","  ",0.);
		/* final carriage return */
		fprintf( ioQQQ, " \n" );

		/* now grains if present */
		if( gv.lgDustOn )
		{
			/* we will first print the total abundances of each element locked up in grains */
			/* initialize print routine for dust*/
			PrtElem("initD","  ",0.);
			for( i=0; i < LIMELM; i++ )
			{
				if( gv.elmSumAbund[i]>SMALLFLOAT )
				{
					/* fill in print buffer with abundances */
					PrtElem("fill",(char*)elementnames.chElementSym[i],
						gv.elmSumAbund[i]/dense.gas_phase[ipHYDROGEN]);
				}
			}
			/* flush the print buffer */
			PrtElem("flus","  ",0.);
			/* final carriage return */
			fprintf( ioQQQ, " \n" );

			/* this is used to store grain density per hydrogen */
			grain_num_rel_hydr[0] = 0.;
			grain_num_rel_hydr[1] = 0.;
			grain_num_rel_hydr[2] = 0.;
			for( nd=0; nd < gv.nBin; nd++ )
			{

				/* >>chng 03 Nov 21  set up abundances of pah by chemical composition and numbers relative
				* to Hydrogen.  This will then be printed in the .out file */
				if( (gv.bin[nd]->matType==MAT_CAR) || (gv.bin[nd]->matType ==MAT_CAR2 ) )
				{
					/* carbonaceous grains */
					grain_num_rel_hydr[0] = (gv.bin[nd]->IntVol/gv.bin[nd]->AvVol)*gv.bin[nd]->dstfactor;
				}
				else if( (gv.bin[nd]->matType ==MAT_SIL) || (gv.bin[nd]->matType ==MAT_SIL2))
				{
					/* silicate grains */
					grain_num_rel_hydr[1] = (gv.bin[nd]->IntVol/gv.bin[nd]->AvVol)*gv.bin[nd]->dstfactor;
				}

				else if( ( gv.bin[nd]->matType ==MAT_PAH) || (gv.bin[nd]->matType ==MAT_PAH2))
				{
					/* PAH grains */
					grain_num_rel_hydr[2] = (gv.bin[nd]->IntVol/gv.bin[nd]->AvVol)*gv.bin[nd]->dstfactor;
				}
				else
					TotalInsanity();

			}

			/* now print total number of grains of each type - there are 3 types */
			fprintf(ioQQQ,"                                                  Number of grains per hydrogen\n");
			fprintf(ioQQQ,"                                       Carbonaceous: %.3f  silicate: %.3f  PAH: %.3f\n\n" ,
			log10( MAX2(1e-30,grain_num_rel_hydr[0] ) ) ,
			log10( MAX2(1e-30,grain_num_rel_hydr[1] ) ) ,
			log10( MAX2(1e-30,grain_num_rel_hydr[2] ) ) );
		}
	}

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

}

/*AbundancesSet print all abundances, both gas phase and grains */
void AbundancesSet(void)
{
	long int i, 
	  nelem;
	double fac;
	static int lgFirstCall=TRUE;
	static int lgElOnOff[LIMELM];

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

	/* if this is the first call to this routine in this coreload, 
	 * save the state of the lgElmOn array, so that it is possible
	 * to turn off elements in later models, but not turn on an
	 * element that was initally turned off.  This is necessary since
	 * the Create... routines that create space for elements will
	 * not be revisited in later models.  You can turn off an initally
	 * enabled element, but not turn a disabled one on.  */

	if( lgFirstCall )
	{
		/* first call - save the inital state of the lgElmtOn vector */
		for( i=0; i<LIMELM; ++i )
		{
			lgElOnOff[i] = dense.lgElmtOn[i];
		}
	}
	lgFirstCall = FALSE;

	/* make sure that initially FALSE elements remain off, while letting 
	 * enabled elements be turned off */
	for( i=ipHYDROGEN; i<LIMELM; ++i )
	{
		dense.lgElmtOn[i] = lgElOnOff[i] && dense.lgElmtOn[i];
	}

	/* rescale so that abundances are H=1 */
	for( i=ipHELIUM; i < LIMELM; i++ )
	{
		abund.solar[i] /= abund.solar[0];
	}
	abund.solar[ipHYDROGEN] = 1.;

	/* set current abundances to "solar" times metals scale factor
	 * and grain depletion factor */
	abund.solar[ipHELIUM] *= abund.depset[1]*abund.ScaleElement[1];

	/* option for density or abundance variations, this flag is TRUE by default,
	 * set in zero, but set false if variations are enabled AND these
	 * are not density variations, but rather abundances */
	if( dense.lgDenFlucOn )
	{
		/* usual case - either density flucuations or none at all */
		fac = 1.;
	}
	else
	{
		/* abundance fluc enabled, set initial value */
		fac = dense.cfirst*cos(dense.flcPhase) + dense.csecnd;
	}

	for( i=ipLITHIUM; i < LIMELM; i++ )
	{
		abund.solar[i] *= (float)(abund.ScaleMetals*abund.depset[i]*
		  abund.ScaleElement[i]*fac);
	}

	/* now fix abundance of any element with element table set */
	if( abund.lgAbTaON )
	{
		for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
		{
			if( abund.lgAbunTabl[nelem] )
			{
				abund.solar[nelem] = (float)(AbundancesTable(radius.Radius,
				  radius.depth,nelem+1));
			}
		}
	}

	/* dense.gas_phase[nelem] contains total abundance of element */
	/* the density of hydrogen itself has already been set at this point -
	 * it is set when commands parsed, most likely by the hden command -
	 * set all heavier elements */
	for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			dense.gas_phase[nelem] = (float)(abund.solar[nelem]*dense.gas_phase[ipHYDROGEN]);
			if( dense.gas_phase[nelem] <= 0. )
			{
				fprintf( ioQQQ, " Abundances must be greater than zero.  Check entered abundance for element%3ld  = %2.2s\n", 
				nelem, elementnames.chElementSym[nelem] );
				puts( "[Stop in AbundancesSet]" );
				cdEXIT(EXIT_FAILURE);
			}
			else if( dense.gas_phase[nelem] < SMALLFLOAT )
			{
				fprintf(ioQQQ," Abundance for %s is %.2e, less than lower limit of %.3e, so turning element off.\n",
					elementnames.chElementSym[nelem],
					dense.gas_phase[nelem],
					SMALLFLOAT );
				dense.lgElmtOn[nelem] = FALSE;
			}
		}
		else
		{
			/* >>chng 04 apr 20, set to zero if element is off */
			dense.gas_phase[nelem] = 0.;
		}
	}

	/* usually computed in PRESUR, but not called before first zone */
	wind.agrav = (float)((6.67e-8*wind.comass)*(SOLAR_MASS/radius.Radius)/
	  radius.Radius);

	if( trace.lgTrace )
	{
		float sumx , sumy , sumz = 0.;

		sumx = dense.gas_phase[ipHYDROGEN]*dense.AtomicWeight[ipHYDROGEN];
		sumy = dense.gas_phase[ipHELIUM]*dense.AtomicWeight[ipHELIUM];

		fprintf( ioQQQ, "\n AbundancesSet sets following densities (cm^-3); \n" );
		for( i=0; i<3 ; i++ )
		{
			for( nelem=i*10; nelem < i*10+10; nelem++ )
			{
				fprintf( ioQQQ, " %2.2s", elementnames.chElementSym[nelem] );
				PrintE82( ioQQQ, dense.gas_phase[nelem] );
				if( nelem>ipHELIUM )
					sumz += dense.gas_phase[nelem]*dense.AtomicWeight[nelem] ;
			}
			fprintf( ioQQQ, " \n" );
		}
		fprintf( ioQQQ, "\n AbundancesSet sets following abundances rel to H; \n" );
		for( i=0; i<3 ; i++ )
		{
			for( nelem=i*10; nelem < i*10+10; nelem++ )
			{
				fprintf( ioQQQ, " %2.2s", elementnames.chElementSym[nelem] );
				PrintE82( ioQQQ, dense.gas_phase[nelem]/dense.gas_phase[ipHYDROGEN] );
			}
			fprintf( ioQQQ, " \n" );
		}
		fprintf( ioQQQ, " \n" );
		fprintf(ioQQQ," Gas-phase mass fractions, X:%.3e Y:%.3e Z:%.3e\n\n",
			sumx/SDIV(sumx+sumy+sumz) , 
			sumy/SDIV(sumx+sumy+sumz) , 
			sumz/SDIV(sumx+sumy+sumz) );
	}

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

/* this is number of elements across one line */
#define	NELEM1LINE	9

/*PrtElem print chemical composition at start of calculation */
static void PrtElem(
  /* the job to do, the options are "init", "fill", "flus" */
  char *chJob, 
  /* label for the element */
  char *chLabl, 
  /* its abundance */
  double abund_prt)
{
	static char chAllLabels[NELEM1LINE][14];/* buffer where elements will be stored*/
	long int i, 
	  noffset;
	static long int nelem;  /* counter for number of elements read in*/

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

	if( strcmp(chJob,"initG") == 0 )
	{
		/* gas phase abundances */
		nelem = 0;
		fprintf( ioQQQ, 
			"                                                  Gas Phase Chemical Composition\n" );
	}
	else if( strcmp(chJob,"initD") == 0 )
	{
		/* abundances in grains */
		nelem = 0;
		fprintf( ioQQQ, 
			"                                                    Grain Chemical Composition\n" );
	}

	else if( strcmp(chJob,"fill") == 0 )
	{
		/* print log of abundance to avoid exponential output */
		abund_prt = log10( abund_prt );
		/* stuff in labels and abundances */
		sprintf( chAllLabels[nelem], "  %2.2s:%8.4f", chLabl, abund_prt );
		if( nelem == NELEM1LINE-1 )
		{
			/* we hit as many as it will hold - print it out and reset*/
			fprintf( ioQQQ, "      " );
			for( i=0; i < NELEM1LINE; i++ )
			{
				fprintf( ioQQQ, "%13.13s", chAllLabels[i] );
			}
			fprintf( ioQQQ, "\n" );
			/* reset counter to zero */
			nelem = 0;
		}
		else
		{
			/* just increment */
			++nelem;
		}
	}

#	if 0
	/* Do this if you want to know about PAH number abundance */
	else if( strcmp(chJob,"fillp") == 0 )
	{
		/* print log of abundance to avoid exponential output */
		abund_prt = log10( abund_prt );

		/* stuff in labels and abundances */
		sprintf( chAllLabels[nelem], "  %2.2s:%8.4f", chLabl, abund_prt );
		if( nelem == NELEM1LINE-1 )
		{
			/* we hit as many as it will hold - print it out and reset*/
			fprintf( ioQQQ, "      " );
			for( i=0; i < NELEM1LINE; i++ )
			{
				fprintf( ioQQQ, "%13.13s", chAllLabels[i] );
			}
			fprintf( ioQQQ, "\n" );
			/* reset counter to zero */
			nelem = 0;
		}
		else
		{
			/* just increment */
			++nelem;
		}
	}
#	endif

	else if( strcmp(chJob,"flus") == 0 )
	{
		/* flush the stack */
		i = NELEM1LINE - (nelem - 2);
		noffset = i/2-1;
		/* make format pretty */
		fprintf( ioQQQ, "      " );

		for(i=0; i < noffset ; i++)
		{
			/* skip out this many fields */
			fprintf( ioQQQ, "             " );
		}

		/* if nelem is even we need to space out another 8 */
		if( !(nelem%2) && nelem > 0)
			fprintf( ioQQQ,"        ");

		for( i=0; i < nelem; i++ )
		{
			fprintf( ioQQQ, "%13.13s", chAllLabels[i] );
		}

		fprintf( ioQQQ, "\n" );
	}
	else
	{
		fprintf( ioQQQ, " PrtElem does not understand job=%4.4s\n", 
		  chJob );
		puts( "[Stop in prtelem]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	return;
}

	
/*AbundancesTable interpolate on table of points to do 'element table' command, */
double AbundancesTable(double r0, 
  double depth, 
  long int iel)
{
	int lgHit;
	long int j;
	double frac, 
	  tababun_v, 
	  x;

#	ifdef DEBUG_FUN
	fputs( "<+>AbundancesTable()\n", debug_fp );
#	endif
	/*interpolate on table of points to do 'element table' command,
	 *based on code by K Volk
	 *each line is log radius and abundance. */
	/* atomic number */

	/* interpolate on radius or depth? */
	if( abund.lgAbTaDepth[iel-1] )
	{
		/* depth key appeared = we want depth */
		x = log10(depth);
	}
	else
	{
		/* use radius */
		x = log10(r0);
	}

	/* this will be reset below, but is here as a safety check */
	tababun_v = -DBL_MAX;

	if( x < abund.AbTabRad[0][iel-1] || x >= abund.AbTabRad[abund.nAbunTabl-1][iel-1] )
	{
		fprintf( ioQQQ, " requested radius outside range of AbundancesTable\n" );
		fprintf( ioQQQ, " radius was%10.2e min, max=%10.2e%10.2e\n", 
		  x, abund.AbTabRad[0][iel-1], abund.AbTabRad[abund.nAbunTabl-1][iel-1] );
		puts( "[Stop in AbundancesTable]" );
		cdEXIT(EXIT_FAILURE);
	}

	else
	{
		lgHit = FALSE;
		j = 1;

		while( !lgHit && j <= abund.nAbunTabl - 1 )
		{
			if( abund.AbTabRad[j-1][iel-1] <= (float)x && 
				abund.AbTabRad[j][iel-1] > (float)x )
			{
				frac = (x - abund.AbTabRad[j-1][iel-1])/(abund.AbTabRad[j][iel-1] - 
				  abund.AbTabRad[j-1][iel-1]);
				tababun_v = abund.AbTabFac[j-1][iel-1] + frac*
				  (abund.AbTabFac[j][iel-1] - abund.AbTabFac[j-1][iel-1]);
				lgHit = TRUE;
			}
			j += 1;
		}

		if( !lgHit )
		{
			fprintf( ioQQQ, " radius outran dlaw table scale, requested=%6.2f largest=%6.2f\n", 
			  x, abund.AbTabRad[abund.nAbunTabl-1][iel-1] );
			puts( "[Stop in AbundancesTable]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* got it, now return value, not log of density */
	tababun_v = pow(10.,tababun_v);

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

