/*ParsePowerlawContinuum parse the power law continuum command */
#include <string.h>
#include "cddefines.h"
#include "incidspec.h"
#include "varypar.h"
#include "input.h"
#include "ffmtread.h"
#include "lgmatch.h"
#include "parse.h"

void ParsePowerlawContinuum(char *chCard)
{
	int lgEOL;
	long int i;

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

	/* power law with cutoff and X-ray continuum */
	IncidSpec.nspec += 1;
	if( IncidSpec.nspec > LIMSPC )
	{
		fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
		puts( "[Stop in ParsePowerlawContinuum]" );
		exit(1);
	}

	strcpy( spnorm.chSpType[IncidSpec.nspec-1], "POWER" );

	/* first parameter is slope of continuum, probably should be negative */
	i = 5;
	IncidSpec.slope[IncidSpec.nspec-1] = (float)FFmtRead(chCard,&i,76,&lgEOL);
	if( lgEOL )
	{
		fprintf( ioQQQ, " There should have been a number on this line.   Sorry.\n" );
		puts( "[Stop in ParsePowerlawContinuum]" );
		exit(1);
	}

	if( IncidSpec.slope[IncidSpec.nspec-1] >= 0. )
	{
		fprintf( ioQQQ, " Is the slope of this power law correct?\n" );
	}

	/* second optional parameter is high energy cut off */
	IncidSpec.cutoff[0][IncidSpec.nspec-1] = (float)FFmtRead(chCard,&i,76,
	  &lgEOL);

	/* no cutoff if eof hit */
	if( lgEOL )
	{
		/* no extra parameters at all, so put in extreme cutoffs */
		IncidSpec.cutoff[0][IncidSpec.nspec-1] = 1e4;
		IncidSpec.cutoff[1][IncidSpec.nspec-1] = 1e-4f;
	}
	else
	{
		/* first cutoff was present, check for second */
		IncidSpec.cutoff[1][IncidSpec.nspec-1] = (float)FFmtRead(chCard,&i,
		  76,&lgEOL);
		if( lgEOL )
			IncidSpec.cutoff[1][IncidSpec.nspec-1] = 1e-4f;
	}

	/* check that energies were entered in the correct order */
	if( IncidSpec.cutoff[1][IncidSpec.nspec-1] > IncidSpec.cutoff[0][IncidSpec.nspec-1] )
	{
		fprintf( ioQQQ, " The optional cutoff energies do not appear to be in the correct order.  Check Hazy.\n" );
		puts( "[Stop in ParsePowerlawContinuum]" );
		exit(1);
	}

	/* check for keywork KELVIN to interprete cutoff energies as degrees */
	if( lgMatch("KELV",chCard) )
	{
		/* temps are log if first le 10 */
		if( IncidSpec.cutoff[0][IncidSpec.nspec-1] <= 10. )
			IncidSpec.cutoff[0][IncidSpec.nspec-1] = (float)pow(10.,IncidSpec.cutoff[0][IncidSpec.nspec-1]);
		if( IncidSpec.cutoff[1][IncidSpec.nspec-1] <= 10. )
			IncidSpec.cutoff[1][IncidSpec.nspec-1] = (float)pow(10.,IncidSpec.cutoff[1][IncidSpec.nspec-1]);
		IncidSpec.cutoff[0][IncidSpec.nspec-1] /= 1.5789e5;
		IncidSpec.cutoff[1][IncidSpec.nspec-1] /= 1.5789e5;
	}

	if( IncidSpec.cutoff[0][IncidSpec.nspec-1] < 0. || IncidSpec.cutoff[1][IncidSpec.nspec-1] < 
	  0. )
	{
		fprintf( ioQQQ, " A negative cutoff energy is not physical.  Sorry.\n" );
		puts( "[Stop in ParsePowerlawContinuum]" );
		exit(1);
	}

	if( IncidSpec.cutoff[1][IncidSpec.nspec-1] == 0. && 
		IncidSpec.slope[IncidSpec.nspec-1] <= -1. )
	{
		fprintf( ioQQQ, " A power-law with this slope, and no low energy cutoff, may have an unphysically large\n brightness temperature in the radio.\n" );
	}

	/* vary option */
	if( VaryPar.lgVarOn )
	{
		/* pointer to where to write */
		VaryPar.nvfpnt[VaryPar.nparm] = input.nRead;
		if( lgMatch("ARYB",chCard) )
		{
			/* this test is for key "varyb", meaning to vary second parameter
			 * the cutoff temperature
			 * this is the number of parameters to feed onto the input line */
			VaryPar.nvarxt[VaryPar.nparm] = 2;
			/* vary ?? */
			sprintf( chVaryPar.chVarFmt[VaryPar.nparm], 
				"POWER LAW %f KELVIN%%f %%f", 
			  IncidSpec.slope[IncidSpec.nspec-1] );
			/* param is linear scale factor */
			VaryPar.vparm[0][VaryPar.nparm] = (float)log10(IncidSpec.cutoff[0][IncidSpec.nspec-1]*
			  1.5789e5);
			VaryPar.vparm[1][VaryPar.nparm] = (float)log10(IncidSpec.cutoff[1][IncidSpec.nspec-1]*
			  1.5789e5);
			VaryPar.vincr[VaryPar.nparm] = 0.2f;
		}
		else if( lgMatch("ARYC",chCard) )
		{
			/* the keyword wsa "varyc"
			 * this is the number of parameters to feed onto the input line */
			VaryPar.nvarxt[VaryPar.nparm] = 1;

			/* vary ?? */
			sprintf( chVaryPar.chVarFmt[VaryPar.nparm], "POWER LAW%f %f KELVIN %%f )", 
			  IncidSpec.slope[IncidSpec.nspec-1], 
			  log10(IncidSpec.cutoff[0][IncidSpec.nspec-1]* 1.5789e5) );

			VaryPar.vparm[0][VaryPar.nparm] = (float)log10(IncidSpec.cutoff[1][IncidSpec.nspec-1]*
			  1.5789e5);
			VaryPar.vincr[VaryPar.nparm] = 0.2f;
		}
		else
		{
			/* vary the first parameter only, but still are two more
			 * this is the number of parameters to feed onto the input line */
			VaryPar.nvarxt[VaryPar.nparm] = 3;
			strcpy( chVaryPar.chVarFmt[VaryPar.nparm], 
				"POWER LAW KELVIN%f %f %f" );
			/* param is slope of the power law */
			VaryPar.vparm[0][VaryPar.nparm] = IncidSpec.slope[IncidSpec.nspec-1];
			VaryPar.vparm[1][VaryPar.nparm] = (float)log10(IncidSpec.cutoff[0][IncidSpec.nspec-1]*
			  1.5789e5);
			VaryPar.vparm[2][VaryPar.nparm] = (float)log10(IncidSpec.cutoff[1][IncidSpec.nspec-1]*
			  1.5789e5);
			VaryPar.vincr[VaryPar.nparm] = 0.2f;
		}
		++VaryPar.nparm;
	}


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