/*ParseInterp parse parameters on interpolate command */
#include <string.h>
#include "cddefines.h"
#include "incidspec.h"
#include "called.h"
#include "bndsok.h"
#include "bounds.h"
#include "egray.h"
#include "trace.h"
#include "bit32.h"
#include "input.h"
#include "ffmtread.h"
#include "readar.h"
#include "cap4.h"
#include "showme.h"
#include "caps.h"
#include "parse.h"

void ParseInterp(char *chCard, 
  int *lgEOF)
{
	char chLab4[5];
	int lgDONE, 
	  lgEOL;
	long int i, 
		k ,
	  j;
	double cmax, 
	  cmin, 
	  fac;

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

	/* 
	 * this sub reads in the "interpolate" command, and has
	 * logic for the "continue" lines as well.
	 * OUTPUT:
	 * TNU is vector of eneriges where the grid is defined
	 * TSLOP initially is vector of log fnu at each freq
	 * converted into slopes here too
	 */

	IncidSpec.nspec += 1;
	if( IncidSpec.nspec >= LIMSPC )
	{
		fprintf( ioQQQ, " Too many spectra entered.  Increase LIMSPC\n" );
		puts( "[Stop in ParseInterp]" );
		exit(1);
	}

	strcpy( spnorm.chSpType[IncidSpec.nspec-1], "INTER" );

	/* read all of the numbers on a line */
	j = 0;

	/* this is flag saying that all numbers are in */
	lgDONE = FALSE;

	/* this flag says we hit end of command stream */
	*lgEOF = FALSE;
	while( !lgDONE && !*lgEOF )
	{
		i = 5;
		lgEOL = FALSE;

		/* keep scanning numbers til we hit eol for current line image */
		while( !lgEOL )
		{
			j = MIN2(j+1,NCELL);
			IncidSpec.tnu[j-1][IncidSpec.nspec-1] = (float)FFmtRead(chCard
			  ,&i,76,&lgEOL);
			IncidSpec.tslop[j-1][IncidSpec.nspec-1] = (float)FFmtRead(chCard
			  ,&i,76,&lgEOL);
		}
		--j;

		/* read a new line, checking for EOF */
		readar(chCard,lgEOF);

		/* option to ignore all lines starting with #, *, or %. */
		while( !*lgEOF && ((chCard[0] == '#' || chCard[0] == '*') || chCard[0] == '%') )
		{
			readar(chCard,lgEOF);
		}

		/* get cap'd first four char of chCard */
		cap4( chLab4 , chCard );
		if( called.lgTalk && (strncmp(chLab4,"CONT",4)==0) )
		{
			fprintf( ioQQQ, "                       * ");
			k=0;
			while( chCard[k]!='\0' )
			{
				fprintf(ioQQQ,"%c",chCard[k]);
				++k;
			}
			while( k<80 )
			{
				fprintf(ioQQQ,"%c",' ');
				++k;
			}
			fprintf( ioQQQ,"*\n"); 
		}

		/* now convert to caps */
		caps(chCard);

		/* is this a continue line? */
		if( strncmp(chCard,"CONT",4) != 0 )
		{
			/* we have a line but next command, not continue */
			lgDONE = TRUE;
		}

		/* this is another way to hit end of input stream - blank lines */
		if( chCard[0] == ' ' )
			*lgEOF = TRUE;
	}

	/* if valid next line, backup one line */
	if( lgDONE )
	{
		input.nRead -= input.iReadWay;
	}

	/* done reading all of the possible lines */
	j -= 1;
	if( j < 1 )
	{
		fprintf( ioQQQ, " There must be at least 2 pairs to interpolate,\n" );
		puts( "[Stop in ParseInterp]" );
		exit(1);
	}

	else if( j > NCELL - 2 )
	{
		fprintf( ioQQQ, " Too many continuum points entered.\n" );
		puts( "[Stop in ParseInterp]" );
		exit(1);
	}

	if( IncidSpec.tnu[0][IncidSpec.nspec-1] == 0. )
	{
		/* special case - if first energy is zero then it is low energy */
		if( IncidSpec.tnu[1][IncidSpec.nspec-1] > 0. )
		{
			/* second energy positive, assume linear Ryd */
			IncidSpec.tnu[0][IncidSpec.nspec-1] = bounds.emm;
		}
		else if( IncidSpec.tnu[1][IncidSpec.nspec-1] < 0. )
		{
			/* second energy negative, assume log of Ryd */
			IncidSpec.tnu[0][IncidSpec.nspec-1] = (float)log10(bounds.emm);
		}
		else
		{
			/* second energy zero, not allowed */
			fprintf( ioQQQ, 
				" An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.  SORRY\n", 
			  IncidSpec.nspec-1 );
			puts( "[Stop in ParseInterp]" );
			exit(1);
		}
	}

	/* 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.517);
		}
	}

	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]);
		}
	}

	else
	{
		/* numbers are linear Rydbergs */
		for( i=0; i < (j + 1); i++ )
		{
			if( IncidSpec.tnu[i][IncidSpec.nspec-1] == 0. )
			{
				fprintf( ioQQQ, " An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.  SORRY\n", 
				  i );
				puts( "[Stop in ParseInterp]" );
				exit(1);
			}
		}
	}

	for( i=0; i < j; i++ )
	{
		/* check that frequencies are monotonically increasing */
		if( IncidSpec.tnu[i+1][IncidSpec.nspec-1] <= IncidSpec.tnu[i][IncidSpec.nspec-1] )
		{
			fprintf( ioQQQ, " The energies MUST be in increasing order.  Energy #%3ld=%10.2e Ryd was greater than or equal to the next one.  Sorry.\n", 
			  i, IncidSpec.tnu[i][IncidSpec.nspec-1] );
			puts( "[Stop in ParseInterp]" );
			exit(1);
		}

		/* TFAC is energy, and TSLOP is slope in f_nu; not photons */
		IncidSpec.tfac[i][IncidSpec.nspec-1] = IncidSpec.tslop[i][IncidSpec.nspec-1];
		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]));
	}
	IncidSpec.tfac[j][IncidSpec.nspec-1] = IncidSpec.tslop[j][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.;
		}
	}

	/* now check that array is defined over all energies */
	if( IncidSpec.tnu[0][IncidSpec.nspec-1] > bounds.emm )
	{
		/* not defined over low energy part of array */
		fprintf( ioQQQ, " Warning:  the input continuum was not defined over the entire energy range. Some energies are set to zero.\n >>>>>>>>>>You are making a big mistake.\n" );
	}

	/* check on IR */
	if( IncidSpec.tnu[0][IncidSpec.nspec-1] > bounds.emm )
	{
		bndsok.lgMMok = FALSE;
	}

	if( IncidSpec.tnu[0][IncidSpec.nspec-1] > 1/36. )
	{
		bndsok.lgHPhtOK = FALSE;
	}

	/* gamma ray, EGAMRY is roughly 100MeV */
	if( IncidSpec.tnu[j][IncidSpec.nspec-1] < bounds.egamry )
	{
		bndsok.lgGamrOK = FALSE;
	}

	/* EnerGammaRay is roughly 100keV; high is gamma ray */
	if( IncidSpec.tnu[j][IncidSpec.nspec-1] < egray.EnerGammaRay )
	{
		bndsok.lgXRayOK = FALSE;
	}

	/* find min and max of continuum */
	cmax = -38.;
	cmin = 28;

	/* tfac can be very large or small */
	for( i=0; i < j; i++ )
	{
		cmax = MAX2(cmax,IncidSpec.tfac[i][IncidSpec.nspec-1]);
		cmin = MIN2(cmin,IncidSpec.tfac[i][IncidSpec.nspec-1]);
	}

	/* check on dynamic range of input continuum */
	if( cmax - cmin > 74. && bit32.lgBit32 )
	{
		fprintf( ioQQQ, " The dynamic range of the specified continuumis too large for a 32-bit computer.  Sorry.\n" );
		puts( "[Stop in ParseInterp]" );
		exit(1);
	}

	if( cmin < -37. )
	{
		fac = -cmin - 37.;
		for( i=0; i < j; i++ )
		{
			IncidSpec.tfac[i][IncidSpec.nspec-1] += (float)fac;
		}
	}

	else if( cmax > 37. )
	{
		fac = cmax - 37.;
		for( i=0; i < j; i++ )
		{
			IncidSpec.tfac[i][IncidSpec.nspec-1] -= (float)fac;
		}
	}

	/* option to print out results at this stage - "trace continuum" */
	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 flux we will interpolate upon passes through unity
	 * at near 1 Ryd.  but first we must find 1 Ryd in the array.
	 * find 1 ryd, j is one less than number of continuum pairs */
	i = 0;
	while( IncidSpec.tnu[i][IncidSpec.nspec-1] <= 1. && i < j-1 )
	{
		i += 1;
	}
	/* i is now the table point where tnu[i-1] is <= 1 ryd,
	 * and tnu[i] > 1 ryd */

	/* following is impossible but sanity check */
	if( i < 1 )
	{
		fprintf( ioQQQ, " ParseInput finds insane i after tnu loop\n");
		ShowMe();
		puts( "[Stop in ParseInterp]" );
		exit(1);
	}

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


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