/*strbst generate abundance set from Fred Hamann's starburst evolution grid */
#include <string.h>
#include "cddefines.h"
#include "solar.h"
#include "solars.h"
#include "varypar.h"
#include "input.h"
#include "strbst.h"
#include "ffmtread.h"
#include "lgmatch.h"
#include "elementnames.h"

void strbst(char *chCard)
{
	int lgDebug, 
	  lgEOL;
	long int i, 
	  j;
	double sqrzed, 
	  zHigh, 
	  zal, 
	  zar, 
	  zc, 
	  zca, 
	  zcl, 
	  zco, 
	  zed, 
	  zed2, 
	  zed3, 
	  zed4, 
	  zedlog, 
	  zfe, 
	  zh, 
	  zhe, 
	  zmg, 
	  zn, 
	  zna, 
	  zne, 
	  zni, 
	  zo, 
	  zs, 
	  zsi;
	/* this is largest stored metallicity */
	static double zLimit = 35.5;

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


	if( lgMatch("TRAC",chCard) )
	{
		lgDebug = TRUE;
	}
	else
	{
		lgDebug = FALSE;
	}

	i = 5;
	/* argument is metallicity  */
	zed = FFmtRead(chCard,&i,76,&lgEOL);
	if( lgDebug )
	{
		/* trace keyword did appear
		 * this will not be used, but must trick the sanity test */
		zHigh = zLimit;

		/* if trace option set, print header now, and init ZED */
		fprintf( ioQQQ, " Abundances relative to H, Z\n" );

		/* this is actual header line */
		fprintf( ioQQQ, "     Z   ");

		/* make line of chemical symbol names */
		for(j=0; j < 30; j++)
		{
			fprintf( ioQQQ, "    %2.2s", elementnames.chElementSym[j]);
		}
		fprintf( ioQQQ, "    \n" );

		zed = 1e-3;
	}
	else if( lgEOL && !lgDebug )
	{
		/* no number and trace keyword did not appear */
		fprintf( ioQQQ, " The metallicity must appear on this line.\n" );
		puts( "[Stop in strbst]" );
		exit(1);
	}
	else
	{
		zHigh = zed;
	}

	/* interpolate on Fred Hamann's set of starburst abundances
	 * scan off keys _log or linear */
	if( lgMatch(" LOG",chCard) )
	{
		zed = pow(10.,zed);
		zHigh = zed;
	}

	else if( lgMatch("LINE",chCard) )
	{
		if( zed <= 0. )
		{
			fprintf( ioQQQ, " Z .le.0 not allowed, Z=%10.2e\n", 
			  zed );
			puts( "[Stop in strbst]" );
			exit(1);
		}
	}

	else
	{
		/* log, linear not specified, neg so log */
		if( zed <= 0. )
		{
			zed = pow(10.,zed);
			zHigh = zed;
		}
	}

	/* following if loop only if trace option is set
	 * >>chng 97 june 17, get rid of go to
	 *999  continue
	 * */
	while( zed <= zHigh )
	{
		if( zed < 1e-3 || zed > zLimit )
		{
			fprintf( ioQQQ, " The metallicity must be between 1E-3 and%10.2e\n", 
			  zLimit );
			puts( "[Stop in strbst]" );
			exit(1);
		}
		zed2 = zed*zed;
		zed3 = zed2*zed;
		zed4 = zed3*zed;
		zedlog = log(zed);
		sqrzed = sqrt(zed);
		/* The value of each element as a function of ZED=[Z] */
		zh = 1.081646723 - 0.04600492*zed + 8.6569e-6*zed2 + 1.922e-5*
		  zed3 - 2.0087e-7*zed4;

		/* helium */
		zhe = 0.864675891 + 0.044423807*zed + 7.10886e-5*zed2 - 5.3242e-5*
		  zed3 + 5.70194e-7*zed4;
		solarCom.solar[1] = (float)zhe;

		/* li, b, bo unchanged */
		solarCom.solar[2] = 1.;
		solarCom.solar[3] = 1.;
		solarCom.solar[4] = 1.;

		/* carbon */
		zc = 0.347489799 + 0.016054107*zed - 0.00185855*zed2 + 5.43015e-5*
		  zed3 - 5.3123e-7*zed4;
		solarCom.solar[5] = (float)zc;

		/* nitrogen */
		zn = -0.06549567 + 0.332073984*zed + 0.009146066*zed2 - 0.00054441*
		  zed3 + 6.16363e-6*zed4;
		if( zn < 0.0 )
		{
			zn = 0.05193*zed;
		}
		if( zed < 0.3 )
		{
			zn = -0.00044731103 + 0.00026453554*zed + 0.52354843*zed2 - 
			  0.41156655*zed3 + 0.1290967*zed4;
			if( zn < 0.0 )
			{
				zn = 0.000344828*zed;
			}
		}
		solarCom.solar[6] = (float)zn;

		/* oxygen */
		zo = 1.450212747 - 0.05379041*zed + 0.000498919*zed2 + 1.09646e-5*
		  zed3 - 1.8147e-7*zed4;
		solarCom.solar[7] = (float)zo;

		/* neon */
		zne = 1.110139023 + 0.002551998*zed - 2.09516e-7*zed3 - 0.00798157*
		  POW2(zedlog);
		solarCom.solar[9] = (float)zne;

		/* fluorine, scale from neon */
		solarCom.solar[8] = solarCom.solar[9];

		/* sodium */
		zna = 1.072721387 - 0.02391599*POW2(zedlog) + .068644972*
		  zedlog + 0.017030935/sqrzed;
		/* this one is slightly negative at very low Z */
		zna = MAX2(1e-12,zna);
		solarCom.solar[10] = (float)zna;

		/* magnesium */
		zmg = 1.147209646 - 7.9491e-7*POW3(zed) - .00264458*POW2(zedlog) - 
		  0.00635552*zedlog;
		solarCom.solar[11] = (float)zmg;

		/* aluminium */
		zal = 1.068116358 - 0.00520227*sqrzed*zedlog - 0.01403851*
		  POW2(zedlog) + 0.066186787*zedlog;
		/* this one is slightly negative at very low Z */
		zal = MAX2(1e-12,zal);
		solarCom.solar[12] = (float)zal;

		/* silicon */
		zsi = 1.067372578 + 0.011818743*zed - 0.00107725*zed2 + 3.66056e-5*
		  zed3 - 3.556e-7*zed4;
		solarCom.solar[13] = (float)zsi;

		/* phosphorus scaled from silicon */
		solarCom.solar[14] = solarCom.solar[13];

		/* sulphur */
		zs = 1.12000;
		solarCom.solar[15] = (float)zs;

		/* chlorine */
		zcl = 1.10000;
		solarCom.solar[16] = (float)zcl;

		/* argon */
		zar = 1.091067724 + 2.51124e-6*zed3 - 0.0039589*sqrzed*zedlog + 
		  0.015686715*zedlog;
		solarCom.solar[17] = (float)zar;

		/* potassium scaled from silicon */
		solarCom.solar[18] = solarCom.solar[13];

		/* calcium */
		zca = 1.077553875 - 0.00888806*zed + 0.001479866*zed2 - 6.5689e-5*
		  zed3 + 1.16935e-6*zed4;
		solarCom.solar[19] = (float)zca;

		/* iron */
		zfe = 0.223713045 + 0.001400746*zed + 0.000624734*zed2 - 3.5629e-5*
		  zed3 + 8.13184e-7*zed4;
		solarCom.solar[25] = (float)zfe;

		/* scandium, scaled from iron */
		solarCom.solar[20] = solarCom.solar[25];

		/* titanium, scaled from iron */
		solarCom.solar[21] = solarCom.solar[25];

		/* vanadium scaled from iron */
		solarCom.solar[22] = solarCom.solar[25];

		/* chromium scaled from iron */
		solarCom.solar[23] = solarCom.solar[25];

		/* manganese scaled from iron */
		solarCom.solar[24] = solarCom.solar[25];

		/* cobalt */
		zco = zfe;
		solarCom.solar[26] = (float)zco;

		/* nickel */
		zni = zfe;
		solarCom.solar[27] = (float)zni;

		/* copper scaled from iron */
		solarCom.solar[28] = solarCom.solar[25];

		/* zinc  scaled from iron */
		solarCom.solar[29] = solarCom.solar[25];

		/* rescale to true abundances */
		solarCom.solar[0] = 1.;
		solarCom.solar[1] = (float)(solarCom.solar[1]*solars.SolarSave[1]/
		  zh);

		for( i=2; i < LIMELM; i++ )
		{
			solarCom.solar[i] = (float)(solarCom.solar[i]*solars.SolarSave[i]*
			  zed/zh);
		}

		if( lgDebug )
		{
			fprintf( ioQQQ, "%10.2e", zed );
			for( i=0; i < LIMELM; i++ )
			{
				fprintf( ioQQQ, "%6.2f", log10(solarCom.solar[i]) );
			}
			fprintf( ioQQQ, "\n" );

			if( zed == zLimit )
			{
				puts( "[Stop in strbst]" );
				exit(1);
			}
		}

		/* this trick makes sure last one is upper limit */
		if( zed < zLimit )
		{
			zed = MIN2(zed*1.5,zLimit);
		}
		else
		{
			zed *= 1.5;
		}
	}

	/* vary option */
	if( VaryPar.lgVarOn )
	{
		/* this is number of parameters to feed onto the input line */
		VaryPar.nvarxt[VaryPar.nparm] = 1;
		strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "ABUNDANCES STARBURST %9.5F LOG" );
		/* following are upper and lower limits to metallicity range */
		VaryPar.varang[VaryPar.nparm][0] = (float)log10(1e-3);
		VaryPar.varang[VaryPar.nparm][1] = (float)log10(zLimit);
		/* pointer to where to write */
		VaryPar.nvfpnt[VaryPar.nparm] = input.nRead;
		/* log of metallicity will be argument */
		VaryPar.vparm[0][VaryPar.nparm] = (float)log10(zed);
		VaryPar.vincr[VaryPar.nparm] = 0.2f;
		++VaryPar.nparm;
	}


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