/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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 "parse.h"

void ParseInterp(char *chCard, 
  int *lgEOF)
{
	char chLab4[5];
	int lgDONE, 
	  lgEOL,
	  lgHit;
	long int i, 
		k ,
	  npairs;
	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
	 */

	if( rfield.nspec >= LIMSPC-1 )
	{
		fprintf( ioQQQ, " Too many spectra entered.  Increase LIMSPC\n" );
		puts( "[Stop in ParseInterp]" );
		cdEXIT(EXIT_FAILURE);
	}
	rfield.lgContMalloc[rfield.nspec] = TRUE;
	if( (rfield.tNuRyd[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.tslop[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.tFluxLog[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) ) ) == NULL )
		BadMalloc();

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

	/* read all of the numbers on a line */
	npairs = 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 */
		/* >>chng 01 aug 10, add check on npairs not exceeding NC_ELL as per 
		 * Ilse van Bemmel bug report */
		/*while( !lgEOL && npairs<rfield.nupper )*/
		while( !lgEOL && npairs<NCELL )
		{
			rfield.tNuRyd[rfield.nspec][npairs] = 
				(float)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
			rfield.tFluxLog[rfield.nspec][npairs] = 
				(float)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
			++npairs;
		}
		/* check if we fell through due to too many points, did not hit EOL */
		if( !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" );
			puts( "[Stop in ParseInterp]" );
			cdEXIT(EXIT_FAILURE);
		}

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

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

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

		/* get cap'd first four char of chCard */
		cap4( chLab4 , chCard );

		/* print the line, but only if it is a continue line */
		if( called.lgTalk && (strncmp(chLab4,"CONT",4)==0) )
		{
			fprintf( ioQQQ, "                       * ");
			/* next logic will print command plus spaces to right justified * */
			k=0;
			while( chCard[k]!='\0' )
			{
				fprintf(ioQQQ,"%c",chCard[k]);
				++k;
			}
			while( k<80 )
			{
				fprintf(ioQQQ,"%c",' ');
				++k;
			}
			fprintf( ioQQQ,"*\n"); 
		}

		/* now convert entire line 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 */
	--npairs;
	if( npairs < 1 )
	{
		fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
		puts( "[Stop in ParseInterp]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 02 apr 19, from > to >=, test did not trigger due to overrun protection
	 * in main loop */
	else if( npairs >= NCELL - 2 )
	{
		fprintf( ioQQQ, " Too many continuum points entered.\n" );
		fprintf( ioQQQ, 
			"The current logic limits the number of possible points to the value of NCELL, which is %i.\nSorry.\n",NCELL );
		fprintf( ioQQQ, " Increase NCELL (in cddefines.h) to more than the number of continuum points.\n" );
		puts( "[Stop in ParseInterp]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( rfield.tNuRyd[rfield.nspec][0] == 0. )
	{
		/* special case - if first energy is zero then it is low energy */
		if( rfield.tNuRyd[rfield.nspec][1] > 0. )
		{
			/* second energy positive, assume linear Ryd */
			rfield.tNuRyd[rfield.nspec][0] = (float)rfield.emm;
		}
		else if( rfield.tNuRyd[rfield.nspec][1] < 0. )
		{
			/* second energy negative, assume log of Ryd */
			rfield.tNuRyd[rfield.nspec][0] = (float)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.nspec );
			puts( "[Stop in ParseInterp]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

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

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

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

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

	for( i=0; i < npairs; i++ )
	{
		/* check that frequencies are monotonically increasing */
		if( rfield.tNuRyd[rfield.nspec][i+1] <= rfield.tNuRyd[rfield.nspec][i] )
		{
			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.tNuRyd[rfield.nspec][i] );
			puts( "[Stop in ParseInterp]" );
			cdEXIT(EXIT_FAILURE);
		}

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

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

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

	/* check on IR */
	if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
	{
		rfield.lgMMok = FALSE;
	}

	if( rfield.tNuRyd[rfield.nspec][0] > 1/36. )
	{
		rfield.lgHPhtOK = FALSE;
	}

	/* gamma ray, EGAMRY is roughly 100MeV */
	if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.egamry )
	{
		rfield.lgGamrOK = FALSE;
	}

	/* EnerGammaRay is roughly 100keV; high is gamma ray */
	if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.EnerGammaRay )
	{
		rfield.lgXRayOK = FALSE;
	}

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

	/* tfac can be very large or small */
	/* >>chng 01 jul 11, from <npairs to <=npairs, [nparis] is highest point in table */
	for( i=0; i <= npairs; i++ )
	{
		cmax = MAX2(cmax,rfield.tFluxLog[rfield.nspec][i]);
		cmin = MIN2(cmin,rfield.tFluxLog[rfield.nspec][i]);
	}

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

	if( cmin < -37. )
	{
		fac = -cmin - 37.;
		/* >>chng 01 jul 11, from <npairs to <=npairs, [nparis] is highest point in table */
		for( i=0; i <= npairs; i++ )
		{
			rfield.tFluxLog[rfield.nspec][i] += (float)fac;
		}
	}

	else if( cmax > 37. )
	{
		fac = cmax - 37.;
		/* >>chng 01 jul 11, from <npairs to <=npairs, [nparis] is highest point in table */
		for( i=0; i <= npairs; i++ )
		{
			rfield.tFluxLog[rfield.nspec][i] -= (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 < (npairs + 1); i++ )
		{
			fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNuRyd[rfield.nspec][i], 
			  rfield.tFluxLog[rfield.nspec][i], rfield.tslop[rfield.nspec][i] );
		}
	}

	/* 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 */
	i = 0;
	while( rfield.tNuRyd[rfield.nspec][i] <= 1. && i < npairs-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]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* do the renormalization so 1 ryd is about unity */
	fac = rfield.tFluxLog[rfield.nspec][i-1];
	for( i=0; i <= npairs; i++ )
	{
		rfield.tFluxLog[rfield.nspec][i] -= (float)fac;
		/*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nspec][i] );*/
	}

	/* finally check that we are within dynamic range of this cpu */
	cmin = log10( FLT_MIN );
	cmax = log10( FLT_MAX );
	lgHit = FALSE;
	for( i=0; i <= npairs; i++ )
	{
		if( rfield.tFluxLog[rfield.nspec][i] <= cmin || 
			rfield.tFluxLog[rfield.nspec][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.tNuRyd[rfield.nspec][i] , 
				rfield.tFluxLog[rfield.nspec][i]);
			lgHit = TRUE;
		}
	}
	if( lgHit )
	{
		fprintf(ioQQQ,"\n The log of the flux given in an interpolate command is outside the range of this cpu.\n");
		fprintf(ioQQQ," I will try to renormalize it to be within the range of this cpu, but I crash, this is a likely reason.\n");
		fprintf(ioQQQ," Note that the interpolate command only is used for the shape of the continuum.\n");
		fprintf(ioQQQ," The order of magnitude of the flux is not used in any way.\n");
		fprintf(ioQQQ," For safety this could be of order unity.\n\n");
	}

	/* zero out remainder of array */
	for( i=npairs; i < rfield.nupper; i++ )
	{
		/* >>chng 05 jul 27, next three indices had been ipspec, changed to nspec
		 * bug caught by Ralf Quast */
		rfield.tFluxLog[rfield.nspec][i] = 0.;
		rfield.tslop[rfield.nspec][i] = 0.;
		rfield.tNuRyd[rfield.nspec][i] = 0.;;
	}

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

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