/*ParseTable parse the table read command */
#include <string.h>
#include <float.h>
#include <limits.h>
#include "cddefines.h"
#include "physconst.h"
#include "incidspec.h"
#include "varypar.h"
#include "bounds.h"
#include "trace.h"
#include "radius.h"
#include "input.h"
#include "ffmtread.h"
#include "lgmatch.h"
#include "getquote.h"
#include "atlas.h"
#include "werner.h"
#include "kurucz79.h"
#include "mihalas.h"
#include "rauch.h"
#include "costar.h"
#include "parse.h"

/* 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	45
double tnurbn[NRUBIN], 
  fnurbn[NRUBIN];

/* 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];

/* Mathews and Ferland generic agn continuum */
#define	NAGN	8
static double tnuagn[NAGN], 
  tslagn[NAGN];
/* this is below and stores values for above vectors */
void ZeroContin(void);

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

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

	long int i, 
	  iabund,
	  j, 
	  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;
	}

	IncidSpec.nspec += 1;

	if( IncidSpec.nspec > LIMSPC )
	{
		fprintf( ioQQQ, " %4ld is too many spectra entered.  Increase LIMSPC\n", 
		  IncidSpec.nspec );
		puts( "[Stop in ParseTable]" );
		exit(1);
	}

	/* set flag telling interpolate */
	strcpy( spnorm.chSpType[IncidSpec.nspec-1], "INTER" );

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

		j = NAGN - 1;

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

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

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

			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" );
				puts( "[Stop in ParseTable]" );
				exit(1);
			}

			if( ConBreak >= IncidSpec.tnu[2][IncidSpec.nspec-1] )
			{
				fprintf( ioQQQ, " The energy of the break cannot be greater than%10.2e Ryd.\n", 
				  IncidSpec.tnu[2][IncidSpec.nspec-1] );
				puts( "[Stop in ParseTable]" );
				exit(1);
			}

			else if( ConBreak <= IncidSpec.tnu[0][IncidSpec.nspec-1] )
			{
				fprintf( ioQQQ, " The energy of the break cannot be less than%10.2e Ryd.\n", 
				  IncidSpec.tnu[0][IncidSpec.nspec-1] );
				puts( "[Stop in ParseTable]" );
				exit(1);
			}

			IncidSpec.tnu[1][IncidSpec.nspec-1] = (float)ConBreak;

			IncidSpec.tslop[1][IncidSpec.nspec-1] = 
				(float)(IncidSpec.tslop[2][IncidSpec.nspec-1] + 
			  log10(IncidSpec.tnu[2][IncidSpec.nspec-1]/IncidSpec.tnu[1][IncidSpec.nspec-1]));

			IncidSpec.tslop[0][IncidSpec.nspec-1] = 
				(float)(IncidSpec.tslop[1][IncidSpec.nspec-1] - 
			  2.5*log10(IncidSpec.tnu[1][IncidSpec.nspec-1]/IncidSpec.tnu[0][IncidSpec.nspec-1]));
		}
	}

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

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

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

	else if( lgMatch(" ISM",chCard) )
	{
		/* local ISM radiation field from Black 1987, Interstellar Processes */
		for( i=0; i < NISM; i++ )
		{
			IncidSpec.tnu[i][IncidSpec.nspec-1] = (float)tnuism[i];
			IncidSpec.tslop[i][IncidSpec.nspec-1] = (float)(fnuism[i] - 
			  tnuism[i]);
		}

		j = NISM - 1;
		i = 5;

		/* optional scale factor to change default luminosity
		 * from observed value
		 * assumed to be log if negative, and linear otherwise */
		scale = FFmtRead(chCard,&i,76,&lgEOL);
		if( scale > 0. )
			scale = log10(scale);

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

		/* 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( IncidSpec.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]" );
			exit(1);
		}

		strcpy( spnorm.chRSpec[*nqh-1], "SQCM" );
		strcpy( spnorm.chSpNorm[*nqh-1], "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 */
		IncidSpec.range[0][*nqh-1] = (float)HIONPOT;

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

		/* set radius to very large value if not already set */
		if( radius.Radius == 0. )
		{
			*ar1 = (float)radius.rdfalt;
			radius.Radius = pow(10.,*ar1);
		}

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

		/* this is low energy for code */
		IncidSpec.tnu[0][IncidSpec.nspec-1] = bounds.emm;
		IncidSpec.tslop[0][IncidSpec.nspec-1] = -5;

		lgLogSet = FALSE;

		/* option to adjust sub-millimeter break */
		brakmm = FFmtRead(chCard,&i,76,&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 = IncidSpec.tnu[0][IncidSpec.nspec-1]*1.001;
		}

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

		IncidSpec.tnu[1][IncidSpec.nspec-1] = (float)brakmm;

		/* this is energy slope between low energy limit and this break */
		slopir = 5./2.;
		IncidSpec.tslop[1][IncidSpec.nspec-1] = 
			(float)(IncidSpec.tslop[0][IncidSpec.nspec-1] + 
		  slopir*log10(IncidSpec.tnu[1][IncidSpec.nspec-1]/IncidSpec.tnu[0][IncidSpec.nspec-1]));

		/* this is high energy limit to code */
		IncidSpec.tnu[3][IncidSpec.nspec-1] = bounds.egamry;

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

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

		else if( brakxr == 0. )
		{
			brakxr = IncidSpec.tnu[3][IncidSpec.nspec-1]/1.001;
		}

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

		IncidSpec.tnu[2][IncidSpec.nspec-1] = (float)brakxr;

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

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

		IncidSpec.tslop[3][IncidSpec.nspec-1] = 
			(float)(IncidSpec.tslop[2][IncidSpec.nspec-1] + 
		  slopxr*log10(IncidSpec.tnu[3][IncidSpec.nspec-1]/IncidSpec.tnu[2][IncidSpec.nspec-1]));

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

	else if( lgMatch("READ",chCard) )
	{
		/* 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 );

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

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

		/* this will be flag saying continuum not read in here */
		IncidSpec.tslop[0][IncidSpec.nspec-1] = 0.;
		IncidSpec.tslop[1][IncidSpec.nspec-1] = 0.;
		/* nothing left to do here, better return
		 * this file is actually read in by routine ReadTable */
		
#		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++ )
		{
			IncidSpec.tnu[i][IncidSpec.nspec-1] = (float)tnurbn[i];
			IncidSpec.tslop[i][IncidSpec.nspec-1] = (float)log10(fnurbn[i]);
		}
		j = NRUBIN - 1;
	}

	else if( lgMatch("STAR",chCard) )
	{
		/* TABLE STAR, numbers are temperature and log(g) */
		i = 5;
		par1 = FFmtRead(chCard,&i,76,&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]" );
			exit(1);
		}

		/* lgEOL is used in following logic to check whether grav set */
		par2 = FFmtRead(chCard,&i,76,&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;
			}

			DoAtlas(&nstar,par1,par2,lgEOL);
			/* set flag saying really just read in continuum exactly as punched */
			strcpy( spnorm.chSpType[IncidSpec.nspec-1], "VOLK " );
			j = 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;

			DoWerner(&nstar,par1,par2);

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

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

		else if( lgMatch("COST",chCard) )
		{
			int nmodid = INT_MAX;
			/* >>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]" );
					exit(1);
				}
			}
			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]" );
					exit(1);
				}
			}
			else
			{
				/* this is Teff, log(g) mode */
				/* Caution: the code in DoCoStar 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;
				}
			}
			DoCoStar(&nstar,intmode,par1,par2,nmodid);

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

			j = 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 */
         DoRauch(&nstar,par1,par2,iabund);

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

			j = 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);
			j = nstar - 1;
		}

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

			mihals(&nstar,par1,par2);
			j = 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]" );
			exit(1);
		}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	/* 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(spnorm.chSpType[IncidSpec.nspec-1],"VOLK ") == 0 )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->ParseTable()\n", debug_fp );
#		endif
		return;
	}

	/* convert from log(Hz) to Ryd if first nu>5 */
	if( IncidSpec.tnu[0][IncidSpec.nspec-1] >= 5. )
	{
		for( i=0; i < (j + 1); i++ )
		{
			IncidSpec.tnu[i][IncidSpec.nspec-1] = 
				(float)pow(10.,IncidSpec.tnu[i][IncidSpec.nspec-1] - 15.51718);
		}
	}

	else if( IncidSpec.tnu[0][IncidSpec.nspec-1] < 0. )
	{
		/* energies entered as logs */
		for( i=0; i < (j + 1); i++ )
		{
			IncidSpec.tnu[i][IncidSpec.nspec-1] = 
				(float)pow(10.,IncidSpec.tnu[i][IncidSpec.nspec-1]);
		}
	}

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

	for( i=0; i < j; i++ )
	{
		IncidSpec.tslop[i][IncidSpec.nspec-1] = 
			(float)((IncidSpec.tslop[i+1][IncidSpec.nspec-1] - 
		  IncidSpec.tslop[i][IncidSpec.nspec-1])/
		  log10(IncidSpec.tnu[i+1][IncidSpec.nspec-1]/
		  IncidSpec.tnu[i][IncidSpec.nspec-1]));
	}

	if( j + 2 < NCELL )
	{
		/* zero out remainder of array */
		for( i=j + 1; i < NCELL; i++ )
		{
			IncidSpec.tnu[i][IncidSpec.nspec-1] = 0.;
		}
	}

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

	/* 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 = TRUE;
	/*following will be reset when hit occurs*/
	fac = -DBL_MAX;
	for( i=0; i < (j + 1); i++ )
	{
		if( IncidSpec.tnu[i][IncidSpec.nspec-1] > 1. )
		{
			fac = IncidSpec.tfac[i][IncidSpec.nspec-1];
			lgHit = FALSE;
		}
	}

	if( lgHit )
	{
		fprintf( ioQQQ, " TABLE fails sanity test in table\n" );
		puts( "[Stop in ParseTable]" );
		exit(1);
	}

	/* do the renormalization */
	for( i=0; i < (j + 1); i++ )
	{
		IncidSpec.tfac[i][IncidSpec.nspec-1] -= (float)fac;
	}

#	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( bounds.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] = bounds.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] = bounds.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

	/* generic AGN continuum taken from 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 Davidson and Fesen, 1985 Ann Rev
	 * 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 */
	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;
	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;
}