/* 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 */
/*ParseAgn parse parameters for the AGN continuum shape command */
#include "cddefines.h"
#include "rfield.h"
#include "parse.h"
#include "physconst.h"

void ParseAgn(char *chCard)
{
	int lgEOL;
	long int i;
	double BigBump, 
	  Ratio, 
	  XRays, 
	  xnu;

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

	/* this radiation field will be something like an AGN */
	strcpy( rfield.chSpType[rfield.nspec], "AGN  " );

	/* there were no numbers on the line - this could be the Kirk option,
		* to use the numbers for the continuum in the atlas paper */
	if( lgMatch("KIRK",chCard) )
	{
		/* million degree cutoff, but in rydbergs */
		rfield.slope[rfield.nspec] = 1e6 /  TE1RYD;

		/* cutoff is second parameter is really alpha ox */
		rfield.cutoff[rfield.nspec][0] = -1.40;

		/* bb slope is third parameter */
		rfield.cutoff[rfield.nspec][1] = -0.50;

		/* slope of X-Ray component is last parameter */
		rfield.cutoff[rfield.nspec][2] = -1.0;
	}
	else
	{
		/* first parameter is temperature of big bump
		* second parameter is alpha ox */
		i = 5;
		/* slope is first parameter is really temperature of bump */
		rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{

			fprintf( ioQQQ, " The big bump temperature should have been on this line.   Sorry.\n" );
			puts( "[Stop in ParseAgn]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( rfield.slope[rfield.nspec] <= 0. )
		{
			fprintf( ioQQQ, " Non positive temperature not allowed.   Sorry.\n" );
			puts( "[Stop in ParseAgn]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* temps are log if first le 10 */
		if( rfield.slope[rfield.nspec] <= 10. )
			rfield.slope[rfield.nspec] = 
			pow(10.,rfield.slope[rfield.nspec]);

		/* want cutoff in ryd not kelvin */
		rfield.slope[rfield.nspec] /= TE1RYD;

		/* cutoff is second parameter is really alpha ox */
		rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " alpha ox should have been on this line.   Sorry.\n" );
			puts( "[Stop in ParseAgn]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( rfield.cutoff[rfield.nspec][0] > 3. || 
			rfield.cutoff[rfield.nspec][0] < -3. )
		{
			fprintf( ioQQQ, " An alpha ox of%10.2e looks funny to me.  Check Hazy to make sure its ok.\n", 
			rfield.cutoff[rfield.nspec][0] );
		}

		if( rfield.cutoff[rfield.nspec][0] >= 0. )
		{
			fprintf( ioQQQ, " The sign of alpha ox is almost certainly incorrect.   Check Hazy.\n" );
		}

		/* bb slope is third parameter */
		rfield.cutoff[rfield.nspec][1] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
		if( lgEOL )
			rfield.cutoff[rfield.nspec][1] = -0.5f;

		/* slope of X-Ray component is last parameter */
		rfield.cutoff[rfield.nspec][2] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
		if( lgEOL )
			rfield.cutoff[rfield.nspec][2] = -1.0f;
	}

	/* 403.3 is ratio of energies where alpha ox defined,
	 * assumed to be 2500A and 2keV */
	Ratio = pow(403.3,rfield.cutoff[rfield.nspec][0] - 1.);

	/* following code must be kept parallel with ffun1 */
	xnu = 0.3645;
	BigBump = pow(xnu,-1. + rfield.cutoff[rfield.nspec][1])*
	  sexp(xnu/rfield.slope[rfield.nspec]);
	xnu = 147.;

	/* XRays = xnu**(-2.) */
	XRays = pow(xnu,rfield.cutoff[rfield.nspec][2] - 1.);
	if( BigBump <= 0. )
	{
		fprintf( ioQQQ, " Big Bump had zero flux at .3645 Ryd.\n" );
		puts( "[Stop in ParseAgn]" );
		cdEXIT(EXIT_FAILURE);
	}
	rfield.cutoff[rfield.nspec][0] = (Ratio/(XRays/BigBump));

	/* lastly increment number of spectra and check that we still have room in the vector */
	++rfield.nspec;
	if( rfield.nspec >= LIMSPC )
	{
		fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
		puts( "[Stop in ParseAgn]" );
		cdEXIT(EXIT_FAILURE);
	}

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

