/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ParseInterp parse parameters on interpolate command */
#include "cddefines.h"
#include "called.h"
#include "rfield.h"
#include "trace.h"
#include "input.h"
#include "parser.h"

void ParseInterp(Parser &p)
{
	DEBUG_ENTRY( "ParseInterp()" );

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

	if( rfield.nShape >= LIMSPC )
	{
		fprintf( ioQQQ, " Too many spectra entered.  Increase LIMSPC\n" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 06 jul 16, fix memory leak when optimizing, PvH */
	if( !rfield.lgContMalloc[rfield.nShape] )
	{
		rfield.tNu[rfield.nShape].resize( NCELL );
		rfield.tslop[rfield.nShape].resize( NCELL );
		rfield.tFluxLog[rfield.nShape].resize( NCELL );
		rfield.lgContMalloc[rfield.nShape] = true;
	}

	strcpy( rfield.chSpType[rfield.nShape], "INTER" );

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

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

	/* this flag says we hit end of command stream */
	p.m_lgEOF = false;
	while( !lgDONE && !p.m_lgEOF )
	{
		/* keep scanning numbers until we hit eol for current line image */
		/* >>chng 01 aug 10, add check on npairs not exceeding NC_ELL as per 
		 * Ilse van Bemmel bug report */
		while( !p.lgEOL() && npairs < NCELL )
		{
			rfield.tNu[rfield.nShape][npairs].set( p.FFmtRead() );
			rfield.tFluxLog[rfield.nShape][npairs] = (realnum)p.FFmtRead();
			++npairs;
		}
		/* check if we fell through due to too many points, did not hit EOL */
		if( !p.lgEOL() )
		{
			fprintf( ioQQQ, "Too many continuum points were entered.\n" );
			fprintf( ioQQQ, 
				"The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL );
			fprintf( ioQQQ, 
				"Increase the value of NCELL in rfield.h.\nSorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}

		/* added too many above, so back off */
		--npairs;

		do
		{
			/* read a new line, checking for EOF */
			p.getline();
			
			/* option to ignore all lines starting with #, *, or %. 
			 * >>chng 06 sep 04 use routine to check for comments */
		}
		while( !p.m_lgEOF && p.isComment());

		/* print the line, but only if it is a continue line */
		if( called.lgTalk && (p.strcmp("CONT")==0) )
		{
			p.echo();
		}

		/* is this a continue line? */
		if( p.strcmp("CONT") != 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( p.last() )
			p.m_lgEOF = true;
	}

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

	/* done reading all of the possible lines */
	if( npairs < 2 )
	{
		fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
		cdEXIT(EXIT_FAILURE);
	}

	if( rfield.tNu[rfield.nShape][0].Ryd() == 0. )
	{
		/* special case - if first energy is zero then it is low energy */
		if( rfield.tNu[rfield.nShape][1].Ryd() > 0. )
		{
			/* second energy positive, assume linear Ryd */
			rfield.tNu[rfield.nShape][0].set(rfield.emm);
		}
		else if( rfield.tNu[rfield.nShape][1].Ryd() < 0. )
		{
			/* second energy negative, assume log of Ryd */
			rfield.tNu[rfield.nShape][0].set(log10(rfield.emm));
		}
		else
		{
			/* second energy zero, not allowed */
			fprintf( ioQQQ, 
				"An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n", 
			  rfield.nShape );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* convert from log(Hz) to Ryd if first nu>5 */
	if( rfield.tNu[rfield.nShape][0].Ryd() >= 5. )
	{
		for( long i=0; i < npairs; i++ )
		{
			rfield.tNu[rfield.nShape][i].set( 
				pow(10.,rfield.tNu[rfield.nShape][i].Ryd())/FR1RYD);
		}
	}
	else if( rfield.tNu[rfield.nShape][0].Ryd() < 0. )
	{
		/* energies entered as logs */
		for( long i=0; i < npairs; i++ )
		{
			rfield.tNu[rfield.nShape][i].set(
				pow(10.,(double)rfield.tNu[rfield.nShape][i].Ryd()));
		}
	}
	else
	{
		/* numbers are linear Rydbergs */
		for( long i=0; i < npairs; i++ )
		{
			if( rfield.tNu[rfield.nShape][i].Ryd() == 0. )
			{
				fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n", 
				  i );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	/* option to print debugging file then exit */
	{
		/* following should be set true to print information */
		enum {DEBUG_LOC=false};
		if( DEBUG_LOC )
		{
			for( long i=0; i < npairs; i++ )
			{
				fprintf(ioQQQ,"%.4e\t%.3e\n",
					rfield.tNu[rfield.nShape][i].Ryd(),
					pow(10.,(double)rfield.tFluxLog[rfield.nShape][i]) * rfield.tNu[rfield.nShape][i].Ryd());
			}
			cdEXIT(EXIT_SUCCESS);
		}
	}

	for( long i=0; i < npairs-1; i++ )
	{
		/* check that frequencies are monotonically increasing */
		if( rfield.tNu[rfield.nShape][i+1].Ryd() <= rfield.tNu[rfield.nShape][i].Ryd() )
		{
			fprintf( ioQQQ, "The energies MUST be in increasing order.  Energy #%3ld=%10.2e Ryd was greater than or equal to the next one.\nSorry.\n", 
				i, rfield.tNu[rfield.nShape][i].Ryd() );
			cdEXIT(EXIT_FAILURE);
		}

		/* TFAC is energy, and TSLOP is slope in f_nu; not photons */
		rfield.tslop[rfield.nShape][i] = (realnum)((rfield.tFluxLog[rfield.nShape][i+1] - 
			rfield.tFluxLog[rfield.nShape][i])/log10(rfield.tNu[rfield.nShape][i+1].Ryd()/
			rfield.tNu[rfield.nShape][i].Ryd()));
	}

	/* now check that array is defined over all energies */
	if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm )
	{
		/* not defined over low energy part of array */
		fprintf( ioQQQ, 
			"\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
		fprintf( ioQQQ, 
			" NOTE You may be making a BIG mistake.\n\n" );
	}

	/* check on IR */
	if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm )
		rfield.lgMMok = false;

	if( rfield.tNu[rfield.nShape][0].Ryd() > 1/36. )
		rfield.lgHPhtOK = false;

	/* gamma ray, EGAMRY is roughly 100MeV */
	if( rfield.tNu[rfield.nShape][npairs-1].Ryd() < rfield.egamry )
		rfield.lgGamrOK = false;

	/* EnerGammaRay is roughly 100keV; high is gamma ray */
	if( rfield.tNu[rfield.nShape][npairs-1].Ryd() < rfield.EnerGammaRay )
		rfield.lgXRayOK = false;

	/* 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, npairs is one less than number of continuum pairs */
	/* If no flux is defined at 1 Ryd, use the nearest endpoint of the supplied spectrum */
	long n;
	for( n=0; n < npairs; n++ )
	{
		if( rfield.tNu[rfield.nShape][n].Ryd() > 1. )
			break;
	}
	/* if present, n is now the table point where rfield.tNuRyd[n-1] is <= 1 ryd,
	 * and rfield.tNuRyd[n] > 1 ryd; max(n-1,0) is the nearest endpoint otherwise */
	realnum fac = rfield.tFluxLog[rfield.nShape][max(n-1,0)];

	for( long i=0; i < npairs; i++ )
	{
		rfield.tFluxLog[rfield.nShape][i] -= fac;
		/*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nShape][i] );*/
	}

	/* option to print out results at this stage - "trace continuum" */
	if( trace.lgConBug && trace.lgTrace )
	{
		fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n",
			 npairs );
		for( long i=0; i < npairs-1; i++ )
		{
			fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n", 
				 i , rfield.tNu[rfield.nShape][i].Ryd(), 
				 rfield.tFluxLog[rfield.nShape][i], rfield.tslop[rfield.nShape][i] );
		}
		fprintf( ioQQQ, "%li\t%.4e\t%.4e\n", 
			 npairs , rfield.tNu[rfield.nShape][npairs].Ryd(), 
			 rfield.tFluxLog[rfield.nShape][npairs]);
	}

	/* finally check that we are within dynamic range of this cpu */
	double cmin = log10( FLT_MIN );
	double cmax = log10( FLT_MAX );
	bool lgHit = false;
	for( long i=0; i < npairs; i++ )
	{
		if( rfield.tFluxLog[rfield.nShape][i] <= cmin || 
		    rfield.tFluxLog[rfield.nShape][i] >= cmax )
		{
			fprintf(ioQQQ,
				" The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
			fprintf(ioQQQ,
				" The frequency is %f and the log of the flux is %f.\n\n",
				rfield.tNu[rfield.nShape][i].Ryd() , 
				rfield.tFluxLog[rfield.nShape][i]);
			lgHit = true;
		}
	}
	if( lgHit )
	{
		fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n");
		fprintf(ioQQQ," NOTE I will try to renormalize it to be within the range of this cpu, but if I crash, this is a likely reason.\n");
		fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n");
		fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n");
		fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n");
	}

	rfield.ncont[rfield.nShape] = npairs;

	/*>>chng 06 may 17, must insure that rNuRyd is sane through NCELL values, not just
	 * nupper values.  This bit when stellar continuum had more cells than cloudy grid,
	 * so npairs is much greater than nupper */
	/* zero out remainder of array */
	rfield.tslop[rfield.nShape][npairs-1] = 0.;
	for( long i=npairs; i < NCELL; i++ )
	{
		/* >>chng 05 jul 27, next three indices had been ipspec, changed to nShape
		 * bug caught by Ralf Quast */
		rfield.tFluxLog[rfield.nShape][i] = 0.;
		rfield.tslop[rfield.nShape][i] = 0.;
		rfield.tNu[rfield.nShape][i].set(0.);
	}

	/* now increment number of continua */
	++rfield.nShape;

	return;
}