/*SetAbundances sets initial abundances after parameters are entered by reading input */
#include "cddefines.h"
#include "prtelem.h"
#include "physconst.h"
#include "abuntabllg.h"
#include "ionfracs.h"
#include "solar.h"
#include "scalem.h"
#include "deplon.h"
#include "dmetal.h"
#include "he.h"
#include "phycon.h"
#include "called.h"
#include "trace.h"
#include "wind.h"
#include "elmton.h"
#include "elementnames.h"
#include "radius.h"
#include "fluct.h"
#include "printe82.h"
#include "tababun.h"
#include "setabundances.h"

void SetAbundances(void)
{
	long int i, 
	  nelem;
	double fac;
	static int lgFirstCall=TRUE;
	static int lgElOnOff[LIMELM];

#	ifdef DEBUG_FUN
	fputs( "<+>SetAbundances()\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] = elmton.lgElmtOn[i];
		}
	}
	lgFirstCall = FALSE;

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

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

	/* set current abundances to "solar" times metals scale factor
	 * and grain depletion factor */
	solarCom.solar[1] *= deplon.depset[1]*scalem.ScaleElement[1];
	if( fluct.lgDenFluc )
	{
		/* usual case - either density flucuations or none at all */
		fac = 1.;
	}
	else
	{
		/* abundance fluc enabled, set initial value */
		fac = fluct.cfirst*cos(fluct.flcPhase) + fluct.csecnd;
	}

	for( i=2; i < LIMELM; i++ )
	{
		solarCom.solar[i] *= (float)(dmetal.ScaleMetals*deplon.depset[i]*
		  scalem.ScaleElement[i]*fac);
	}

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

	if( called.lgTalk )
	{
		PrtElem("init","  ",0.);/* initialize print routine*/
		for( i=0; i < LIMELM; i++ )
		{
			if( elmton.lgElmtOn[i] )
			{
				/* fill in print buffer with abundances */
				PrtElem("fill",(char*)elementnames.chElementSym[i],
				  solarCom.solar[i]);
			}
		}
		PrtElem("flus","  ",0.); /* flush the print buffer */
		fprintf( ioQQQ, " \n" );
	}

	/* xIonFracs(nelem,0) contains total abundance of element */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		xIonFracs[0][nelem] = (float)(solarCom.solar[nelem]*phycon.hden);
		if( xIonFracs[0][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 SetAbundances]" );
			exit(1);
		}
	}

	he.ahe = (float)(solarCom.solar[1]*phycon.hden);
	he.heiii = he.ahe;
	he.heii = 0.;
	he.hei = 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 )
	{
		fprintf( ioQQQ, " SetAbundances sets following densities; \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, xIonFracs[0][nelem] );
			}
			fprintf( ioQQQ, " \n" );
		}
	}


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