/* 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 */
/*ParseTable parse the table read command */
/*lines_table envoked by table lines command, check if we can find all lines in a given list */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "optimize.h"
#include "rfield.h"
#include "trace.h"
#include "radius.h"
#include "input.h"
#include "stars.h"
#include "lines.h"
#include "parse.h"

/* these will become the label and wl for a possible list of lines,
 * obtained when tables lines used */
static char **chLabel ;
static float *wl;
static long int nLINE_TABLE = 0;
static char chLINE_LIST[FILENAME_PATH_LENGTH];

/* values for these are established in zerocontinuum below -
 * they have file scope and are not used elsewhere */
#define	NCRAB	10
static double tnucrb[NCRAB], 
  fnucrb[NCRAB];

/* Bob Rubin's fixed theta 1 Ori C continuum */
#define	NRUBIN	56
double tnurbn[NRUBIN] = {1.05E-08,1.05E-07,1.05E-06,1.04E-05,1.00E-04,1.00E-03,1.00E-02,3.01E-02,1.00E-01,
	1.50E-01,2.50E-01,4.01E-01,6.01E-01,9.8E-01,9.96E-01,1.00E+00,1.02445,1.07266,1.12563,1.18411,1.23881,
	1.29328,1.35881,1.42463,1.48981,1.55326,1.6166,1.68845,1.76698,1.8019,1.808,1.84567,1.9317,2.04891,2.14533,
	2.19702,2.27941,2.37438,2.43137,2.49798,2.56113,2.59762,2.66597,2.80543,2.95069,3.02911,3.11182,3.22178,
	3.3155,3.42768,3.50678,3.56157,3.61811,3.75211,3.89643,4.}, 
	fnurbn[NRUBIN] = {1.56E-20,1.55E-17,1.54E-14,1.53E-11,1.35E-08,1.34E-05,1.35E-02,3.62E-01,1.27E+01,
	3.90E+01,1.48E+02,4.52E+02,1.02E+03,2.27E+03,8.69E+02,8.04E+02,6.58E+02,6.13E+02,6.49E+02,6.06E+02,
	6.30E+02,5.53E+02,5.55E+02,5.24E+02,4.86E+02,4.49E+02,4.42E+02,3.82E+02,3.91E+02,2.91E+02,2.61E+02,
	1.32E+02,1.20E+02,1.16E+02,9.56E+01,9.94E+01,9.10E+01,4.85E+01,7.53E+01,9.53E+01,8.52E+01,5.76E+01,
	6.72E+01,5.20E+01,8.09E+00,3.75E+00,5.57E+00,3.80E+00,2.73E+00,2.22E+00,3.16E-01,1.26E-01,1.39E-01,
	6.15E-02,3.22E-02,7.98E-03};

/* stores numbers for table cooling flow */
#define	NCFL	40
static double cfle[NCFL], 
  cflf[NCFL];

/* Brad Peterson's AKN 120 continuum */
#define	NKN120	11
static double tnuakn[NKN120], 
  fnuakn[NKN120];

/* Black's ISM cotninuum, with He hole filled in */
#define	NISM	23
static double tnuism[NISM], 
  fnuism[NISM];

/* z=2 background,
 * >>refer	continuum	background	Haardt, Francesco, & Madau, Piero, 1996, 
 * >>refercon	ApJ, 461, 20 */
#define	NHM96	14
static double tnuHM96[NHM96]={-8,-1.722735683,-0.351545683,-0.222905683,-0.133385683,
/*                                    ================= mod these two to prevent degeneracy */
-0.127655683,-0.004575683,0.297544317,0.476753,0.476756,0.588704317,
0.661374317,1.500814317,2.245164317};
/*-0.127655683,-0.004575683,0.297544317,0.476754317,0.476754317,0.588704317,*/
/*log J in the units of (erg cm^{-2} s^{-1} Hz^{-1} sr^{-1})*/
static double fnuHM96[NHM96]={-32.53342863,-19.9789,-20.4204,-20.4443,-20.5756,-20.7546,
-21.2796,-21.6256,-21.8404,-21.4823,-22.2102,-22.9263,-23.32,-24.2865};

/* Mathews and Ferland generic agn continuum */
#define	NAGN	8
static double tnuagn[NAGN], 
  tslagn[NAGN];

#define NDRAINE	15
static double tnudrn[NDRAINE] , tsldrn[NDRAINE];

/* this is below and stores values for above vectors */
void ZeroContin(void);

