/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*ffun evaluate total flux for sum of all continuum sources */
/*ffun1 derive flux at a specific energy, for one continuum */
/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
/*ReadSB99 called to work the table starburst command  */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "called.h"
#include "stars.h"
#include "ipoint.h"
#include "bit32.h"
#include "opacity.h"
#include "continuum.h"
/*lint -e668 passing null pointer */

/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
static void ReadTable(FILE * io);
/*ReadTlusty called by TABLE TLUSTY to read in continuum from a TLUSTY file */
static void ReadTlusty(FILE * io);

/*ReadSB99 called to work the table starburst command  */
static void ReadSB99(FILE * io , double age);
/* these will contain vectors with sb continuum */
static double *sb_energy , *sb_flux;
long int nsp_points;

double ffun(double anu)
{
	double ffun_v;
	static int lgWarn = FALSE;

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

	/* FFUN returns the sum of photons per unit time, area, energy,
	 * for all continuua in the calculation
	 * we loop over all nspec continuum sources - ipspec points to each
	 * */
	ffun_v = 0.;
	for( rfield.ipspec=0; rfield.ipspec < rfield.nspec; rfield.ipspec++ )
	{
		ffun_v += ffun1(anu)*rfield.spfac[rfield.ipspec];
	}

	if( (ffun_v > 1e35 && (!lgWarn)) && bit32.lgBit32 )
	{
		lgWarn = TRUE;
		fprintf( ioQQQ, " FFUN:  The net continuum is very intense for a 32 bit cpu.\n" );
		fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
	}


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

/*ffun1 derive flux at a specific energy, for one continuum */
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( rfield.ipspec >= 0);
	ASSERT( rfield.ipspec < rfield.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 */
	ASSERT( xnu >= rfield.emm*0.99 );
	ASSERT( xnu <= rfield.egamry*1.01 );
	/*end sanity check */

	strcpy( chKey, rfield.chSpType[rfield.ipspec] );

	if( strcmp(chKey,"AGN  ") == 0 )
	{
		/* power law with cutoff at both ends
		 * nomenclature all screwed up - slope is cutoff energy in Ryd,
		 * cutoff[1][i] is ratio of two continua from alpha ox
		 * cutoff[2][i] is slope of bb temp */
		ffun1_v = pow(xnu,-1. + rfield.cutoff[1][rfield.ipspec])*
		  sexp(xnu/rfield.slope[rfield.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,rfield.cutoff[2][rfield.ipspec] - 
				  1.)*rfield.cutoff[0][rfield.ipspec]*sexp(1./
				  xnu);
			}
			else
			{
				ffun1_v += pow(7350.,rfield.cutoff[2][rfield.ipspec] - 
				  1.)*rfield.cutoff[0][rfield.ipspec]/
				  POW3(xnu/7350.);
			}
		}

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

	}
	else if( strcmp(chKey,"BLACK") == 0 )
	{
		const double db_log = log(DBL_MAX);

		/* black body */
		fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
		/* >>>chng 00 apr 13 from 80 to log(dbl_max) */
		if( fac > db_log )
		{
			ffun1_v = 0.;
		}
		else if( fac > 1.e-5 )
		{
			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 >= rfield.tNuRyd[0][rfield.ipspec] )
		{
			/* loop starts at second array element, [1], since want to 
			 * find next continuum energy greater than desired point */
			for( i=1; i < rfield.nupper; i++ )
			{
				if( xnu < rfield.tNuRyd[i][rfield.ipspec] )
				{
					/* points to interpolate upon are in f_nu */
					y = rfield.tFluxLog[i-1][rfield.ipspec] + 
					  rfield.tslop[i-1][rfield.ipspec]*
					  log10(xnu/rfield.tNuRyd[i-1][rfield.ipspec]);
					ASSERT( y >= MIN2( rfield.tFluxLog[i-1][rfield.ipspec] , rfield.tFluxLog[i][rfield.ipspec] )) ;
					ASSERT( y <= MAX2( rfield.tFluxLog[i-1][rfield.ipspec] , rfield.tFluxLog[i][rfield.ipspec] )) ;

					/* return value is photon density, div by energy */
					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/rfield.slope[rfield.ipspec];
		ffun1_v = sexp(fac)/pow(xnu,1.2);

	}
	else if( strcmp(chKey,"LASER") == 0 )
	{
#		ifdef BIG
#			undef BIG
#		endif
#		define	BIG	1e10
#		ifdef SMALL
#			undef SMALL
#		endif
#		define	SMALL	1e-25
		/* a laser, mostly one frequency */
		/* >>chng 01 jul 01, was hard-wired 0.05 rel frac, change to optional
		 * second parameter, with default of 0.05 */
		/*if( xnu > 0.95*rfield.slope[rfield.ipspec] && xnu < 
		  1.05*rfield.slope[rfield.ipspec] )*/
		if( xnu > (1.-rfield.cutoff[0][rfield.ipspec])*rfield.slope[rfield.ipspec] && 
			xnu < (1.+rfield.cutoff[0][rfield.ipspec])*rfield.slope[rfield.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( rfield.tslop[0][rfield.ipspec] == 0. )
		{
			/* >>chng 01 nov 01, array index for below was bogus, only worked for a single continuum */
			ReadTable(rfield.ioTableRead[rfield.ipspec]);
			rfield.tslop[0][rfield.ipspec] = 1.;
		}
		/* use array of values read in on TABLE READ command */
		if( xnu >= rfield.egamry )
		{
			ffun1_v = 0.;
		}
		else
		{
			i = ipoint(xnu);
			if( i > rfield.nupper || i < 1 )
			{
				ffun1_v = 0.;
			}
			else
			{
				ffun1_v = rfield.ConTabRead[i-1]/rfield.anu[i-1]/rfield.anu[i-1];
			}
		}
	}

	else if( strcmp(chKey,"SB99 ") == 0 )
	{
		static long int ipsb=0;
		/* starburst 99 table read
		 * if this is first time called then we need to read in file */
		if( rfield.tslop[1][rfield.ipspec] == 0. )
		{
			/*fprintf(ioQQQ," bug bug age is %f\n", 
				rfield.tslop[0][rfield.ipspec]);*/
			/* >>chng 01 nov 01, array index for below was bogus, only worked for a single continuum */
			/* the file was opened but unread when the table starburst command was parsed,
			 * this call will read the file then close it */
			ReadSB99(
				/* this is the file handle for the already opened starburst file */
				rfield.ioTableRead[rfield.ipspec],
				/* this is the age of the starburst */
				rfield.tslop[0][rfield.ipspec]);
			rfield.tslop[1][rfield.ipspec] = 1.;
		}

		/* use array of values read in on TABLE READ command */
		if( xnu < sb_energy[0] || xnu> sb_energy[nsp_points-1] )
		{
			ffun1_v = 0.;
		}
		else
		{
			/* we want sb_energy[ipsb] <= xnu < sb_energy[ipsb+1] */
			if( sb_energy[ipsb]>xnu )
			{
				/* need to back down */
				while( sb_energy[ipsb]>xnu && ipsb>0 )
					--ipsb;
			}
			if( xnu>=sb_energy[ipsb+1] )
			{
				/* most often xnu will be bigger than it was before,
				 * want to make sure sb_energy[ipsb]<xnu */
				while( xnu>=sb_energy[ipsb+1] && ipsb<nsp_points-1 )
					++ipsb;
			}
			ASSERT( sb_energy[ipsb] <= xnu && xnu < sb_energy[ipsb+1] );

			/* linear interpolation between cells */
			ffun1_v = sb_flux[ipsb] + (sb_flux[ipsb+1]-sb_flux[ipsb])*
				(xnu-sb_energy[ipsb])/(sb_energy[ipsb+1]-sb_energy[ipsb]);
		}

	}
	else if( strcmp(chKey,"TLUST") == 0 )
	{
		/* if this is first time called then we need to read in file */
		if( rfield.tslop[0][rfield.ipspec] == 0. )
		{
			/* >>chng 01 nov 01, array index for below was bogus, only worked for a single continuum */
			ReadTlusty(rfield.ioTableRead[rfield.ipspec]);
			rfield.tslop[0][rfield.ipspec] = 1.;
		}
		/* use array of values read in on TABLE TLUSTY command */
		if( xnu >= rfield.egamry )
		{
			ffun1_v = 0.;
		}
		else
		{
			i = ipoint(xnu);
			if( i > rfield.nupper || i < 1 )
			{
			  ffun1_v = 0.;
			}
			else
			{
			  ffun1_v = rfield.ConTabRead[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 >= rfield.egamry )
		{
			ffun1_v = 0.;
		}
		else
		{
			i = ipoint(xnu);
			if( i > rfield.nupper )
			{
				fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" );
				fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n", 
				  i, rfield.nupper );
				puts( "[Stop in ffun1]" );
				cdEXIT(EXIT_FAILURE);
			}
			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 = rfield.tslop[i-1][rfield.ipspec]/ rfield.anu[i-1];
			}
		}
	}
	else
	{
		fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n", 
		  chKey , rfield.ipspec);
		puts( "[Stop in ffun1]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( (ffun1_v > 1e35 && (!lgWarn)) && bit32.lgBit32 )
	{
		lgWarn = TRUE;
		fprintf( ioQQQ, " FFUN1:  Continuum %ld is very intense for a 32 bit cpu.\n", 
		  rfield.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 );
}

/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
static void ReadTable(FILE * io)
{
	char chLine[INPUT_LINE_LENGTH];
	long int i, 
	  nPoints;
	double Differ, 
	  EnerLast;

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

	/* make sure the file handle points somewhere */
	ASSERT( io != NULL );

	/* read in first line of header */
	if( NULL==fgets( chLine , (int)sizeof(chLine) , io ) )
	{
		fprintf( ioQQQ, " error 1 reading input continuum file.\n" );
		puts( "[Stop in ReadTable]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read in the file of numbers */
	i = 0;
	/* keep reading until we hit eol or run out of room in the array */
	while( (NULL!=fgets(chLine, (int)sizeof(chLine),io)) && (i<rfield.nupper) )
	{
		sscanf( chLine, "%f%f ", &opac.tmn[i], &rfield.ConTabRead[i] );
		++i;
	}
	/* put pointer at last good value */
	nPoints = i - 1;

	/* check that sane number of values entered */
	if( nPoints < 1 )
	{
		fprintf( ioQQQ, " ReadTable, file for TABLE READ has too few points, =%5ld\n", 
		  nPoints );
		puts( "[Stop in ReadTable]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check on units of energy scale, convert to Rydbergs if necessary
	 * tmn is scale read in from punch file, anu is correct scale
	 * EnerLast is energy of last point in rydbergs */
	EnerLast = opac.tmn[nPoints];
	if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 )
	{
		/* first energy does not appear to have been in Rydbergs */
		if( opac.tmn[0] < 0. )
		{
			fprintf( ioQQQ, " ReadTable:  the first energy in table read file is not positive.  Something is wrong.\n" );
			puts( "[Stop in ReadTable]" );
			cdEXIT(EXIT_FAILURE);
		}
		else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0]*1e-4)/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in microns */
			EnerLast = RYDLAM/opac.tmn[nPoints]/1e4;
		}
		else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0])/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in Angstroms */
			EnerLast = RYDLAM/opac.tmn[nPoints];
		}
		else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD*1e-3)/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in keV */
			EnerLast = opac.tmn[nPoints]/EVRYD/1e-3;
		}
		else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD)/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in eV */
			EnerLast = opac.tmn[nPoints]/EVRYD;
		}
		else
		{
			fprintf( ioQQQ, " ReadTable:  the energy scale in the table read file makes no sense to me.\n" );
			puts( "[Stop in ReadTable]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* now check now the energies of the highest points agree */
	Differ = fabs(EnerLast-rfield.anu[nPoints])/rfield.anu[nPoints];
	if( Differ > 0.001 )
	{
		fprintf( ioQQQ, " ReadTable: The energy mesh of the file read in by the TABLE READ command does not agree with this version of Cloudy.\n" );
		fprintf( ioQQQ, " ReadTable: Was the file generated by an older version of the code?\n" );
		fprintf( ioQQQ, " ReadTable: Number of points read in=%5ld\n", 
		  nPoints );
		fprintf( ioQQQ, " ReadTable: input, internal energies=%12.4e%12.4e\n", 
		  opac.tmn[nPoints], rfield.anu[nPoints] );
		puts( "[Stop in ReadTable]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* make sure rest of array has valid zeros */
	for( i=nPoints + 1; i < rfield.nupper; i++ )
	{
		rfield.ConTabRead[i] = 0.;
	}

	fclose(io);

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

/*ReadSB99 called to work the table starburst command  */
static void ReadSB99(FILE * io , double age)
{
	char chLine[INPUT_LINE_LENGTH];
	long int i, 
	  nCont;
	double fmax ,
	  t1, 
	  t2 , 
	  agehi, 
	  agelo , 
	  ageold , 
	  fac;
	double *wl , *flo , *fhi , *finterp ;
	int lgReadHiDone;

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

	/* make sure the file handle points somewhere */
	ASSERT( io != NULL );

	/* line just before first numbers starts with TIME *
	 * keep reading until we hit eol find TIME */
	while( (NULL!=fgets(chLine, (int)sizeof(chLine),io))  )
	{
		if( strncmp( chLine , " TIME" , 5 )==0 )
			break;
	}

	/* now count number of points in stellar continuum */
	if( (NULL==fgets(chLine, (int)sizeof(chLine),io))  )
	{
		fprintf( ioQQQ, " ReadSB99, file fail 2 \n" );
		puts( "[Stop in ReadSB99]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* get first age */
	sscanf( chLine, "%le ", &t1 );
	nsp_points = 1;
	/* now count how many wavelength cells have this same age */
	while( (NULL!=fgets(chLine, (int)sizeof(chLine),io))  )
	{
		sscanf( chLine, "%le ", &t2 );
		if( fabs(t2 -t1)/t1 > 1e-4 )
			break;
		++nsp_points;
	}
	/* now have the number of continuum points */
	/* create vectors to store wl and flux */
	if( (wl = (double*)MALLOC(sizeof(double)*(unsigned)nsp_points ))==NULL )
		BadMalloc();
	if( (flo = (double*)MALLOC(sizeof(double)*(unsigned)nsp_points ))==NULL )
		BadMalloc();
	if( (fhi = (double*)MALLOC(sizeof(double)*(unsigned)nsp_points ))==NULL )
		BadMalloc();
	if( (finterp = (double*)MALLOC(sizeof(double)*(unsigned)nsp_points ))==NULL )
		BadMalloc();

	/* these two are file static */
	if( (sb_flux = (double*)MALLOC(sizeof(double)*(unsigned)nsp_points ))==NULL )
		BadMalloc();
	if( (sb_energy = (double*)MALLOC(sizeof(double)*(unsigned)nsp_points ))==NULL )
		BadMalloc();

	/* now rewind the file so we can read it a second time*/
	if( fseek( io , 0 , SEEK_SET ) != 0 )
	{
		fprintf( ioQQQ, " ReadSB99 could not rewind.\n");
		puts( "[Stop in ReadSB99]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read down to TIME again */
	while( (NULL!=fgets(chLine, (int)sizeof(chLine),io))  )
	{
		if( strncmp( chLine , " TIME" , 5 )==0 )
			break;
	}

	ageold = 0;
	lgReadHiDone = FALSE;
	nCont = 0;
	agelo = 0.;
	agehi = 0.;

	while( !lgReadHiDone )
	{
		/* read in nsp_points at a time - this should be one age block */
		for( i=0; i<nsp_points; ++i )
		{
			if( (NULL==fgets(chLine, (int)sizeof(chLine),io))  )
			{
				if( i==0 )
				{
					fprintf( ioQQQ, " ReadSB99 could not read first continuum point in next continuum while searching for correct age.\n");
					fprintf( ioQQQ, " ReadSB99 requested age beyond range of ages in data file.\n");
				}
				else
				{
					fprintf( ioQQQ, " ReadSB99 could not read continuum point - could data file be corrupted?\n");
				}
				puts( "[Stop in ReadSB99]" );
				cdEXIT(EXIT_FAILURE);
			}
			sscanf( chLine, "%le %lf %lf ", &agehi , &wl[i] , &fhi[i] );
			/* if first point set baseline time to this value */
			if( i==0 )
			{
				if( ageold > agehi )
				{
					fprintf( ioQQQ, " ReadSB99 the burst ages were not in increasing order.\n");
					puts( "[Stop in ReadSB99]" );
					cdEXIT(EXIT_FAILURE);
				}
				t1 = agehi;
			}
			if( fabs(agehi -t1)/t1 > 1e-4 )
			{
				fprintf( ioQQQ, " ReadSB99 error reading spectrum - time changed within spectrum.\n");
				puts( "[Stop in ReadSB99]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* now have the block, see if time of continuum is less than requested age */
		if( agehi < age   )
		{
			/* block we just read in is still too young - copy to agelo */
			agelo = agehi;
			/* copy fhi continuum in flo */
			for( i=0; i<nsp_points; ++i )
			{
				flo[i] = fhi[i];
			}
		}
		else if( agehi>= age )
		{
			/* have greater age, but first make sure lesser age already read in */
			if( agelo==0. )
			{
				fprintf( ioQQQ, " ReadSB99 error interpolating on spectrum - specified age less than first age.\n");
				puts( "[Stop in ReadSB99]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* we have bracketed the target - now have greater age and are done */
			lgReadHiDone = TRUE;
		}
		++nCont;
	}

	if( called.lgTalk )
		fprintf( ioQQQ, 
		" ReadSB99 interpolating betweeen continuum %li, age=%.3e yr, and continuum age %.3e yr for requested age %.3e yr.\n\n" ,
		nCont-1, agelo , agehi , age );

	/* now interpolate, linear in log log space */
	fac = (log10(age) - log10(agelo)) / (log10(agehi)-log10(agelo));
	if( fac < 0. || fac > 1. )
		TotalInsanity();

	fmax = 0.;
	for( i=0; i<nsp_points; ++i )
	{
		finterp[i] = flo[i] + fac*(fhi[i]-flo[i]);
		fmax = MAX2( fmax , finterp[i] );
	}

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			for( i=0; i<nsp_points; ++i )
			{
				fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\t%.3e\n",
					wl[i] , flo[i] , fhi[i] , finterp[i] );
			}
		}
	}

	/* now convert wl and flux in cloudy units */
	for( i=0; i<nsp_points; ++i )
	{
		/* energy in Ryd in increasing order */
		sb_energy[i] = RYDLAM / wl[nsp_points - 1 - i]; 
		/* original units were erg /s /A , convert to nu F nu */
		sb_flux[i] = wl[nsp_points - 1 - i] * pow(10. , (finterp[nsp_points - 1 - i]-fmax) );
		/* now convert to photons per ryd, div by two energies */
		sb_flux[i] /= POW2( sb_energy[i] );
		/*fprintf(ioQQQ,"%.3e\t%.3e\n", sb_energy[i] , sb_flux[i] );*/
	}

	/* no longer need flo or fhi */
	free( flo );
	free( fhi );
	free( wl );
	free( finterp );

	fclose(io);

#	ifdef DEBUG_FUN
	fputs( " <->ReadSB99()\n", debug_fp );
#	endif
	return;
}
/*lint +e668 passing null pointer */

#define NTLUSTY 200000

/*ReadTlusty called by TABLE TLUSTY to read in continuum from a TLUSTY spectrum file (Added 5 May 2003 by K. Volk) */
static void ReadTlusty(FILE * ipLUSTY)
{

  double *freqtlusty,*fnutlusty;
  double tempvalue,offset;
  long i,j,k,noffset;
  float *starener, *starflux, *cloudyflux;
  float edges[3] = {0.9994678f,1.8071406f,3.9996377f};
  long nedges = 3;
  char tlline[80];


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

	/* make sure the file handle points somewhere */
	ASSERT( ipLUSTY != NULL );

	/* Read the TLUSTY file.  First, allocate the arrays needed. */

	freqtlusty=(double *) malloc(NTLUSTY*sizeof(double));
	if (freqtlusty == NULL) 
	{
		printf( " not enough memory to allocate freqtlusty in ReadTlusty\n" );
		puts( "[Stop in ReadTlusty]" );
		cdEXIT(EXIT_FAILURE);
	}
	fnutlusty=(double *) malloc(NTLUSTY*sizeof(double));
	if (freqtlusty == NULL) 
	{
		printf( " not enough memory to allocate fnutlusty in ReadTlusty\n" );
		puts( "[Stop in ReadTlusty]" );
		cdEXIT(EXIT_FAILURE);
	}

	j = -1000;
	for (i=0; i < NTLUSTY; i++)
	{
		if (fgets(tlline,80,ipLUSTY) != NULL)
		{
			for (k=0; k < 80; k++)
			{
				if (tlline[k]=='D')
					tlline[k]='e';
			}
			k = sscanf(tlline,"%le %le\n",&freqtlusty[i],&fnutlusty[i]);
			if(k != 2 || fnutlusty[i] <= 0. || freqtlusty[i] <= 0.)
			{
				fprintf( ioQQQ, " I could not read TLUSTY file properly.  Bad values on line %li\n", i+1 );
				puts( "[Stop in ReadTlusty]" );
				fclose(ipLUSTY);
				cdEXIT(EXIT_FAILURE);
			}
		}
		else
		{
			j = i-1;
			if (j < 100000)
			{
				fprintf( ioQQQ, " I could not read TLUSTY file properly.  Too few values were read in.\n");
				puts( "[Stop in ReadTlusty]" );
				fclose(ipLUSTY);
				cdEXIT(EXIT_FAILURE);
			}
			else
			{
				for (i=0; i <= j; i++)
				{
					freqtlusty[i]=freqtlusty[i]/3.289841960367e15;
				}
				/* set the values in ascending energy order */
				for (i=0; i < j/2; i++)
				{
					tempvalue=freqtlusty[i];
					freqtlusty[i] = freqtlusty[j-i];
					freqtlusty[j-i] = tempvalue;
					tempvalue = fnutlusty[i];
					fnutlusty[i] = fnutlusty[j-i];
					fnutlusty[j-i] = tempvalue;
				}
				/* The following is to check for duplicate energy values 
				   in the TLUSTY spectrum: will crash RebinAtmosphere. When
				   found, they are changed to slightly separated values. */
				for (i=0; i < j; i++)
				{
					if (freqtlusty[i] == freqtlusty[i+1])
					{
						if (i == 0)
						{
							tempvalue = freqtlusty[3]-freqtlusty[2];
						}
						else
						{
							if (i == j-1)
							{
								tempvalue = freqtlusty[i]-freqtlusty[i-1];
							}
							else
							{
								tempvalue = (freqtlusty[i+1]-freqtlusty[i-1])/2.;
							}
						}
						offset = tempvalue/200.;
						for (noffset=1;  noffset > 0; noffset++)
						{
							tempvalue = freqtlusty[i]+offset;
							if ((float)tempvalue == (float)freqtlusty[i])
							{
								offset = offset*1.5;
							}
							else
							{
								noffset =  -10;
							}
						}
						freqtlusty[i] = freqtlusty[i]-offset;
						freqtlusty[i+1] = freqtlusty[i+1]+offset;
					}
				}
				i = NTLUSTY+1;
			}
		}
	}

	/* Now, change values to FLOAT and call RebinAtmosphere */

	starener = (float *) malloc(NTLUSTY*sizeof(float));
	if (starener == NULL)
	{
		printf( " not enough memory to allocate starener in ReadTlusty\n" );
		puts( "[Stop in ReadTlusty]" );
		cdEXIT(EXIT_FAILURE);
	}
	starflux = (float *) malloc(NTLUSTY*sizeof(float));
	if (starflux == NULL)
	{
		printf( " not enough memory to allocate starflux in ReadTlusty\n" );
		puts( "[Stop in ReadTlusty]" );
		cdEXIT(EXIT_FAILURE);
	}
	cloudyflux = (float *) malloc(NTLUSTY*sizeof(float));
	if (cloudyflux == NULL)
	{
		printf( " not enough memory to allocate cloudyflux in ReadTlusty\n" );
		puts( "[Stop in ReadTlusty]" );
		cdEXIT(EXIT_FAILURE);
	}

	ASSERT( j > -10 );
	tempvalue = -1000;
	for (i=0; i < j; i++)
	{
		if (freqtlusty[i] < 1. && freqtlusty[i+1] >= 1.)
		{
			tempvalue = fnutlusty[i];
			i = j+1;
		}
	}
	/* This should not happen...but check for a zero value */
	ASSERT( tempvalue > 0.);

	for (i=0; i <= j; i++)
	{
		fnutlusty[i] = fnutlusty[i]/tempvalue;
	}

	for (i=0; i <= j; i++)
	{
		starener[i] = (float)freqtlusty[i];
		starflux[i] = (float)fnutlusty[i];
	}
	RebinAtmosphere(j,starener,starflux,cloudyflux,nedges,edges);

	for (i = 0; i < rfield.nupper; i++)
	{
		rfield.ConTabRead[i] = cloudyflux[i];
	}


	fclose(ipLUSTY);

	/* free the malloced values */
	free(freqtlusty);
	free(fnutlusty);
	free(starener);
	free(starflux);
	free(cloudyflux);

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