/* 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 */
/*ParseBlackbody parse parameters off black body command */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "radius.h"
#include "parse.h"

void ParseBlackbody(char *chCard,	/* input command line, already changed to caps */
  long int *nqh,				/* counter for which continuum source this is */
  float *ar1)					/* optional area that might be set here */
{
	int lgEOL;
	long int i;
	double a, 
	  dil, 
	  rlogl;

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

	/* type is blackbody */
	strcpy( rfield.chSpType[rfield.nspec], "BLACK" );

	/* these two are not used for this continuum shape */
	rfield.cutoff[rfield.nspec][0] = 0.;
	rfield.cutoff[rfield.nspec][1] = 0.;

	/* get the temperature */
	i = 5;
	rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* there must be at least one number, the temperature */
	if( lgEOL )
	{
		fprintf( ioQQQ, " There must be at least 1 number on a BLACKBODY command line.  Sorry.\n" );
		puts( "[Stop in ParseBlackbody]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this is the temperature - make sure its linear in the end
	 * there are two keys, LINEAR and LOG, that could be here,
	 * else choose which is here by which side of 10 */
	if( (rfield.slope[rfield.nspec] <= 10. && !lgMatch("LINE",chCard)) || 
		lgMatch(" LOG",chCard) )
	{
		/* log option */
		rfield.slope[rfield.nspec] = pow(10.,rfield.slope[rfield.nspec]);
	}

	/* check that temp is not too low - could happen if log misused */
	if( rfield.slope[rfield.nspec] < 1e4 )
	{
		fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n", 
		  rfield.slope[rfield.nspec] );
	}

	/* now check that temp not too low - would peak below low
	 * energy limit of the code
	 * factor is temperature of 1 Ryd, egamry is high energy limit of code */
	if( rfield.slope[rfield.nspec]/TE1RYD < rfield.emm )
	{
		fprintf( ioQQQ, " This temperature is very low - the blackbody will have significant flux low the low energy limit of the code, presently %10.2e Ryd.\n", 
		  rfield.emm );
		fprintf( ioQQQ, " Was this intended?\n" );
	}

	/* now check that temp not too high - would extend beyond high
	 * energy limit of the code
	 * factor is temperature of 1 Ryd, egamry is high energy limit of code */
	if( rfield.slope[rfield.nspec]/TE1RYD*2. > rfield.egamry )
	{
		fprintf( ioQQQ, " This temperature is very high - the blackbody will have significant flux above the high energy limit of the code,%10.2e Ryd.\n", 
		  rfield.egamry );
		fprintf( ioQQQ, " Was this intended?\n" );
	}

	/* increment continuum indices */
	++rfield.nspec;
	if( rfield.nspec >= LIMSPC )
	{
		fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
		puts( "[Stop in ParseBackgrd]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* also possible to input log(total luminosity)=real log(l) */
	a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	if( lgEOL )
	{
		/* there was not a second number on the line; check if LTE or STE */
		if( lgMatch(" LTE",chCard) || lgMatch("LTE ",chCard) ||
		    lgMatch(" STE",chCard) || lgMatch("STE ",chCard) )
		{
			/* set energy density to the STE - strick thermodynamic equilibrium - value */
			strcpy( rfield.chSpNorm[*nqh], "LUMI" );

			/* need nspec-1 since was incremented above */
			a = log10(rfield.slope[rfield.nspec-1]);
			rlogl = log10(2.99792e10*7.56464e-15) + 4.*a;

			/* set radius to very large value if not already set */
			/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
			if( !radius.lgRadiusKnown )
			{
				*ar1 = (float)radius.rdfalt;
				radius.Radius = pow(10.,radius.rdfalt);
			}
			strcpy( rfield.chRSpec[*nqh], "SQCM" );
			rfield.range[*nqh][0] = rfield.emm;
			rfield.range[*nqh][1] = rfield.egamry;
			rfield.totpow[*nqh] = rlogl;
			*nqh = MIN2(*nqh+1,LIMSPC);

			/* check that stack of shape and luminosity specifications
			 * is parallel, stop if not - this happens is background comes
			 * BETWEEN another set of shape and luminosity commands */
			if( rfield.nspec != *nqh )
			{
				fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete eachcontinuum specification before starting another.\n" );
				fprintf( ioQQQ, " Sorry.\n" );
				puts( "[Stop in ParseBlackbody]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->ParseBlackbody()\n", debug_fp );
#		endif
		return;
	}

	strcpy( rfield.chSpNorm[*nqh], "LUMI" );

	/* a second number was entered, what was it? */
	if( lgMatch("LUMI",chCard) )
	{
		rlogl = a;
		strcpy( rfield.chRSpec[*nqh], "4 PI" );
	}

	else if( lgMatch("RADI",chCard) )
	{
		/* radius was entered, convert to total lumin */
		rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nspec-1]);
		strcpy( rfield.chRSpec[*nqh], "4 PI" );
	}

	else if( lgMatch("DENS",chCard) )
	{
		/* number was temperature to deduce energy density
		 * number is linear if greater than 10, or if LINEAR appears on line
		 * want number to be log of temperature at end of this */
		if( lgMatch("LINE",chCard) || a > 10. )
		{
			a = log10(a);
		}
		rlogl = log10(2.99792e10*7.56464e-15) + 4.*a;
		/* set radius to very large value if not already set */
		/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
		if( !radius.lgRadiusKnown )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.,radius.rdfalt);
		}
		strcpy( rfield.chRSpec[*nqh], "SQCM" );
	}

	else if( lgMatch("DILU",chCard) )
	{
		/* number is dilution factor, if negative then its log */
		if( a < 0. )
		{
			dil = a;
		}
		else
		{
			dil = log10(a);
		}
		if( dil > 0. )
		{
			fprintf( ioQQQ, " Is a dilution factor greater than one physical?\n" );
		}

		/* this is LTE energy density */
		a = log10(rfield.slope[rfield.nspec-1]);
		rlogl = log10(2.99792e10*7.56464e-15) + 4.*a;

		/* add on dilution factor */
		rlogl += dil;

		/* set radius to very large value if not already set */
		/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
		if( !radius.lgRadiusKnown )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.,radius.rdfalt);
		}
		strcpy( rfield.chRSpec[*nqh], "SQCM" );

	}
	else
	{
		fprintf( ioQQQ, " I dont understand what the second number was- keywords are LUMINOSITY, RADIUS, DILUTION, and DENSITY.\n" );
		puts( "[Stop in ParseBlackbody]" );
		cdEXIT(EXIT_FAILURE);
	}

	rfield.range[*nqh][0] = rfield.emm;
	rfield.range[*nqh][1] = rfield.egamry;
	rfield.totpow[*nqh] = rlogl;
	/* the second number will be some sort of luminosity */
	*nqh = MIN2(*nqh+1,LIMSPC);

	/* check that stack of shape and luminosity specifications
	 * is parallel, stop if not - this happens is background comes
	 * BETWEEN another set of shape and luminosity commands */
	if( rfield.nspec != *nqh )
	{
		fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete eachcontinuum specification before starting another.\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in ParseBlackbody]" );
		cdEXIT(EXIT_FAILURE);
	}

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

