/*ffun1 derive flux at a specific energy, for one continuum */
#include <string.h>
#include "cddefines.h"
#	ifdef BIG
#undef BIG
#	endif
#define	BIG	1e10
#	ifdef SMALL
#undef SMALL
#	endif
#define	SMALL	1e-25
#include "showme.h"
#include "physconst.h"
#include "incidspec.h"
#include "rfield.h"
#include "cnread.h"
#include "bounds.h"
#include "bit32.h"
#include "sexp.h"
#include "ipoint.h"
#include "readtable.h"
#include "ffun1.h"

double ffun1(double xnu)
{
	char chKey[6];
	long int i;
	double fac, 
	  ffun1_v, 
	  y;
	static int lgWarn = FALSE;

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


	/* confirm that pointer is within range */
	assert( IncidSpec.ipspec >= 0);
	assert( IncidSpec.ipspec < IncidSpec.nspec );

	/* FFUN1 returns photons per unit time, area, energy, for one continuum
	 * ipspec is the pointer to the continuum source, in the order
	 * entered on the command lines
	 *
	 *begin sanity check */
	if( xnu < bounds.emm*0.99 )
	{
		fprintf( ioQQQ, " fun1 called with impossible energy%11.3e\n", 
		  xnu );
		fprintf( ioQQQ, " the low energy limit was %11.3e\n", 
		  bounds.emm );
		ShowMe();
		puts( "[Stop in ffun1]" );
		exit(1);
	}
	else if( xnu > bounds.egamry*1.01 )
	{
		fprintf( ioQQQ, " fun1 called with impossible energy%11.3e\n", 
		  xnu );
		fprintf( ioQQQ, " the high energy limit was %11.3e\n", 
		  bounds.egamry );
		ShowMe();
		puts( "[Stop in ffun1]" );
		exit(1);
	}
	/*end sanity check */
	strcpy( chKey, spnorm.chSpType[IncidSpec.ipspec] );

	if( strcmp(chKey,"AGN  ") == 0 )
	{
		/* power law with cutoff at both ends
		 * nomenclature all screwed up - slope is cutoff energy in Ryd,
		 * cutoff(i,1) is ratio of two continua from alpha ox
		 * cutoff(i,2) is slope of bb temp */
		ffun1_v = pow(xnu,-1. + IncidSpec.cutoff[1][IncidSpec.ipspec])*
		  sexp(xnu/IncidSpec.slope[IncidSpec.ipspec])*sexp(0.01/
		  xnu);
		/* only add on x-ray for energies above 0.1 Ryd */
		if( xnu > 0.1 )
		{
			if( xnu < 7350. )
			{
				/* cutoff is actually normalization constant
				 * below 100keV continuum is nu-1 */
				ffun1_v += pow(xnu,IncidSpec.cutoff[2][IncidSpec.ipspec] - 
				  1.)*IncidSpec.cutoff[0][IncidSpec.ipspec]*sexp(1./
				  xnu);
			}
			else
			{
				ffun1_v += pow(7350.,IncidSpec.cutoff[2][IncidSpec.ipspec] - 
				  1.)*IncidSpec.cutoff[0][IncidSpec.ipspec]/
				  POW3(xnu/7350.);
			}
		}

	}
	else if( strcmp(chKey,"POWER") == 0 )
	{
		/* power law with cutoff at both ends */
		ffun1_v = pow(xnu,-1. + IncidSpec.slope[IncidSpec.ipspec])*
		  sexp(xnu/IncidSpec.cutoff[0][IncidSpec.ipspec]+IncidSpec.cutoff[1][IncidSpec.ipspec]/
		  xnu);

	}
	else if( strcmp(chKey,"BLACK") == 0 )
	{
		/* black body */
		fac = TE1RYD*xnu/IncidSpec.slope[IncidSpec.ipspec];
		if( fac >= 80. )
		{
			ffun1_v = 0.;
		}
		else if( fac < 80. && fac > 0.01 )
		{
			ffun1_v = xnu*xnu/(exp(fac) - 1.);
		}
		else
		{
			ffun1_v = xnu*xnu/(fac*(1. + fac/2.));
		}

	}
	else if( strcmp(chKey,"INTER") == 0 )
	{
		/* interpolate on input spectrum */
		if( xnu >= IncidSpec.tnu[0][IncidSpec.ipspec] )
		{
			for( i=1; i < NCELL; i++ )
			{
				if( xnu < IncidSpec.tnu[i][IncidSpec.ipspec] )
				{
					y = IncidSpec.tfac[i-1][IncidSpec.ipspec] + 
					  IncidSpec.tslop[i-1][IncidSpec.ipspec]*
					  log10(xnu/IncidSpec.tnu[i-1][IncidSpec.ipspec]);
					ffun1_v = (pow(10.,y))/xnu;
					
#	ifdef DEBUG_FUN
					fputs( " <->ffun1()\n", debug_fp );
#	endif
					return( ffun1_v );
				}
			}
			/* energy above highest in table */
			ffun1_v = 0.;
		}
		else
		{
			/* energy below lowest on table */
			ffun1_v = 0.;
		}

	}
	else if( strcmp(chKey,"BREMS") == 0 )
	{
		/* brems continuum, rough gaunt factor */
		fac = TE1RYD*xnu/IncidSpec.slope[IncidSpec.ipspec];
		ffun1_v = sexp(fac)/pow(xnu,1.2);

	}
	else if( strcmp(chKey,"LASER") == 0 )
	{
		/* a laser, mostly one frequency */
		if( xnu > 0.95*IncidSpec.slope[IncidSpec.ipspec] && xnu < 
		  1.05*IncidSpec.slope[IncidSpec.ipspec] )
		{
			ffun1_v = BIG;
		}
		else
		{
			ffun1_v = SMALL;
		}

	}
	else if( strcmp(chKey,"READ ") == 0 )
	{
		/* if this is first time called then we need to read in file */
		if( IncidSpec.tslop[0][IncidSpec.ipspec] == 0. )
		{
			ReadTable(spnorm.ioTableRead[IncidSpec.nspec-1]);
			IncidSpec.tslop[0][IncidSpec.ipspec] = 1.;
		}
		/* use array of values read in on TABLE READ command */
		if( xnu >= bounds.egamry )
		{
			ffun1_v = 0.;
		}
		else
		{
			i = ipoint(xnu);
			if( i > rfield.nupper || i < 1 )
			{
				ffun1_v = 0.;
			}
			else
			{
				ffun1_v = cnread.ConTabRead[i-1]/rfield.anu[i-1]/rfield.anu[i-1];
			}
		}

	}
	else if( strcmp(chKey,"VOLK ") == 0 )
	{
		/* use array of values read in from Kevin Volk's rebinning of
		 * large atlas grids */
		if( xnu >= bounds.egamry )
		{
			ffun1_v = 0.;
		}
		else
		{
			i = ipoint(xnu);
			if( i > NCELL )
			{
				fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" );
				fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n", 
				  i, NCELL );
				puts( "[Stop in ffun1]" );
				exit(1);
			}
			if( i > rfield.nupper || i < 1 )
			{
				ffun1_v = 0.;
			}
			else
			{
				/* bug fixed Jul 9 93: FFUN1 = TSLOP(IPSPEC,I) / ANU(I) / ANU(I)
				 *   i has value 939 */
				ffun1_v = IncidSpec.tslop[i-1][IncidSpec.ipspec]/
				  rfield.anu[i-1];
			}
		}
	}
	else
	{
		fprintf( ioQQQ, " ffun1: I do not understand continuum label %5.5s\n", 
		  chKey );
		puts( "[Stop in ffun1]" );
		exit(1);
	}

	if( (ffun1_v > 1e35 && (!lgWarn)) && bit32.lgBit32 )
	{
		lgWarn = TRUE;
		fprintf( ioQQQ, " FFUN1:  Continuum%4ld is very intense for a 32 bit cpu.\n", 
		  IncidSpec.ipspec );
		fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
	}


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