/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  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"
#include "parser.h"

void ParseAgn(Parser &p)
{
	double BigBump, 
	  Ratio, 
	  XRays, 
	  xnu;

	DEBUG_ENTRY( "ParseAgn()" );

	/* this radiation field will be something like an AGN */
	strcpy( rfield.chSpType[rfield.nShape], "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( p.nMatch("KIRK") )
	{
		/* million degree cutoff, but in rydbergs */
		rfield.slope[rfield.nShape] = 1e6 /  TE1RYD;

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

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

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

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

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

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

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

		/* cutoff is second parameter is really alpha ox */
		rfield.cutoff[rfield.nShape][0] = p.FFmtRead();
		if( p.lgEOL() )
		{
			fprintf( ioQQQ, " alpha ox should have been on this line.   Sorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}

		if( rfield.cutoff[rfield.nShape][0] > 3. || 
			rfield.cutoff[rfield.nShape][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.nShape][0] );
		}

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

		/* bb slope is third parameter */
		rfield.cutoff[rfield.nShape][1] = p.FFmtRead();
		if( p.lgEOL() )
			rfield.cutoff[rfield.nShape][1] = -0.5f;

		/* slope of X-Ray component is last parameter */
		rfield.cutoff[rfield.nShape][2] = p.FFmtRead();
		if( p.lgEOL() )
			rfield.cutoff[rfield.nShape][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.nShape][0] - 1.);

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

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

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