/* 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 */
/*ParseDriveCmnd parse the drive command - drive calls to various subs */
/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
/*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
/*dgaunt drive gaunt factor routines by letting user query values */
#include "cddefines.h"
#include "trace.h"
#include "hydroeinsta.h"
#include "phycon.h"
#include "tfidle.h"
#include "atmdat.h"  
#include "taulines.h"
#include "abund.h"
#include "helike.h" 
#include "rt.h"
#include "continuum.h"
#include "parse.h"
#include "gaunt.h"
#include "physconst.h"

/*dgaunt drive gaunt factor routines by letting user query values */
static void dgaunt(void);

/*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
static void DrvHyas(void);

/* drive escape probability routines */
void DrvEscP( void );

/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
static void DrvCaseBHS(void);

void ParseDriveCmnd(char *chCard )
{
	int lgEOL;
	long int n, 
	  i;
	double fac, 
	  zed;

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

	/* NB evolve all following names to style DrvSomething */

	/* option to drive cloudy, which one? */
	if( lgMatch("FFMT",chCard) )
	{
		/* free format parser */
		char chInput[INPUT_LINE_LENGTH];
		fprintf( ioQQQ, " FFmtRead ParseDriveCmnd entered.  Enter number.\n" );
		lgEOL = FALSE;
		while( !lgEOL )
		{
			if( fgets( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
			{
				fprintf( ioQQQ, " ParseDriveCmnd.dat error getting magic number\n" );
				puts( "[Stop in ParseDriveCmnd]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " FFmtRead hit the EOL with no value, return=%10.2e\n", 
				  fac );
				break ;
			}
			else if( fac == 0. )
			{
				break;
			}
			else
			{
				fprintf( ioQQQ, " FFmtRead returned with value%11.4e\n", 
				  fac );
			}
			fprintf( ioQQQ, " Enter 0 to stop, or another value.\n" );
		}
		fprintf( ioQQQ, " FFmtRead ParseDriveCmnd exits.\n" );
	}

	else if( lgMatch("CASE",chCard) )
	{
		/* option to interpolate on Hummer and Storey case b hydrogenic emission routines */
		DrvCaseBHS( );
	}

	else if( lgMatch("CDLI",chCard) )
	{
		/* drive cdLine to check that it finds all the right lines, routine is in lines.c */
		trace.lgDrv_cdLine = TRUE;
	}

	else if( lgMatch("ESCA",chCard) )
	{
		/* option to drive escape probability routines */
		DrvEscP( );
	}

	else if( lgMatch("HELI",chCard) )
	{
		fprintf(ioQQQ," output will go to output.txt\n");
		/* option to drive Ryan Porter he-like helium routine,
		 * routine in helike.c */
		HeCreate();
	}

	else if( lgMatch("HYAS",chCard) )
	{
		/* option to drive Jason's hydrogen transition probabilities */
		DrvHyas();
	}

	else if( lgMatch("GAUN",chCard) )
	{
		/* drive gaunt factor routine */
		dgaunt();
	}

	else if( lgMatch("POIN",chCard) )
	{
		/* later on, check cell pointers, centers, widths */
		fprintf( ioQQQ, " Define entire model first, later I will ask for energies.\n" );
		trace.lgPtrace = TRUE;
	}

	else if( lgMatch("PUMP",chCard) )
	{
		char chInput[INPUT_LINE_LENGTH];
		lgEOL = FALSE;
		fprintf( ioQQQ, " Continuum pump ParseDriveCmnd entered - Enter log tau\n" );
		while( !lgEOL )
		{
			if( fgets( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
			{
				fprintf( ioQQQ, " ParseDriveCmnd.dat error getting magic number\n" );
				puts( "[Stop in ParseDriveCmnd]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* get tau */
			i = 1;
			fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
				break;
			fac = pow(10.,fac);
			fprintf( ioQQQ, " Tau =%11.4e\n", fac );
			TauDummy.TauIn = (float)fac;
			fac = ContPump(&TauDummy);
			fprintf( ioQQQ, " ContPump =%11.4e\n", fac );
			fprintf( ioQQQ, " Enter null to stop, or another value.\n" );
		}
		fprintf( ioQQQ, " ContPump ParseDriveCmnd exits.\n" );
	}

	else if( lgMatch("STAR",chCard) )
	{
		char chInput[INPUT_LINE_LENGTH];
		/* get starburst abundances */
		for( i=0; i < 40; i++ )
		{
			zed = (float)(i+1)/4. + 0.01;
			sprintf( chInput, "starburst, zed=%10.4f", zed );
			abund_starburst(chInput);
			fprintf( ioQQQ, "%8.1e", zed );
			for(n=0; n < LIMELM; n++)
				fprintf( ioQQQ, "%8.1e", abund.solar[n] );
			fprintf( ioQQQ, "\n" );
		}
	}

	else
	{
		fprintf( ioQQQ, 
			" Unrecognized key; keys are FFmtRead, CASE, GAUNT, HYAS, POINT, MOLE, STAR, "
			"PUMP, and ESCApe.  Sorry.\n" );
		puts( "[Stop in ParseDriveCmnd]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*DrvEscP user queries escape probability routines, which return values */

void DrvEscP( void )
{
	char chCard[INPUT_LINE_LENGTH];
	int lgEOL ;
	long i;
	double tau;

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

	/* this routine is enterd with the command escape probability, and
	 * drives the escape probability routine to compare answers */
	fprintf( ioQQQ, " Enter the log of the one-sided optical depth; line with no number to stop.\n" );

	lgEOL = FALSE;
	while( !lgEOL )
	{
		if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
		{
			break;
		}

		i = 1;
		tau = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			break;
		}

		tau = pow(10.,tau);
		fprintf( ioQQQ, "tau was %e\n", tau );
		fprintf( ioQQQ, " ESCINC=%11.3e\n", esc_PRD_1side(tau,1e-4) );
		fprintf( ioQQQ, " ESCCOM=%11.3e\n", esc_CRDwing_1side(tau,1e-4 ) );
		fprintf( ioQQQ, " ESCA0K2=%11.3e\n", esca0k2(tau) );

	}

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

/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
static void DrvCaseBHS(void)
{
	char chCard[INPUT_LINE_LENGTH];
	int lgEOL,
		lgHit;
	long int i, 
	  n1, 
	  nelem ,
	  n2;
	double Temperature,
		Density;

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

	/* this routine is entered with the command DRIVE CASEB, and
	 * interpolates on the Hummer & Storey case b data set */

	/* read in some external data files, but only if this is first call */
	fprintf(ioQQQ," I will get needed H data files. This will take a second.\n");
	atmdat_readin();

	{
		/* following should be set true to print input lines */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			double xLyman , alpha;
			long int ipHi;
			nelem = 2;
			Temperature = 2e4;
			Density = 1e2;
			for( ipHi=3; ipHi<25; ++ipHi )
			{
				double photons = (1./POW2(ipHi-1.)-1./POW2((double)ipHi) ) /(1.-1./ipHi/ipHi );
				xLyman = atmdat_HS_caseB(1,ipHi, nelem,Temperature , Density , 'A' ) ;
				alpha = atmdat_HS_caseB(ipHi-1,ipHi, nelem,Temperature , Density , 'A' );
				fprintf(ioQQQ,"%li\t%.3e\t%.3e\n", ipHi, xLyman/alpha*photons, photons );
			}
			cdEXIT(EXIT_SUCCESS);
		}
	}

	/* first get the charge, only H and He at present */
	lgHit = FALSE;
	nelem = 0;
	while( !lgHit )
	{
		fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" );
		if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
		{
			fprintf( ioQQQ, " error getting species \n" );
			puts( "[Stop in DrvCaseBHS]" );
			cdEXIT(EXIT_FAILURE);
		}

		i = 1;
		nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL || nelem< 1 || nelem > 2 )
		{
			fprintf( ioQQQ, " must be either 1 or 2!\n" );
		}
		else
		{
			lgHit = TRUE;
		}
	}

	fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n");
	fprintf(ioQQQ," The density is always a log.\n");
	fprintf(ioQQQ," The order of the quantum numbers do not matter.\n");
	fprintf(ioQQQ," The smallest must not be smaller than 2,\n");
	fprintf(ioQQQ," and the largest must not be larger than 25.\n");
	fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n");
	fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n");

	lgHit = TRUE;
	/* this is always true */
	while( lgHit )
	{
		fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" );
		if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
		{
			fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
			puts( "[Stop in DrvCaseBHS]" );
			cdEXIT(EXIT_FAILURE);
		}

		i = 1;
		Temperature = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting temperature!\n" );
			break;
		}

		/* log if less than 10 */
		if( Temperature < 10. )
		{
			Temperature = pow(10., Temperature );
		}
		fprintf(ioQQQ," Temperature is %g\n", Temperature );

		Density = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting density!\n" );
			break;
		}
		Density = pow(10., Density );
		fprintf(ioQQQ," Density is %g\n", Density );

		/* these quantum numbers can be in any order */
		n1 = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting quantum number!\n" );
			break;
		}

		n2 = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting quantum number!\n" );
			break;
		}

		if( MAX2( n1 , n2 ) > 25 )
		{
			fprintf( ioQQQ," The limits of the HS tables are 2 <= n <= 25.  Sorry.\n");
			break;
		}

		fprintf( ioQQQ, 
			" 4pJ(%ld,%ld)/n_e n_p=%11.3e\n", 
			 n1, n2, 
			 atmdat_HS_caseB(n1,n2, nelem,Temperature , Density , 'B' ) );
		/* this is check that we were in bounds */

	}

	fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
	puts( "[Stop in DrvCaseBHS]" );
	cdEXIT(EXIT_FAILURE);

}

/*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
static void DrvHyas(void)
{
	char chCard[INPUT_LINE_LENGTH];
	int lgEOL;
	long int i, 
	  ihi, 
	  low;

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

	/* this routine is entered with the command DRIVE HYAS, and
	 * drives Jason's hydrogen einstein A routines */

	ihi = 1;
	/* ihi never lt 1 */
	while( ihi != 0 )
	{
		fprintf( ioQQQ, " Enter two quantum numbers, null line to stop.\n" );
		if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
		{
			fprintf( ioQQQ, " error getting drvhyas \n" );
			puts( "[Stop in drvhyas]" );
			cdEXIT(EXIT_FAILURE);
		}

		i = 1;
		low = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
			break;

		ihi = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " must be two numbers!\n" );
			break;
		}
		fprintf( ioQQQ, " A(%3ld,%3ld)=%11.3e\n", 
			MIN2(low,ihi), MAX2(low,ihi), HydroEinstA(low,ihi) );

	}
	fprintf( ioQQQ, " Driver exits, enter next line.\n" );


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

/*dgaunt drive gaunt factor routines by letting user query values */
static void dgaunt(void)
{
	char chCard[INPUT_LINE_LENGTH];
	int lgEOL, inputflag;
	long int i, 
	  ierror;
	float enerlin[1],
	  SaveTemp;
	double z,mygaunt=0.;
	double loggamma2, logu;
	
#	ifdef DEBUG_FUN
	fputs( "<+>dgaunt()\n", debug_fp );
#	endif

	SaveTemp = phycon.te;

	/* this routine is entered with the command DRIVE GAUNT, and
	 * drives the gaunt factor routine to check range
	 * */
	fprintf( ioQQQ, " Enter 0 to input temp, energy, and net charge, or 1 for u and gamma**2.\n" );
	if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
	{
		fprintf( ioQQQ, " dgaunt error getting magic number\n" );
		puts( "[Stop in dgaunt]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	inputflag = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	
	if( inputflag == 0 )
	{
		fprintf( ioQQQ, " Enter the temperature (log if <=10), energy (Ryd), and net charge. Null line to stop.\n" );
		/* >>chng 96 july 07, got rid of statement labels replacing with do while
		 * */
		ierror = 0;
		while( ierror == 0 )
		{
			if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
			{
				fprintf( ioQQQ, " dgaunt error getting magic number\n" );
				puts( "[Stop in dgaunt]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			phycon.alogte = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* the line may be trash but ierror will pick it up  */
			if( lgEOL  )
			{
				fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
				break;
			}
			/* numbers less than or equal to 10 are the log of the temperature */
			if( phycon.alogte > 10. )
			{
				phycon.te = phycon.alogte;
				phycon.alogte = (float)log10((double)phycon.te);
			}
			else
			{
				phycon.te = (float)pow(10.,(double)phycon.alogte);
			}
			tfidle(FALSE);

			enerlin[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL || enerlin[0] == 0. )
			{
				fprintf( ioQQQ, " Sorry, but there should be two more numbers, energy and charge.\n" );
			}

			z = (double)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL || z == 0. )
			{
				fprintf( ioQQQ, " Sorry, but there should be a third number, charge.\n" );
			}

			/* This is non-thermally averaged gaunt factors.	*/
			mygaunt = CalcThermAveGaunt( (double)phycon.te, z, enerlin[0] );

			fprintf( ioQQQ, " Using my  routine, Gff= \t" );
			fprintf( ioQQQ, "%11.3e\n", mygaunt );

		}
	}
	else
	{
		/* this routine is entered with the command DRIVE GAUNT, and
		 * drives the gaunt factor routine to check range
		 * */
		fprintf( ioQQQ, " Enter log u and log gamma2. Null line to stop.\n" );
		ierror = 0;
		while( ierror == 0 )
		{
			if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
			{
				fprintf( ioQQQ, " dgaunt error getting magic number\n" );
				puts( "[Stop in dgaunt]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			logu = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* the line may be trash but ierror will pick it up  */
			if( lgEOL  )
			{
				fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
				break;
			}

			loggamma2 = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL  )
			{
				fprintf( ioQQQ, " Sorry, but there should be another numbers, log gamma2.\n" );
			}

			/* This is my attempt to calculate non-thermally averaged gaunt factors.	*/
			mygaunt = CalcThermAveGaunt( TE1RYD/pow(10.,loggamma2), 1., pow(10.,logu-loggamma2) );

			phycon.te = (float)(TE1RYD/pow(10.,loggamma2));
			tfidle(FALSE);

			fprintf( ioQQQ, " Using my  routine, Gaunt factor is" );
			fprintf( ioQQQ, "%11.3e\n", mygaunt );
		}
	}

	phycon.te = SaveTemp;
	tfidle(FALSE);

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