void ParseTable(long int *nqh, 
  char *chCard,
  float *ar1)
{
	char chFile[FILENAME_PATH_LENGTH_2];	/*file name for table read */

	int intmode = 0,
	  lgEOL, 
	  lgHit, 
	  lgLogSet;
	static int lgCalled=FALSE;

	long int i, 
	  iabund=-1,
	  j, 
	  ncont,
	  nstar;

	double alpha, 
	  brakmm, 
	  brakxr, 
	  ConBreak, 
	  fac, 
	  par1, 
	  par2, 
	  scale, 
	  slopir, 
	  slopxr;

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

	/* if first call then set up values for table */
	if( !lgCalled )
	{
		ZeroContin();
		lgCalled = TRUE;
	}

	if( rfield.nspec >= LIMSPC )
	{
		fprintf( ioQQQ, " %ld is too many spectra entered.  Increase LIMSPC\n Sorry.\n", 
		  rfield.nspec );
		puts( "[Stop in ParseTable]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* set flag telling interpolate */
	strcpy( rfield.chSpType[rfield.nspec], "INTER" );

	/* when adding more keys also change the comment at the end */
	if( lgMatch(" AGN",chCard) )
	{
		/* do Mathews and Ferland generic AGN continuum */
		ASSERT( NAGN < NCELL);
		for( i=0; i < NAGN; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = tnuagn[i];
			rfield.tslop[i][rfield.nspec] = tslagn[i];
		}

		ncont = NAGN - 1;

		/* optional keywork, to adjust IR cutoff */
		if( lgMatch("BREA",chCard) )
		{
			i = 4;
			ConBreak = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

			if( lgEOL )
			{
				/* no break, set to low energy limit of code */
				if( lgMatch(" NO ",chCard) )
				{
					ConBreak = rfield.emm*1.01f;
				}
				else
				{
					fprintf( ioQQQ, " There must be a number for the break.\n Sorry.\n" );
					puts( "[Stop in ParseTable]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			if( ConBreak == 0. )
			{
				fprintf( ioQQQ, " The break must be greater than 0.2 Ryd.\n Sorry.\n" );
				puts( "[Stop in ParseTable]" );
				cdEXIT(EXIT_FAILURE);
			}

			if( lgMatch("MICR",chCard) )
			{
				/*  optional keyword, ``microns'', convert to Rydbergs */
				ConBreak = 0.0912/ConBreak;
			}

			if( ConBreak < 0. )
			{
				/*  option to enter break as LOG10 */
				ConBreak = pow(10.,ConBreak);
			}

			else if( ConBreak == 0. )
			{
				fprintf( ioQQQ, " An energy of 0 is not allowed.\n Sorry.\n" );
				puts( "[Stop in ParseTable]" );
				cdEXIT(EXIT_FAILURE);
			}

			if( ConBreak >= rfield.tNuRyd[2][rfield.nspec] )
			{
				fprintf( ioQQQ, " The energy of the break cannot be greater than%10.2e Ryd.\n Sorry.\n", 
				  rfield.tNuRyd[2][rfield.nspec] );
				puts( "[Stop in ParseTable]" );
				cdEXIT(EXIT_FAILURE);
			}

			else if( ConBreak <= rfield.tNuRyd[0][rfield.nspec] )
			{
				fprintf( ioQQQ, " The energy of the break cannot be less than%10.2e Ryd.\n Sorry.\n", 
				  rfield.tNuRyd[0][rfield.nspec] );
				puts( "[Stop in ParseTable]" );
				cdEXIT(EXIT_FAILURE);
			}

			rfield.tNuRyd[1][rfield.nspec] = ConBreak;

			rfield.tslop[1][rfield.nspec] = 
				(rfield.tslop[2][rfield.nspec] + 
			  log10(rfield.tNuRyd[2][rfield.nspec]/rfield.tNuRyd[1][rfield.nspec]));

			rfield.tslop[0][rfield.nspec] = 
				(rfield.tslop[1][rfield.nspec] - 
			  2.5*log10(rfield.tNuRyd[1][rfield.nspec]/rfield.tNuRyd[0][rfield.nspec]));
		}
	}

	else if( lgMatch("AKN1",chCard) )
	{
		/* AKN120 continuum from Brad Peterson's plot */
		ASSERT( NKN120 < NCELL );
		for( i=0; i < NKN120; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = tnuakn[i];
			rfield.tslop[i][rfield.nspec] = log10(fnuakn[i]);
		}
		ncont = NKN120 - 1;
	}

	else if( lgMatch("CRAB",chCard) )
	{
		ASSERT( NCRAB < NCELL );
		/* Crab Nebula continuum from Davidson and Fesen 1985 Ann Rev article */
		for( i=0; i < NCRAB; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = tnucrb[i];
			rfield.tslop[i][rfield.nspec] = log10(fnucrb[i]);
		}
		ncont = NCRAB - 1;
	}

	else if( lgMatch("COOL",chCard) )
	{
		ASSERT( NCFL < NCELL );
		/* cooling flow from Johnstone et al. 1992, MN 255, 431. */
		for( i=0; i < NCFL; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = cfle[i];
			rfield.tslop[i][rfield.nspec] = cflf[i];
		}
		ncont = NCFL - 1;
	}

	else if( (i=lgMatch("HM96",chCard))>0 )
	{
		/* this command does not include the cmb - do that separately with the fireball command */
#		if 0
		strcpy( rfield.chSpType[rfield.nspec], "BLACK" );
		/* >>chng 03 may 23, CMB temp from 2.756 to 2.725 */
		rfield.slope[rfield.nspec] = (CMB_TEMP*(1. + z));
		rfield.cutoff[0][rfield.nspec] = 0.;
		rfield.cutoff[1][rfield.nspec] = 0.;
		strcpy( rfield.chSpNorm[*nqh], "LUMI" );
		a = log10(rfield.slope[rfield.nspec]);
		rlogl = log10(2.99792e10*7.56464e-15) + 4.*a;
		strcpy( rfield.chRSpec[*nqh], "SQCM" );
		rfield.range[0][*nqh] = rfield.emm;
		rfield.range[1][*nqh] = rfield.egamry;
		rfield.totpow[*nqh] = rlogl;
		++rfield.nspec;
		++*nqh;
#		endif

		/* set flag telling interpolate */
		strcpy( rfield.chSpType[rfield.nspec], "INTER" );

		ASSERT( NHM96 < NCELL );
		/* z=2 background,
		* >>refer	continuum	background	Haardt, Francesco, & Madau, Piero, 1996, 
		* >>refercon	ApJ, 461, 20 */
		for( j=0; j < NHM96; j++ )
		{
			/* frequency was stored as log of ryd */
			rfield.tNuRyd[j][rfield.nspec] = pow( 10. , tnuHM96[j] );
			rfield.tslop[j][rfield.nspec] = fnuHM96[j];
		}

		ncont = NHM96 - 1;

		/* optional scale factor to change default luminosity
		 * from observed value
		 * assumed to be log if negative, and linear otherwise */
		/* note did not reset i, left after above 96 */
		i += 4;
		scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( scale > 0. )
			scale = log10(scale);

		/* this also sets continuum intensity*/
		if( *nqh >= LIMSPC )
		{
			fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 
			  *nqh );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* check that stack of shape and luminosity specifications
		 * is parallel, stop if not - this happens is background comes
		 * BETWEEN another set of shape and luminosity commands */
		if( rfield.nspec != *nqh )
		{
			fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		strcpy( rfield.chRSpec[*nqh], "SQCM" );
		strcpy( rfield.chSpNorm[*nqh], "FLUX" );

		/* this will be flux density at some frequency on the table.  the numbers
		 * are per Hz and sr so must multiply by 4 pi
		 * [2] is not special, could have been any within array*/
		rfield.range[0][*nqh] = pow(10. , tnuHM96[2] );

		/* interpolted from Black 1987 */
		rfield.totpow[*nqh] = (fnuHM96[2] + log10(PI4) + scale);

		/* set radius to very large value if not already set */
		/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
		if( !radius.lgRadiusKnown )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.f,radius.rdfalt);
		}
		++*nqh;

		/* set radius to very large value if not already set */
		/* set radius to very large value if not already set */
		/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
		if( !radius.lgRadiusKnown )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.f,radius.rdfalt);
		}
	}
	else if( lgMatch(" ISM",chCard) )
	{
		ASSERT( NISM < NCELL );
		/* local ISM radiation field from Black 1987, Interstellar Processes */
		/* >>chng 04 mar 16, rm cmb rom field so that it can be used at
		 * any redshift */
		rfield.tNuRyd[0][rfield.nspec] = 6.;
		rfield.tslop[0][rfield.nspec] = -21.21 - 6.;
		for( i=6; i < NISM; i++ )
		{
			rfield.tNuRyd[i-5][rfield.nspec] = tnuism[i];
			rfield.tslop[i-5][rfield.nspec] = (fnuism[i] - 
			  tnuism[i]);
		}

		ncont = NISM - 1 -5;
		i = 5;

		/* optional scale factor to change default luminosity
		 * from observed value
		 * want final number to be log 
		 * assumed to be log if negative, and linear otherwise unless log option is present */
		scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( scale > 0. && !lgMatch(" LOG",chCard) )
			scale = log10(scale);

		/* this also sets continuum intensity*/
		if( *nqh >= LIMSPC )
		{
			fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 
			  *nqh );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* check that stack of shape and luminosity specifications
		 * is parallel, stop if not - this happens is background comes
		 * BETWEEN another set of shape and luminosity commands */
		if( rfield.nspec != *nqh )
		{
			fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		strcpy( rfield.chRSpec[*nqh], "SQCM" );
		strcpy( rfield.chSpNorm[*nqh], "FLUX" );

		/* this will be flux density at 1 Ryd
		 * >>chng 96 dec 18, from 1 ryd to H mass rydberg
		 * >>chng 97 jan 10, had HLevNIonRyd but not defined yet */
		rfield.range[0][*nqh] = HIONPOT;

		/* interpolted from Black 1987 */
		rfield.totpow[*nqh] = (-18.517 + scale);

		/* set radius to very large value if not already set */
		/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
		if( !radius.lgRadiusKnown )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.f,radius.rdfalt);
		}
		*nqh += 1;

		if( optimize.lgVarOn )
		{
			optimize.nvarxt[optimize.nparm] = 1;
			strcpy( optimize.chVarFmt[optimize.nparm], "TABLE ISM LOG %f");
			/*  pointer to where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;
			/*  the scale factor */
			optimize.vparm[0][optimize.nparm] = (float)scale;
			optimize.vincr[optimize.nparm] = 0.2f;
			++optimize.nparm;
		}
	}
	else if( lgMatch("DRAI",chCard) )
	{
		ASSERT( NDRAINE < NCELL );
		/* local ISM radiation field from Draine 1978 */
		for( i=0; i < NDRAINE; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = tnudrn[i];
			rfield.tslop[i][rfield.nspec] = tsldrn[i];
		}

		ncont = NDRAINE - 1;
		i = 5;

		/* optional scale factor to change default luminosity
		 * from observed value
		 * assumed to be log if negative, and linear otherwise unless log option is present */
		scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( scale > 0. && !lgMatch(" LOG",chCard) )
			scale = log10(scale);

		/* this also sets continuum intensity*/
		if( *nqh >= LIMSPC )
		{
			fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 
			  *nqh );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* check that stack of shape and luminosity specifications
		 * is parallel, stop if not - this happens is background comes
		 * BETWEEN another set of shape and luminosity commands */
		if( rfield.nspec != *nqh )
		{
			fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
			fprintf( ioQQQ, " Sorry.\n" );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		strcpy( rfield.chRSpec[*nqh], "SQCM" );
		strcpy( rfield.chSpNorm[*nqh], "FLUX" );

		/* continuum normalization given by flux density at first point */
		rfield.range[0][*nqh] = tnudrn[0];

		/* this is f_nu at this first point */
		rfield.totpow[*nqh] = tsldrn[0] + scale;

		/* set radius to very large value if not already set */
		/* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
		if( !radius.lgRadiusKnown )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.f,radius.rdfalt);
		}
		*nqh += 1;
	}

	else if( lgMatch("LINE",chCard) )
	{
		static int lgLines_Malloc = FALSE;

		/* this is not a continuum source - it is to read a table of lines */
		if( lgLines_Malloc )
		{
			fprintf(ioQQQ," sorry, only one table line per input stream\n");
			cdEXIT(EXIT_FAILURE);
		}

		/* get file name within double quotes, if not present will use default */
		GetQuote( chLINE_LIST , chCard , FALSE );
		if( chLINE_LIST[0] == '\0' )
			strcpy( chLINE_LIST , "BLRLineList.dat" );

		lgLines_Malloc = TRUE;
		/* this flag says we did not find a continuum source - nearly all table XXX
		 * commands deal with continuum sources */
		ncont = -10;

		/* actually get the lines, and malloc the space in the arrays */
		nLINE_TABLE = cdGetLineList( chLINE_LIST , &chLabel , &wl);
	}

	else if( lgMatch("POWE",chCard) )
	{
		/* simple power law continuum between 10 mincron and 50 kev
		 *  option to read in any slope for the intermediate continuum */
		i = 5;
		alpha = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* default (no number on line) is f_nu propto nu^-1 */
		if( lgEOL )
			alpha = -1.;

		/* this is low energy for code */
		rfield.tNuRyd[0][rfield.nspec] = rfield.emm;
		/* and the value of the flux at this point (f_nu units)*/
		rfield.tslop[0][rfield.nspec] = -5;

		lgLogSet = FALSE;

		/* option to adjust sub-millimeter break */
		brakmm = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* default is 10 microns */
		if( lgEOL )
		{
			lgLogSet = TRUE;
			brakmm = 9.115e-3;
		}

		else if( brakmm == 0. )
		{
			/* if second number on line is zero then set lower limit to
			 * low-energy limit of the code.  Also set linear mode,
			 * so that last number will also be linear. */
			lgLogSet = FALSE;
			brakmm = rfield.tNuRyd[0][rfield.nspec]*1.001;
		}

		else if( brakmm < 0. )
		{
			/* if number is negative then this and next are logs */
			lgLogSet = TRUE;
			brakmm = pow(10.,brakmm);
		}

		/* optional microns keyword - convert to Rydbergs */
		if( lgMatch("MICR",chCard) )
			brakmm = RYDLAM / (1e4*brakmm) ;

		rfield.tNuRyd[1][rfield.nspec] = brakmm;

		/* check whether this is a reasonable mm break */
		if( brakmm > 1. )
			fprintf(ioQQQ,
			" Check the order of parameters on this table powerlaw command - the low-energy break of %f Ryd seems high to me.\n",
			brakmm );

		/* this is spectral slope, in F_nu units, between the low energy limit 
		 * and the break that may have been adjusted above 
		 * this is the slope appropriate for self-absorbed synchrotron, see eq 6.54, p.190
		 *>>refer	continuum	synchrotron	Rybicki, G. B., & Lightman, A.P. 1979, 
		 *>>refercon	Radiative Processes in Astrophysics (New York: Wiley)*/
		slopir = 5./2.;

		/* now extrapolate a flux at this energy using the flux entered for
		 * the first point, and this slope */
		rfield.tslop[1][rfield.nspec] = 
		  (rfield.tslop[0][rfield.nspec] + 
		  slopir*log10(rfield.tNuRyd[1][rfield.nspec]/rfield.tNuRyd[0][rfield.nspec]));

		/* this is energy of the high energy limit to code */
		rfield.tNuRyd[3][rfield.nspec] = rfield.egamry;

		/* option to adjust hard X-ray break */
		brakxr = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* default is 50 keV */
		if( lgEOL )
		{
			brakxr = 3676.;
		}

		else if( brakxr == 0. )
		{
			brakxr = rfield.tNuRyd[3][rfield.nspec]/1.001;
		}

		else if( lgLogSet )
		{
			/* first number was negative this is a logs */
			brakxr = pow(10.,brakxr);
		}

		/* note that this second cutoff does not have the micron keyword */
		rfield.tNuRyd[2][rfield.nspec] = brakxr;

		/* >>chng 03 jul 19, check that upper energy is greater than lower energy,
		 * quit if this is not the case */
		if( brakmm >= brakxr )
		{
			fprintf( ioQQQ, " HELP!!  The lower energy for the power law, %f, is greater than the upper energy, %f. This is not possible.\n Sorry.\n",
				brakmm , brakxr );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* alpha was first option on line, is slope of mid-range */
		rfield.tslop[2][rfield.nspec] = 
		  (rfield.tslop[1][rfield.nspec] + 
		  alpha*log10(rfield.tNuRyd[2][rfield.nspec]/rfield.tNuRyd[1][rfield.nspec]));

		/* high energy range is nu^-2 */
		slopxr = -2.;

		rfield.tslop[3][rfield.nspec] = 
			(rfield.tslop[2][rfield.nspec] + 
		  slopxr*log10(rfield.tNuRyd[3][rfield.nspec]/rfield.tNuRyd[2][rfield.nspec]));

		/* following is number of portions of continuum */
		ncont = 3;
	}

	else if( lgMatch("READ",chCard) )
	{
		/* make sure not more than one table read command */
		if( rfield.lgTableRead )
		{
			fprintf(ioQQQ," Oops, there are more than one TABLE READ (or TLUSTY) command in this input stream.  I can only deal with a single TABLE READ.\n Sorry.\n");
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}
		rfield.lgTableRead = TRUE;

		/* set up eventual read of table of points previously punched by code 
		 * get file name within double quotes, return as null terminated string
		 * also blank out original, chCard  version of name so do not trigger */
		GetQuote( chFile , chCard , TRUE );

		/* now open this file */
		if( NULL==(rfield.ioTableRead[rfield.nspec]=fopen( chFile,"r")))
		{
			fprintf( ioQQQ, " I could not open file %s.\n",chFile );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* set flag saying really just read in continuum exactly as punched */
		strcpy( rfield.chSpType[rfield.nspec], "READ " );

		/* this will be flag saying continuum not read in here */
		rfield.tslop[0][rfield.nspec] = 0.;
		rfield.tslop[1][rfield.nspec] = 0.;
		/* nothing left to do here, so return
		 * the file will actually read in by routine ReadTable (in ffun.c)
		 * called by ffun when continuum is being defined */

		/* number of spectra shapes that have been specified */
		++rfield.nspec;
		
#		ifdef DEBUG_FUN
		fputs( " <->ParseTable()\n", debug_fp );
#		endif
		return;

	}
	else if( lgMatch("TLUSTY",chCard) )
	{
		/* make sure not more than one table read command */
		if( rfield.lgTableRead )
		{
			fprintf(ioQQQ," Oops, there are more than one TABLE READ (or TLUSTY) command in this input stream.  I can only deal with a single TABLE READ.\n Sorry.\n");
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}
		rfield.lgTableRead = TRUE;

		/* set up eventual read of a TLUSTY spectrum file (added 5 May 2003 K Volk)
		 * get file name within double quotes, return as null terminated string
		 * also blank out original, chCard  version of name so do not trigger */
		GetQuote( chFile , chCard , TRUE );

		/* now open this file */
		if( NULL==(rfield.ioTableRead[rfield.nspec]=fopen( chFile,"r")))
		{
			fprintf( ioQQQ, " I could not open file %s.\n",chFile );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* set flag saying read TLUSTY file */
		strcpy( rfield.chSpType[rfield.nspec], "TLUST" );

		/* this will be flag saying continuum not read in here */
		rfield.tslop[0][rfield.nspec] = 0.;
		rfield.tslop[1][rfield.nspec] = 0.;
		/* nothing left to do here, so return
		 * the file will actually read in by routine ReadTlusty
		 * called by ffun when continuum is being defined */

		/* number of spectra shapes that have been specified */
		++rfield.nspec;
		
#		ifdef DEBUG_FUN
		fputs( " <->ParseTable()\n", debug_fp );
#		endif
		return;

	}
	else if( lgMatch("RUBI",chCard) )
	{
		/* do Rubin modified theta 1 ori c */
		for( i=0; i < NRUBIN; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = tnurbn[i];
			rfield.tslop[i][rfield.nspec] = log10(fnurbn[i] /tnurbn[i] );
		}
		ncont = NRUBIN - 1;
	}

	else if( lgMatch("STAR",chCard) && lgMatch("BURST",chCard) )
	{
		/* set up eventual read of starburst 99 predictions 
		 * get file name within double quotes, return as null terminated string
		 * also blank out original, chCard  version of name so do not trigger */
		GetQuote( chFile , chCard , TRUE );

		/* now open this file */
		if( NULL==(rfield.ioTableRead[rfield.nspec]=fopen( chFile,"r")))
		{
			fprintf( ioQQQ, " I could not open file %s.\n",chFile );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* set flag saying want to read in output from starburst 99 code */
		strcpy( rfield.chSpType[rfield.nspec], "SB99 " );

		/* this will be flag saying continuum not read in here */
		rfield.tslop[0][rfield.nspec] = 0.;
		rfield.tslop[1][rfield.nspec] = 0.;

		/* get linear age of the starburst in years */
		i = 5;
		rfield.tslop[0][rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* keyword log forces interpretation as log of age */
		if( lgMatch(" LOG",chCard) )
		{
			rfield.tslop[0][rfield.nspec] = pow(10., rfield.tslop[0][rfield.nspec] );
		}

		if( optimize.lgVarOn )
		{
			optimize.nvarxt[optimize.nparm] = 1;
			strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STARBURST %f");
			sprintf( optimize.chVarFmt[optimize.nparm], "TABLE STARBURST %%f LOG AGE FILE=\"%s\" ", 
			  chFile );
			/*  pointer to where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;
			/*  log of age */
			optimize.vparm[0][optimize.nparm] = (float)log10(rfield.tslop[0][rfield.nspec]);
			optimize.vincr[optimize.nparm] = 0.3f;
			++optimize.nparm;
		}

		/* nothing left to do here, so return
		 * the file will actually read in by routine ReadTable
		 * called by ffun when continuum is being defined */

		/* number of spectra shapes that have been specified */
		++rfield.nspec;
		
#		ifdef DEBUG_FUN
		fputs( " <->ParseTable()\n", debug_fp );
#		endif
		return;

	}

	else if( lgMatch("STAR",chCard) )
	{
		/* TABLE STAR, numbers are temperature and log(g) */
		i = 5;
		par1 = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* option for log if less than or equal to 10 */
		/* Caution: For CoStar models this implicitly assumes that
		 * the minimum ZAMS mass and youngest age is greater than 10.
		 * As of June 1999 this is the case. However, this is not
		 * guaranteed for possible future upgrades. */
		if( par1 <= 10. )
			par1 = pow(10.,par1);

		if( lgEOL )
		{
			fprintf( ioQQQ, " The stellar temperature MUST be entered.\n" );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* lgEOL is used in following logic to check whether grav set */
		par2 = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgMatch("ATLA",chCard) )
		{
			/* this sub-branch added by Kevin Volk, to read in large
			 * 1991 grid of Kurucz atmospheres
			 * default log(g) is 5 since all T_eff values have a model
			 * for this log(g) */
			if( lgEOL )
				par2 = 5.0;

			if( lgMatch("TRAC",chCard) )
			{
				lgEOL = TRUE;
			}
			else
			{
				lgEOL = FALSE;
			}

			AtlasInterpolate(&nstar,par1,par2,lgEOL);
			/* set flag saying really just read in continuum exactly as punched */
			strcpy( rfield.chSpType[rfield.nspec], "VOLK " );
			ncont = nstar - 1;
		}

		else if( lgMatch("WERN",chCard) )
		{
			/* this subbranch added by Kevin Volk, to read in large
			 * grid of hot PN atmospheres computer by Klaus Werner.
			 * default log(g) is 8 since all T_eff values have a model
			 * for this log(g) */
			if( lgEOL )
				par2 = 8.0;

			WernerInterpolate(&nstar,par1,par2);

			/* set flag saying really just read in continuum exactly as punched */
			strcpy( rfield.chSpType[rfield.nspec], "VOLK " );

			ncont = nstar - 1;
			/* >>chng 97 aug 23, K Volk added Rauch stellar atmospheres */
		}

		else if( lgMatch("COST",chCard) )
		{
			int nmodid = INT_MAX;

			/* there are two abundance sets, solar and halo.  If halo keyword
			 * appears use halo, else use solar */
			if( lgMatch("HALO",chCard) )
			{
				iabund = 1;
			}
			else
			{
				iabund = 0;
			}

			/* >>chng 99 apr 30, added CoStar stellar atmospheres */
			/* this subbranch reads in CoStar atmospheres, no log(g),
			 * but second parameter is age sequence, a number between 1 and 7,
			 * default is 1 */
			if( lgMatch("INDE",chCard) )
			{
				/* this is Teff, index mode */
				intmode = 1;
				if( lgEOL )
					par2 = 1.;
				/* now make sure that second parameter is within allowed range -
				 * we do not have enough information at this time to verify temperature */
				nmodid = (int)par2;
				if( nmodid < 1 || nmodid > 7 )
				{
					fprintf( ioQQQ, " The second number must be the id sequence number, 1 to 7.\n" );
					fprintf( ioQQQ, " reset to 1.\n" );
					par2 = 1.;
					nmodid = 1;
				}
			}
			else if( lgMatch("ZAMS",chCard) )
			{
				/* this is M_ZAMS, age mode */
				intmode = 3;
				if( lgEOL )
				{
					fprintf( ioQQQ, " A second number, the age of the star, must be present.\n" );
					puts( "[Stop in ParseTable]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else if( lgMatch(" AGE",chCard) )
			{
				/* this is age, M_ZAMS mode */
				intmode = 4;
				if( lgEOL )
				{
					fprintf( ioQQQ, " A second number, the ZAMS mass of the star, must be present.\n" );
					puts( "[Stop in ParseTable]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else
			{
				/* this is Teff, log(g) mode */
				/* Caution: the code in CoStarInterpolate implicitly assumes that the dependence
				 * of log(g) with age is strictly monotomic. As of June 1999 this is the case. */
				intmode = 2;
				/* default is to use ZAMS models, i.e. use index 1 */
				if( lgEOL )
				{
					intmode = 1;
					par2 = 1.;
					nmodid = 1;
				}
			}
			CoStarInterpolate(&nstar,intmode,par1,par2,nmodid,(int)iabund);

			/* set flag saying really just read in continuum exactly as punched */
			strcpy( rfield.chSpType[rfield.nspec], "VOLK " );

			ncont = nstar - 1;
		}

		else if( lgMatch("RAUC",chCard) )
		{
			/* default surface gravity for these very hot stars is 8.0 */
			if( lgEOL )
			{
				par2 = 8.0;
			}

			/* there are two abundance sets, solar and halo.  If halo keyword
			 * appears use halo, else use solar */
			if( lgMatch("HALO",chCard) )
			{
				iabund = 1;
			}
			else
			{
				iabund = 0;
			}

			/* go get the atmosphere */
			if( lgMatch(" OLD",chCard) )
			{
				RauchInterpolate(&nstar,par1,par2,iabund);
			}
			else
			{
				RauchInterpolate02(&nstar,par1,par2,iabund);
			}

			/* set flag saying really just read in continuum exactly as punched */
			strcpy( rfield.chSpType[rfield.nspec], "VOLK " );

			ncont = nstar - 1;
			/* End of block (K Volk) */
		}

		else if( lgMatch("KURU",chCard) )
		{
			/* only log(g) = 4 supported now in Kurucz */
			if( lgEOL )
				par2 = 4.;

			kurucz79(&nstar,par1,par2);
			ncont = nstar - 1;
		}

		else if( lgMatch("MIHA",chCard) )
		{
			/* only log(g) = 4 supported now in Mihalas */
			if( lgEOL )
				par2 = 4.;

			mihals(&nstar,par1,par2);
			ncont = nstar - 1;
		}

		else
		{
			fprintf( ioQQQ, " There must be a second key on TABLE STAR;  the options are ATLAS, KURUCZ, MIHALAS, RAUCH, and WERNER.\n" );
			puts( "[Stop in ParseTable]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* vary option */
		if( optimize.lgVarOn )
		{
			optimize.vincr[optimize.nparm] = (float)0.1;

			if( lgMatch("ATLA",chCard) )
			{
				/* this is number of parameters to feed onto the input line */
				optimize.nvarxt[optimize.nparm] = 2;
				strcpy( optimize.chVarFmt[optimize.nparm], 
					"TABLE STAR ATLAS %f log(g)=%f" );

				/* following are upper and lower limits to temperature range */
				optimize.varang[optimize.nparm][0] = (float)log10(3500.);
				optimize.varang[optimize.nparm][1] = (float)log10(50000.);
			}

			else if( lgMatch("WERN",chCard) )
			{
				/* this is number of parameters to feed onto the input line */
				optimize.nvarxt[optimize.nparm] = 2;
				strcpy( optimize.chVarFmt[optimize.nparm], 
					"TABLE STAR WERNER %f log(g)=%f" );

				/* following are upper and lower limits to temperature range */
				optimize.varang[optimize.nparm][0] = (float)log10(80000.);
				optimize.varang[optimize.nparm][1] = (float)log10(200000.);
			}

			else if( lgMatch("RAUC",chCard) )
			{
				/* this is number of parameters to feed onto the input line */
				optimize.nvarxt[optimize.nparm] = 2;

				/* iabund was set when command read in, for different abunances.
				 * 0 is for solar, and 1 is for halo */
				ASSERT( iabund==0 || iabund==1 );
				if( iabund )
				{
					/* the halo case */
					strcpy( optimize.chVarFmt[optimize.nparm], 
						"TABLE STAR RAUCH HALO %f log(g)=%f" );
				}
				else
				{
					/* the solar case */
					strcpy( optimize.chVarFmt[optimize.nparm], 
						"TABLE STAR RAUCH SOLAR %f log(g)=%f" );
				}

				/* following are upper and lower limits to temperature range */
				optimize.varang[optimize.nparm][0] = (float)log10(50000.);
				optimize.varang[optimize.nparm][1] = (float)log10(500000.);
			}

			else if( lgMatch("COST",chCard) )
			{
				/* this is number of parameters to feed onto the input line */
				optimize.nvarxt[optimize.nparm] = 2;
				if( intmode == 1 )
				{
					strcpy( optimize.chVarFmt[optimize.nparm], 
						"TABLE STAR COSTAR TEFF= %f INDEX= %f" );
				}
				else if( intmode == 2 )
				{
					strcpy( optimize.chVarFmt[optimize.nparm], 
						"TABLE STAR COSTAR TEFF= %f LOG(G)= %f" );
				}
				else if( intmode == 3 )
				{
					strcpy( optimize.chVarFmt[optimize.nparm], 
						"TABLE STAR COSTAR ZAMS MASS= %f TIME= %f" );
				}
				else if( intmode == 4 )
				{
					strcpy( optimize.chVarFmt[optimize.nparm], 
						"TABLE STAR COSTAR AGE= %f MASS= %f" );
					optimize.vincr[optimize.nparm] = (float)0.5;
				}
				/* the limits for the parameter to be varied have already been set in CoStarInterpolate */
			}

			else if( lgMatch("KURU",chCard) )
			{
				/* this is number of parameters to feed onto the input line */
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], 
					"TABLE STAR KURUCZ %f" );

				/* following are upper and lower limits to temperature range */
				optimize.varang[optimize.nparm][0] = (float)log10(30000.);
				optimize.varang[optimize.nparm][1] = (float)log10(50000.);
			}

			else if( lgMatch("MIHA",chCard) )
			{
				/* this is number of parameters to feed onto the input line */
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], 
					"TABLE STAR MIHALAS %f" );

				/* following are upper and lower limits to temperature range */
				optimize.varang[optimize.nparm][0] = (float)log10(30000.);
				optimize.varang[optimize.nparm][1] = (float)log10(55000.);
			}

			/* pointer to where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;

			/* log of temp will be pointer */
			optimize.vparm[0][optimize.nparm] = (float)log10(par1);
			optimize.vparm[1][optimize.nparm] = (float)par2;

			/* finally increment this counter */
			++optimize.nparm;
		}
	}

	else
	{
		fprintf( ioQQQ, " There MUST be a keyword on this line. The keys are: AGN, AKN120, _ISM, COOL, CRAB, POWERlaw, RUBIN, STAR and TLUSTY.  Sorry.\n" );
		puts( "[Stop in ParseTable]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* table star werner and table star atlas are special
	 * cases put in by Kevin Volk - they are not really tables
	 * at all, so return if sptype is Volk */
	if( strcmp(rfield.chSpType[rfield.nspec],"VOLK ") == 0 )
	{ 
		++rfield.nspec;
#		ifdef DEBUG_FUN
		fputs( " <->ParseTable()\n", debug_fp );
#		endif
		return;
	}

	/* ncont only positive when real continuum entered */
	/* convert from log(Hz) to Ryd if first nu>5 */
	if( rfield.tNuRyd[0][rfield.nspec] >= 5. )
	{
		for( i=0; i < (ncont + 1); i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = 
				pow(10.,rfield.tNuRyd[i][rfield.nspec] - 15.51718);
		}
	}

	else if( rfield.tNuRyd[0][rfield.nspec] < 0. )
	{
		/* energies entered as logs */
		for( i=0; i < (ncont + 1); i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = 
				pow(10.,rfield.tNuRyd[i][rfield.nspec]);
		}
	}

	/* tfac will be log(fnu) at that spot, read into tslop */
	for( i=0; i < (ncont + 1); i++ )
	{
		rfield.tFluxLog[i][rfield.nspec] = rfield.tslop[i][rfield.nspec];
	}

	for( i=0; i < ncont; i++ )
	{
		rfield.tslop[i][rfield.nspec] = 
			((rfield.tslop[i+1][rfield.nspec] - 
		  rfield.tslop[i][rfield.nspec])/
		  log10(rfield.tNuRyd[i+1][rfield.nspec]/
		  rfield.tNuRyd[i][rfield.nspec]));
	}

	if( ncont > 0 && (ncont + 2 < rfield.nupper) )
	{
		/* zero out remainder of array */
		for( i=ncont + 1; i < rfield.nupper; i++ )
		{
			rfield.tNuRyd[i][rfield.nspec] = 0.;
		}
	}

	if( trace.lgConBug && trace.lgTrace )
	{
		fprintf( ioQQQ, " Table for this continuum; TNU, TFAC, TSLOP\n" );
		for( i=0; i < (ncont + 1); i++ )
		{
			fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNuRyd[i][rfield.nspec], 
			  rfield.tFluxLog[i][rfield.nspec], rfield.tslop[i][rfield.nspec] );
		}
	}

	/* renormalize continuum so that quantity passes through unity
	 * at 1 Ryd */
	/* lgHit will be set false when we get a hit in following loop,
	 * then test made to make sure it happened*/
	lgHit = FALSE;
	/*following will be reset when hit occurs*/
	fac = -DBL_MAX;
	/* >>chng 04 mar 16, chn loop so breaks when hit, previously had cycled
	 * until last point reached, so last pnt always used */
	for( i=0; i < (ncont + 1) && !lgHit; i++ )
	{
		if( rfield.tNuRyd[i][rfield.nspec] > 1. )
		{
			fac = rfield.tFluxLog[i][rfield.nspec];
			lgHit = TRUE;
		}
	}

	if( ncont > 0 && lgHit )
	{
		/* do the renormalization */
		for( i=0; i < (ncont + 1); i++ )
		{
			rfield.tFluxLog[i][rfield.nspec] -= fac;
		}
	}

	if( ncont > 0 )
		++rfield.nspec;

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

/*ZeroContin store sets of built in continua */

/* this allows the low energy point of any built in array to be reset to the
 * current low energy point in the code - nothing need be done if this is reset 
 * tnu is array of energies, [0] is first, and we want it to be lower
 * fluxog is flux at tnu, and may or may not be log 
 * lgLog says whether it is */
void resetBltin( double *tnu , double *fluxlog , int lgLog )
{
	/* this will mult low energy bounds of code and go into element[0] 
	 * ensures that energy range is fully covered */
	const double RESETFACTOR = 0.98 ;
	double power;
	/* this makes sure we are called after emm is defined */
	ASSERT( rfield.emm  > 0. );

	if( lgLog )
	{
		/* continuum comes in as log of flux */
		/* this is current power-law slope of low-energy continuum */
		power = (fluxlog[1] - fluxlog[0] ) / log10( tnu[1]/tnu[0] );
		/* this will be new low energy bounds to this continuum */
		tnu[0] = rfield.emm*RESETFACTOR;
		fluxlog[0] = fluxlog[1] + power * log10( tnu[0] /tnu[1] );
	}
	else
	{
		/* continuum comes in as linear flux */
		/* this is current power-law slope of low-energy continuum */
		power = log10( fluxlog[1]/fluxlog[0]) / log10( tnu[1]/tnu[0] );
		/* this will be new low energy bounds to this continuum */
		tnu[0] = rfield.emm*RESETFACTOR;
		fluxlog[0] = log10(fluxlog[1]) + power * log10( tnu[0] /tnu[1] );
		/* flux is not really log, we want linear */
		fluxlog[0] = pow(10. , fluxlog[0]);
	}
	/*fprintf(ioQQQ," power is %f lgLog is %i\n", power, lgLog );*/
	return;
}

void ZeroContin(void)
{

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

	/* Draine 1978 ISM continuum */
	/* freq in ryd */
	tnudrn[0] = 0.3676;
	tnudrn[1] = 0.4144;
	tnudrn[2] = 0.4558;
	tnudrn[3] = 0.5064;
	tnudrn[4] = 0.5698;
	tnudrn[5] = 0.6511;
	tnudrn[6] = 0.7012;
	tnudrn[7] = 0.7597;
	tnudrn[8] = 0.8220;
	tnudrn[9] = 0.9116;
	tnudrn[10] = 0.9120;
	tnudrn[11] = 0.9306;
	tnudrn[12] = 0.9600;
	tnudrn[13] = 0.9806;
	tnudrn[14] = 0.9999;

	/* f_nu from equation 23 of 
	 * >>refer	ism	field	Draine, B.T. & Bertoldi, F. 1996, ApJ, 468, 269 */
	{
		int i;
		i= 0;
		tsldrn[i] = -17.8063;
		++i;tsldrn[i] = -17.7575;
		++i;tsldrn[i] = -17.7268;
		++i;tsldrn[i] = -17.7036;
		++i;tsldrn[i] = -17.6953;
		++i;tsldrn[i] = -17.7182;
		++i;tsldrn[i] = -17.7524;
		++i;tsldrn[i] = -17.8154;
		++i;tsldrn[i] = -17.9176;
		++i;tsldrn[i] = -18.1675;
		++i;tsldrn[i] = -18.1690;
		++i;tsldrn[i] = -18.2477;
		++i;tsldrn[i] = -18.4075;
		++i;tsldrn[i] = -18.5624;
		++i;tsldrn[i] = -18.7722;
	}

#	if 0
	tsldrn[0] = -15.8344;
	tsldrn[1] = -16.0454;
	tsldrn[2] = -16.2217;
	tsldrn[3] = -16.4273;
	tsldrn[4] = -16.6747;
	tsldrn[5] = -16.9875;
	tsldrn[6] = -17.1827;
	tsldrn[7] = -17.4195;
	tsldrn[8] = -17.6929;
	tsldrn[9] = -18.1675;
	tsldrn[10] = -18.1700;
	tsldrn[11] = -18.2925;
	tsldrn[12] = -18.5199;
	tsldrn[13] = -18.7209;
	tsldrn[14] = -18.9730;
#	endif

	/* generic AGN continuum taken from 
	 * >>refer	AGN	cont	Mathews and Ferland ApJ Dec15 '87
	 * except self absorption at 10 microns, nu**-2.5 below that */
	tnuagn[0] = 1e-5;
	tnuagn[1] = 9.12e-3;
	tnuagn[2] = .206;
	tnuagn[3] = 1.743;
	tnuagn[4] = 4.13;
	tnuagn[5] = 26.84;
	tnuagn[6] = 7353.;
	tnuagn[7] = 7.4e6;

	tslagn[0] = -3.388;
	tslagn[1] = 4.0115;
	tslagn[2] = 2.6576;
	tslagn[3] = 2.194;
	tslagn[4] = 1.819;
	tslagn[5] = -.6192;
	tslagn[6] = -2.326;
	tslagn[7] = -7.34;
	resetBltin( tnuagn , tslagn , TRUE );

	/* Crab Nebula continuum from 
	 *>>refer	continuum	CrabNebula	Davidson, K., and Fesen, 1985, ARAA,  
	 * second vector is F_nu as seen from Earth */
	tnucrb[0] = 1.0e-5;
	tnucrb[1] = 5.2e-4;
	tnucrb[2] = 1.5e-3;
	tnucrb[3] = 0.11;
	tnucrb[4] = 0.73;
	tnucrb[5] = 7.3;
	tnucrb[6] = 73.;
	tnucrb[7] = 7300.;
	tnucrb[8] = 1.5e6;
	tnucrb[9] = 7.4e6;

	fnucrb[0] = 3.77e-21;
	fnucrb[1] = 1.38e-21;
	fnucrb[2] = 2.10e-21;
	fnucrb[3] = 4.92e-23;
	fnucrb[4] = 1.90e-23;
	fnucrb[5] = 2.24e-24;
	fnucrb[6] = 6.42e-26;
	fnucrb[7] = 4.02e-28;
	fnucrb[8] = 2.08e-31;
	fnucrb[9] = 1.66e-32;
	resetBltin( tnucrb , fnucrb , FALSE );

	/* Bob Rubin's theta 1 Ori C continuum, modified from Kurucz to
	 * get NeIII right */
	/* >>chng 02 may 02, revise continuum while working with Bob Rubin on NII revisit */
#	if 0
	tnurbn[0] = 1e-5;
	tnurbn[1] = 1.0;
	tnurbn[2] = 1.02445;
	tnurbn[3] = 1.07266;
	tnurbn[4] = 1.12563;
	tnurbn[5] = 1.18411;
	tnurbn[6] = 1.23881;
	tnurbn[7] = 1.29328;
	tnurbn[8] = 1.35881;
	tnurbn[9] = 1.42463;
	tnurbn[10] = 1.48981;
	tnurbn[11] = 1.55326;
	tnurbn[12] = 1.6166;
	tnurbn[13] = 1.68845;
	tnurbn[14] = 1.76698;
	tnurbn[15] = 1.8019;
	tnurbn[16] = 1.8078;
	tnurbn[17] = 1.8082;
	tnurbn[18] = 1.84567;
	tnurbn[19] = 1.9317;
	tnurbn[20] = 2.04891;
	tnurbn[21] = 2.14533;
	tnurbn[22] = 2.19702;
	tnurbn[23] = 2.27941;
	tnurbn[24] = 2.37438;
	tnurbn[25] = 2.43137;
	tnurbn[26] = 2.49798;
	tnurbn[27] = 2.56113;
	tnurbn[28] = 2.59762;
	tnurbn[29] = 2.66597;
	tnurbn[30] = 2.80543;
	tnurbn[31] = 2.95069;
	tnurbn[32] = 3.02911;
	tnurbn[33] = 3.11182;
	tnurbn[34] = 3.22178;
	tnurbn[35] = 3.3155;
	tnurbn[36] = 3.42768;
	tnurbn[37] = 3.50678;
	tnurbn[38] = 3.56157;
	tnurbn[39] = 3.61811;
	tnurbn[40] = 3.75211;
	tnurbn[41] = 3.89643;
	tnurbn[42] = 4.0;
	tnurbn[43] = 4.1;
	tnurbn[44] = 1e5;
	/* Emergent fluxes are in units of erg/(cm^2 s Hz) and are what Kurucz 
	 * (and model atmosphere types) call pi*F_nu.    */
	fnurbn[0] = .181419e-11;
	fnurbn[1] = .181419e-01;
	fnurbn[2] = .145003e-01;
	fnurbn[3] = .129057e-01;
	fnurbn[4] = .130125e-01;
	fnurbn[5] = .115540e-01;
	fnurbn[6] = .114839e-01;
	fnurbn[7] = .964582e-02;
	fnurbn[8] = .922259e-02;
	fnurbn[9] = .829783e-02;
	fnurbn[10] = .735761e-02;
	fnurbn[11] = .651893e-02;
	fnurbn[12] = .617587e-02;
	fnurbn[13] = .510496e-02;
	fnurbn[14] = .499538e-02;
	fnurbn[15] = .364324e-02;
	fnurbn[16] = .325440e-02;
	fnurbn[17] = .325440e-02;
	fnurbn[18] = .162006e-02;
	fnurbn[19] = .139964e-02;
	fnurbn[20] = .127473e-02;
	fnurbn[21] = .100621e-02;
	fnurbn[22] = .102095e-02;
	fnurbn[23] = .900858e-03;
	fnurbn[24] = .460721e-03;
	fnurbn[25] = .699557e-03;
	fnurbn[26] = .861525e-03;
	fnurbn[27] = .751004e-03;
	fnurbn[28] = .500242e-03;
	fnurbn[29] = .568666e-03;
	fnurbn[30] = .418121e-03;
	fnurbn[31] = .619107e-04;
	fnurbn[32] = .279322e-04;
	fnurbn[33] = .404350e-04;
	fnurbn[34] = .266258e-04;
	fnurbn[35] = .186196e-04;
	fnurbn[36] = .146275e-04;
	fnurbn[37] = .203668e-05;
	fnurbn[38] = .799454e-06;
	fnurbn[39] = .868859e-06;
	fnurbn[40] = .370180e-06;
	fnurbn[41] = .186293e-06;
	fnurbn[42] = .450348e-07;
	fnurbn[43] = 1e-15;
	fnurbn[44] = 1e-20;
#	endif
	resetBltin( tnurbn , fnurbn , FALSE );

	/* cooling flow continuum from Johnstone et al. MNRAS 255, 431. */
	cfle[0] = 0.0000100;
	cflf[0] = -0.8046910;
	cfle[1] = 0.7354023;
	cflf[1] = -0.8046910;
	cfle[2] = 1.4708046;
	cflf[2] = -0.7436830;
	cfle[3] = 2.2062068;
	cflf[3] = -0.6818757;
	cfle[4] = 2.9416091;
	cflf[4] = -0.7168990;
	cfle[5] = 3.6770115;
	cflf[5] = -0.8068384;
	cfle[6] = 4.4124136;
	cflf[6] = -0.6722584;
	cfle[7] = 5.1478162;
	cflf[7] = -0.7626385;
	cfle[8] = 5.8832183;
	cflf[8] = -1.0396487;
	cfle[9] = 6.6186204;
	cflf[9] = -0.7972314;
	cfle[10] = 7.3540230;
	cflf[10] = -0.9883416;
	cfle[11] = 14.7080460;
	cflf[11] = -1.1675659;
	cfle[12] = 22.0620689;
	cflf[12] = -1.1985949;
	cfle[13] = 29.4160919;
	cflf[13] = -1.2263466;
	cfle[14] = 36.7701149;
	cflf[14] = -1.2918345;
	cfle[15] = 44.1241379;
	cflf[15] = -1.3510833;
	cfle[16] = 51.4781609;
	cflf[16] = -1.2715496;
	cfle[17] = 58.8321838;
	cflf[17] = -1.1098027;
	cfle[18] = 66.1862030;
	cflf[18] = -1.4315782;
	cfle[19] = 73.5402298;
	cflf[19] = -1.1327956;
	cfle[20] = 147.080459;
	cflf[20] = -1.6869649;
	cfle[21] = 220.620681;
	cflf[21] = -2.0239367;
	cfle[22] = 294.160919;
	cflf[22] = -2.2130392;
	cfle[23] = 367.701141;
	cflf[23] = -2.3773901;
	cfle[24] = 441.241363;
	cflf[24] = -2.5326197;
	cfle[25] = 514.7816162;
	cflf[25] = -2.5292389;
	cfle[26] = 588.3218384;
	cflf[26] = -2.8230250;
	cfle[27] = 661.8620605;
	cflf[27] = -2.9502323;
	cfle[28] = 735.4022827;
	cflf[28] = -3.0774822;
	cfle[29] = 1470.8045654;
	cflf[29] = -4.2239799;
	cfle[30] = 2206.2067871;
	cflf[30] = -5.2547927;
	cfle[31] = 2941.6091309;
	cflf[31] = -6.2353640;
	cfle[32] = 3677.0114746;
	cflf[32] = -7.1898708;
	cfle[33] = 4412.4135742;
	cflf[33] = -8.1292381;
	cfle[34] = 5147.8159180;
	cflf[34] = -9.0594845;
	cfle[35] = 5883.2182617;
	cflf[35] = -9.9830370;
	cfle[36] = 6618.6206055;
	cflf[36] = -10.9028034;
	cfle[37] = 7354.0229492;
	cflf[37] = -11.8188877;
	cfle[38] = 7400.0000000;
	cflf[38] = -30.0000000;
	cfle[39] = 10000000.0000000;
	cflf[39] = -30.0000000;
	resetBltin( cfle , cflf , TRUE );

	/* AKN120 continuum from Brad Peterson's plot
	 * second vector is F_nu*1E10 as seen from Earth */
	tnuakn[0] = 1e-5;
	tnuakn[1] = 1.9e-5;
	tnuakn[2] = 3.0e-4;
	tnuakn[3] = 2.4e-2;
	tnuakn[4] = 0.15;
	tnuakn[5] = 0.30;
	tnuakn[6] = 0.76;
	tnuakn[7] = 2.0;
	tnuakn[8] = 76.0;
	tnuakn[9] = 760.;
	tnuakn[10] = 7.4e6;
	fnuakn[0] = 1.5e-16;
	fnuakn[1] = 1.6e-16;
	fnuakn[2] = 1.4e-13;
	fnuakn[3] = 8.0e-15;
	fnuakn[4] = 1.6e-15;
	fnuakn[5] = 1.8e-15;
	fnuakn[6] = 7.1e-16;
	fnuakn[7] = 7.9e-17;
	fnuakn[8] = 1.1e-18;
	fnuakn[9] = 7.1e-20;
	fnuakn[10] = 1.3e-24;
	resetBltin( fnuakn , fnuakn , FALSE );

	/* interstellar radiation field from Black 1987, "Interstellar Processes"
	 * table of log nu, log nu*fnu taken from his figure 2 */
	/* >>chng 99 jun 14 energy range lowered to 1e-8 ryd */
	tnuism[0] = 6.00;
	/*tnuism[0] = 9.00;*/
	tnuism[1] = 10.72;
	tnuism[2] = 11.00;
	tnuism[3] = 11.23;
	tnuism[4] = 11.47;
	tnuism[5] = 11.55;
	tnuism[6] = 11.85;
	tnuism[7] = 12.26;
	tnuism[8] = 12.54;
	tnuism[9] = 12.71;
	tnuism[10] = 13.10;
	tnuism[11] = 13.64;
	tnuism[12] = 14.14;
	tnuism[13] = 14.38;
	tnuism[14] = 14.63;
	tnuism[15] = 14.93;
	tnuism[16] = 15.08;
	tnuism[17] = 15.36;
	tnuism[18] = 15.43;
	tnuism[19] = 16.25;
	tnuism[20] = 17.09;
	tnuism[21] = 18.00;
	tnuism[22] = 23.00;
	/* these are log nu*Fnu */
	fnuism[0] = -16.708;
	/*fnuism[0] = -7.97;*/
	fnuism[1] = -2.96;
	fnuism[2] = -2.47;
	fnuism[3] = -2.09;
	fnuism[4] = -2.11;
	fnuism[5] = -2.34;
	fnuism[6] = -3.66;
	fnuism[7] = -2.72;
	fnuism[8] = -2.45;
	fnuism[9] = -2.57;
	fnuism[10] = -3.85;
	fnuism[11] = -3.34;
	fnuism[12] = -2.30;
	fnuism[13] = -1.79;
	fnuism[14] = -1.79;
	fnuism[15] = -2.34;
	fnuism[16] = -2.72;
	fnuism[17] = -2.55;
	fnuism[18] = -2.62;
	fnuism[19] = -5.68;
	fnuism[20] = -6.45;
	fnuism[21] = -6.30;
	fnuism[22] = -11.3;

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

/*lines_table envoked by table lines command, check if we can find all lines in a given list */
int lines_table(void)
{
	long int n,
		miss;

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

	if( !nLINE_TABLE )
		return 0;

	fprintf( ioQQQ , "lines_table checking lines within data table %s\n", chLINE_LIST );
		miss = 0;

	for( n=0; n<nLINE_TABLE; ++n )
	{
		double relative , absolute;
		if( (cdLine( chLabel[n], wl[n] , &relative , &absolute ))<=0 )
		{
			++miss;
			fprintf(ioQQQ,"lines_table in parse_table.c did not find line %4s ",chLabel[n]);
			cdPrtWL(ioQQQ,wl[n]);
			fprintf(ioQQQ,"\n");
		}
	}
	if( miss )
	{
		/* this is so that we pick up problem automatically */
		fprintf(  ioQQQ , "  BOTCHED ASSERTS!!!   Botched Asserts!!! lines_table could not find a total of %li lines\n\n", miss );
	}
	else
	{
		fprintf(  ioQQQ , "lines_table found all lines\n\n" );
	}

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

}
