/* 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 */
/*ParseStop parse the stop command */
#include "cddefines.h"
#include "physconst.h"
#include "optimize.h"
#include "rfield.h"
#include "radius.h"
#include "geometry.h"
#include "stopcalc.h"
#include "input.h"
#include "parse.h"

void ParseStop(char *chCard)
{
	int lgEOL;

	long int i, 
	  j;

	double a, 
	  effcol, 
	  tread;

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

	/* first scan off a generic number, all subcommands use at least one */
	i = 5;
	a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* >>chng 04 nov 09, introduce off option */
	if( lgEOL && !lgMatch(" OFF",chCard) )
	{
		NoNumb(chCard);
	}

	if( lgMatch("TEMP",chCard) )
	{
		/* >>chng 04 nov 09, introduce off option to disable this stopping criterion */
		if( lgEOL && lgMatch(" OFF",chCard) )
		{
			/* this is ending temperature */
			StopCalc.tend = -1.f;
			/* >>chng 05 aug 15, do not reset StopCalc.TeLowest since this is lowest
			 * possible temperature code can handle, set to 2.7 K */
			/*StopCalc.TeLowest = -1.f;*/
		}
		else
		{
			/* lowest electron temperature allowed before stopping
			 * assumed to be the log of the temperature if <10
			 * optional keyword LINEAR forces linear */
			if( a <= 10. && !lgMatch("LINE",chCard) )
			{
				tread = pow(10.,a);
			}
			else
			{
				tread = a;
			}

			/* tread is linear temperature*/
			if( tread < 3. )
			{
				fprintf( ioQQQ, " Temperatures below 3K not allowed. Reset to 3K.  I am doing this myself.\n" );
				tread = 3.;
			}

			if( lgMatch("EXCE",chCard) )
			{
				/* option for this to be highest allowed temperature,
				 * stop temperate exceeds */
				StopCalc.T2High = (float)tread;
			}
			else
			{
				/* this is ending temperature */
				StopCalc.tend = (float)tread;
				/* >>chng 05 aug 15, do not set TeLowest - this is lower temp bound,
				 * set to 2.7K - cannot go below this */
				/*StopCalc.TeLowest = (float)MIN2(StopCalc.TeLowest,StopCalc.tend);*/
			}
		}
	}

	/* stop at 21cm line center optical depth  */
	else if( lgMatch("OPTI",chCard) && lgMatch("21CM",chCard) )
	{
		/* default is for number to be log of optical depth */
		int lgLOG = TRUE;
		if( lgMatch("LINE",chCard) )
		{
			/* force it to be linear not log */
			lgLOG = FALSE;
		}
		i = 4;
		j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( j!=21 )
		{
			fprintf( ioQQQ, " First number on STOP 21CM OPTICAL DEPTH command must be 21\n" );
			puts( "[Stop  in ParseStop]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* now get the next number, which is the optical depth */
		a = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
 
		/* tau must be a log, and second number is energy where tau specified */
		if( lgLOG )
		{
			StopCalc.tauend = (float)pow(10.,a);
		}
		else
		{
			StopCalc.tauend = (float)a;
		}
		/* this flag says that 21cm line optical depth is the stop quantity */
		StopCalc.lgStop21cm = TRUE;
	}
	/* stop optical depth at some energy */
	else if( lgMatch("OPTI",chCard) )
	{
		/* default is for number to be log of optical depth */
		int lgLOG = TRUE;
		if( lgMatch("LINE",chCard) )
		{
			/* force it to be linear not log */
			lgLOG = FALSE;
		}

		if( a > 37. && lgBit32 )
		{
			fprintf( ioQQQ, " optical depth too big for a 32 bit cpu?\n" );
			puts( "[Stop  in ParseStop]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* tau must be a log, and second number is energy where tau specified */
		if( lgLOG )
		{
			StopCalc.tauend = (float)pow(10.,a);
		}
		else
		{
			StopCalc.tauend = (float)a;
		}

		/* energy where tau specified */
		StopCalc.taunu = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		if( lgEOL )
		{
			if( lgMatch("LYMA",chCard) )
			{
				/* largest Lyman limit optical depth */
				StopCalc.taunu = 1.;
			}
			else if( lgMatch("BALM",chCard) )
			{
				/* stop at this balmer continuum optical depth */
				StopCalc.taunu = 0.25;
			}
			else
			{
				fprintf( ioQQQ, " There must be a second number, the energy in Ryd.  Sorry.\n" );
				puts( "[Stop  in ParseStop]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else
		{
			/* if second number is negative then log of energy in rydbergs */
			if( StopCalc.taunu < 0. )
			{
				StopCalc.taunu = (float)pow(10.f,StopCalc.taunu);
			}

			/* check that energy is within bounds of code */
			if( StopCalc.taunu < rfield.emm || StopCalc.taunu > rfield.egamry )
			{
				fprintf( ioQQQ, " The energy must be in the range %10.2e to %10.2e.  It was %10.2e. Sorry.\n", 
				  rfield.emm, rfield.egamry, StopCalc.taunu );
				puts( "[Stop  in ParseStop]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	/* stop optical depth at extinction in V filter */
	else if( lgMatch(" AV ",chCard) )
	{
		/* default is for number to be A_V, log if negative */
		if( a<=0. )
		{
			a = pow(10.,a);
		}
		/* A_V can be for either point or extended source, difference is (1-g) mult scat opacity
		 * if keyword point occurs then for point source, otherwise for extended source */
		if( lgMatch("EXTE" , chCard ) )
		{
			StopCalc.AV_extended = (float)a;
		}
		else
		{
			/* default is point, as measured in ism work */
			StopCalc.AV_point = (float)a;
		}
	}

	/* stop when a fraction of molecules frozen out on grain surfaces is reached */
	else if( lgMatch("MOLE",chCard) && lgMatch("DEPL",chCard) )
	{
		if( a <= 0. )
		{
			StopCalc.StopDepleteFrac = (float)pow(10.,a);
		}
		else
		{
			StopCalc.StopDepleteFrac = (float)a;
		}
	}

	/* stop at a given computed mass */
	else if( lgMatch("MASS",chCard) )
	{
		/* number of log of mass in gm if inner radius is specified,
		 * mass per unit area, gm cm-2 if not
		 * leave it as a log since dare not deal with linear cm on 32-bit processor */
		StopCalc.xMass = (float)a;
		/* NB 0 is sentinal for not set, if a is zero we must reset it */
		if( StopCalc.xMass == 0 )
			StopCalc.xMass = SMALLFLOAT;
	}

	/* stop thickness command, also stop depth, this must come after stop
	 * optical depth, since do not want to trigger on depth in optical depth */
	else if( lgMatch("THIC",chCard) || lgMatch("DEPT",chCard) )
	{
		/* decide whether linear or log,
		 * this branch if key linear appears */
		if( lgMatch("LINE",chCard) )
		{
			radius.router[0] = a;
		}
		else
		{
			/* this is default branch, no key, so log */
			if( a > 37. && lgBit32 )
			{
				fprintf( ioQQQ, " thickness too large for a 32 bit cpu?\n" );
				puts( "[Stop  in ParseStop]" );
				cdEXIT(EXIT_FAILURE);
			}
			radius.router[0] = pow(10.,a);
		}

		/* optional parsec unit */
		if( lgMatch("PARS",chCard) )
		{
			radius.router[0] *= PARSEC;
		}

		/* can stop at different thickness on each iteration */
		for( j=1; j < ITRDIM; j++ )
		{
			a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				radius.router[j] = radius.router[j-1];
			}
			else
			{
				if( lgMatch("LINE",chCard) )
				{
					radius.router[j] = a;
				}
				else
				{
					if( a > 37. && lgBit32 )
					{
						fprintf( ioQQQ, " thickness too large for a 32 bit cpu?\n" );
						puts( "[Stop  in ParseStop]" );
						cdEXIT(EXIT_FAILURE);
					}
					radius.router[j] = pow(10.,a);
				}
				if( lgMatch("PARS",chCard) )
				{
					radius.router[j] *= PARSEC;
				}
			}
		}

		/* vary option */
		if( optimize.lgVarOn )
		{
			optimize.nvarxt[optimize.nparm] = 1;
			strcpy( optimize.chVarFmt[optimize.nparm], "STOP THICKNESS %f" );

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

			/* log of temp will be pointer */
			optimize.vparm[0][optimize.nparm] = (float)log10(radius.router[0]);
			optimize.vincr[optimize.nparm] = 0.5;

			++optimize.nparm;
		}
	}

	/* stop at a particular zone, for each iteration */
	else if( lgMatch("ZONE",chCard) )
	{
		/* stop after computing this zone */
		/* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
		 * as stop zone 1, bug caught by Joop Schaye */
		geometry.nend[0] = (long)MAX2(1.,a);
		geometry.lgZoneSet = TRUE;

		/* this tells code that we intend to stop at this zone, so caution not generated*/
		geometry.lgEndDflt = FALSE;

		for( j=1; j < ITRDIM; j++ )
		{
			geometry.nend[j] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* if eol on this iteration, set to previous.  In most cases
			 * all will be equal to the first */
			if( lgEOL )
			{
				geometry.nend[j] = geometry.nend[j-1];
			}
			else
			{
				/* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
				 * as stop zone 1, bug caught by Joop Schaye */
				geometry.nend[j] = MAX2( 1 , geometry.nend[j] );
			}
		}
	}

	/* stop at this electron fraction, relative to hydrogen */
	else if( lgMatch("EFRA",chCard) )
	{
		if( a <= 0. )
		{
			StopCalc.StopElecFrac = (float)pow(10.,a);
		}
		else
		{
			StopCalc.StopElecFrac = (float)a;
		}
	}

	/* stop at a hydrogen molecular fraction, relative to total hydrogen,
	 * this is 2H_2 / H_total*/
	else if( lgMatch("MFRA",chCard) )
	{
		if( a <= 0. )
		{
			StopCalc.StopH2MoleFrac = (float)pow(10.,a);
		}
		else
		{
			StopCalc.StopH2MoleFrac = (float)a;
		}
	}

	/* stop at a ionized hydrogen fraction, relative to total hydrogen,
	 * this is H+ / H_total */
	else if( lgMatch("PFRA",chCard) )
	{
		if( a <= 0. )
		{
			StopCalc.StopHPlusFrac = (float)pow(10.,a);
		}
		else
		{
			StopCalc.StopHPlusFrac = (float)a;
		}
	}

	/* stop at a particular column density */
	else if( lgMatch("COLU",chCard) )
	{
		/*  stop at an effective column density */
		if( lgMatch("EFFE",chCard) )
		{
			/* actually stop at certain optical depth at 1keV */
			effcol = pow(10.,a);
			StopCalc.tauend = (float)(effcol*2.14e-22);
			StopCalc.taunu = 73.5;
			/* vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP EFFECTIVE COLUMN DENSITY %f" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/* log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(effcol);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		else if( lgMatch("IONI",chCard) )
		{
			/*  this is ionized column */
			if( a > 37. && lgBit32 )
			{
				fprintf( ioQQQ, " column too big for a 32 bit cpu?\n" );
				puts( "[Stop  in ParseStop]" );
				cdEXIT(EXIT_FAILURE);
			}

			StopCalc.colpls = (float)pow(10.,a);

			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP IONIZED COLUMN DENSITY %f" );
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.colpls);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		/*  stop at a neutral column */
		else if( lgMatch("NEUT",chCard) )
		{
			StopCalc.colnut = (float)pow(10.,a);
			
			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP NEUTRAL COLUMN DENSITY %f");
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.colnut);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		/* >>chng 03 apr 15, add this option 
		 *  stop at a molecular hydrogen column density, input parameter
		 * is log of the H2 column density */
		else if( lgMatch(" H2 ",chCard) )
		{
			/* this command has a 2 in the H2 label - must not parse the two by
			 * accident.  Get the first number off the line image, and confirm that
			 * it is a 2 */
			i = 5;
			j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( j != 2 )
			{
				fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
				fprintf( ioQQQ, " The first number I encounter should be the 2 in H2.\n Sorry.\n" );
				puts( "[Stop in ParseStop]" );
				cdEXIT(EXIT_FAILURE);
			}
			a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			StopCalc.col_h2 = (float)pow(10.,a);
		
			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP H2 COLUMN DENSITY %f");
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_h2);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		else if( lgMatch("ATOM",chCard) )
		{
			StopCalc.col_h2_nut = (float)pow(10.,a);
			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP ATOMIC COLUMN DENSITY %f");
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_h2_nut);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		else if( lgMatch("H/TS",chCard) )
		{
			/* >> 05 jan 09, add stop integrated n(H0) / Tspin */
			StopCalc.col_H0_ov_Tspin = (float)pow(10.,a);
			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP H/TSPIN COLUMN DENSITY %f");
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_H0_ov_Tspin);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		else if( lgMatch(" CO ",chCard) )
		{
			/* chng << 03 Oct. 27--Nick Abel, add this option */
			/* stop at a carbon monoxide column density */
			StopCalc.col_monoxco = (float)pow(10.,a);
			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP CO COLUMN DENSITY %f");
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_monoxco);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}

		/* fall through default is total hydrogen column density */
		else
		{
			/*  both HII and HI */
			if( a > 37. && lgBit32 )
			{
				fprintf( ioQQQ, " column too big for a 32 bit cpu?\n" );
				puts( "[Stop  in ParseStop]" );
				cdEXIT(EXIT_FAILURE);
			}

			StopCalc.HColStop = (float)pow(10.,a);

			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP COLUMN DENSITY %f" );
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.HColStop);
				optimize.vincr[optimize.nparm] = 0.5;
				++optimize.nparm;
			}
		}
	}

	/* stop when electron density falls below this value, linear or log */
	else if( lgMatch("EDEN",chCard) )
	{
		/* stop if electron density falls below this value
		 * LINEAR option */
		if( lgMatch("LINE",chCard) )
		{
			StopCalc.StopElecDensity = (float)a;
		}
		else
		{
			StopCalc.StopElecDensity = (float)pow(10.,a);
		}
	}

	/* stop at a particular line ratio - this must come last since many commands
	 * have linear option - don't want to trigger on that */
	else if( lgMatch("LINE",chCard) )
	{
		char chLabel[5];
		/* first line wavelength, then intensity relative to Hbeta for stop
		 * if third number is entered, it is wl of line in denominator */

		/* get label for the line - must do this first so we clear the label string before
		 * trying to read the wavelength */
		GetQuote( chLabel , chCard , TRUE );

		/* copy first four char of label into caps, and null terminate*/
		cap4( StopCalc.chStopLabel[StopCalc.nstpl], chLabel);

		/* now must rescan the line since label was probably picked in in the first read */

		i = 5;
		/* save line wavelength */
		StopCalc.StopLineWl[StopCalc.nstpl] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);

		/* line was entered, look for possible micron or cm label */
		if( input.chCARDCAPS[i-1] == 'M' )
		{
			/* microns */
			StopCalc.StopLineWl[StopCalc.nstpl] *= 1e4;
		}
		else if( input.chCARDCAPS[i-1] == 'C' )
		{
			/* microns */
			StopCalc.StopLineWl[StopCalc.nstpl] *= 1e8;
		}

		/* get the error assocated with 4 significant figures */
		a = log10( StopCalc.StopLineWl[StopCalc.nstpl] );
		if( StopCalc.StopLineWl[StopCalc.nstpl]>=1. )
		{
			a -= floor(a);
		}
		else
		{
			a += floor(a);
		}
		a = pow(10.,a);
		/* a is now 1 - 10, so this corresponds to 4 sig fig */
		StopCalc.ErrorLineWl[StopCalc.nstpl] = (float)(0.00105/a);

		/* get relative intensity */
		StopCalc.stpint[StopCalc.nstpl] = 
			(float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);

		if( lgEOL )
		{
			fprintf( ioQQQ, " There MUST be a relative intensity  entered.  Sorry.\n" );
			puts( "[Stop  in ParseStop]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* check for second line - may not be there - save index for scanner
		 * and check if any more numbers on the line */

		j = i;
		a = FFmtRead(chCard,&j,INPUT_LINE_LENGTH, &lgEOL);

		if( lgEOL )
		{
			/* hit the eol, no more numbers */
			StopCalc.StopLineWl2[StopCalc.nstpl] = 1;
		}
		else
		{
			/* get label for the line - must do this first so we clear the label string before
			 * trying to read the wavelength */
			GetQuote( chLabel , chCard , TRUE );

			/* copy first four char of label into caps, and null terminate*/
			cap4( StopCalc.chStopLabel2[StopCalc.nstpl], chLabel);

			/* wavelength of second line, may be absent and so zero -
			 * we will use Hbeta if not specified */
			StopCalc.StopLineWl2[StopCalc.nstpl] = 
				(float)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);

			/* line was entered, look for possible micron or cm label */
			if( input.chCARDCAPS[i-1] == 'M' )
			{
				/* microns */
				StopCalc.StopLineWl2[StopCalc.nstpl] *= 1e4;
			}
			else if( input.chCARDCAPS[i-1] == 'C' )
			{
				/* microns */
				StopCalc.StopLineWl2[StopCalc.nstpl] *= 1e8;
			}

			/* get the error assocated with 4 significant figures */
			if( StopCalc.StopLineWl2[StopCalc.nstpl] > 0. )
			{
				a = log10( StopCalc.StopLineWl2[StopCalc.nstpl] );
				if( StopCalc.StopLineWl2[StopCalc.nstpl]>=1. )
				{
					a -= floor(a);
				}
				else
				{
					a += floor(a);
				}
				a = pow(10.,a);
				/* a is now 1 - 10, so this corresponds to 4 sig fig */
				StopCalc.ErrorLineWl2[StopCalc.nstpl] = (float)(0.00105/a);
			}
			else
			{
				StopCalc.ErrorLineWl2[StopCalc.nstpl] = 1e-4f;
			}
		}
		/* increment number of stop lines intered, this number is used as test
		 * for whether any stop line ever set */
		StopCalc.nstpl = MIN2(StopCalc.nstpl+1,MXSTPL-1);
	}

	/* oops! no keyword that we could find */
	else
	{
		fprintf( ioQQQ, " keyword error for STOP line, line image follows;\n" );
		fprintf( ioQQQ, "==%s==\n" , chCard);
		puts( "[Stop  in ParseStop]" );
		cdEXIT(EXIT_FAILURE);
	}

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

