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

void ParseStop(Parser &p)
{
	long int j;

	double effcol, 
	  tread;

	DEBUG_ENTRY( "ParseStop()" );

	/* time option, for stopping time dependent calculations, used to stop
	 * iterations rather than zones.  Only some stop commands have this option */
	bool lgStopZone = true;
	if( p.nMatch("TIME") )
		lgStopZone = false;

	if( p.nMatch("TEMP") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("temperature");
		}

		/* off option disables this stopping criterion */
		if( p.lgEOL() && p.nMatch(" OFF") )
		{
			/* this is special case for ending temperature - do not use -
			 * but will still stop if Te falls below TeLowest, the lowest
			 * possible temperature */
			if( lgStopZone )
				StopCalc.TempLoStopZone = -1.f;
			else
				StopCalc.TempLoStopIteration = -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. && !p.nMatch("LINE") )
			{
				tread = pow(10.,a);
			}
			else
			{
				tread = a;
			}

			/* tread is linear temperature*/
			if( tread < phycon.TEMP_LIMIT_LOW )
			{
				fprintf( ioQQQ, 
					" Temperatures below %.2e K not allowed. Reset to lowest value."
					"   I am doing this myself.\n" ,
					phycon.TEMP_LIMIT_LOW );
				/* set slightly off extreme limit for safety */
				tread = phycon.TEMP_LIMIT_LOW*1.01;
			}
			else if( tread > phycon.TEMP_LIMIT_HIGH )
			{
				fprintf( ioQQQ, 
					" Temperatures is above %.2e K not allowed. Reset to highest value."
					"   I am doing this myself.\n" ,
					phycon.TEMP_LIMIT_HIGH);
				/* set slightly off extreme limit for safety */
				tread = phycon.TEMP_LIMIT_HIGH*0.99;
			}

			if( p.nMatch("EXCE") )
			{
				/* option for this to be highest allowed temperature,
				 * stop temperate exceeds */
				if( lgStopZone )
					StopCalc.TempHiStopZone = (realnum)tread;
				else
					StopCalc.TempHiStopIteration = (realnum)tread;
			}
			else
			{
				/* this is ending temperature - we stop if kinetic temperature
				 * falls below this */
				if( lgStopZone )
					StopCalc.TempLoStopZone = (realnum)tread;
				else
					StopCalc.TempLoStopIteration = (realnum)tread;
			}
		}
	}

	/* stop at 21cm line center optical depth  */
	else if( p.nMatch("OPTI") && p.nMatch("21CM") )
	{
		/* default is for number to be log of optical depth */
		bool lgLOG = true;
		if( p.nMatch("LINE") )
		{
			/* force it to be linear not log */
			lgLOG = false;
		}
		j = (long int)p.FFmtRead();
		if( j!=21 )
		{
			fprintf( ioQQQ, " First number on STOP 21CM OPTICAL DEPTH command must be 21\n" );
			cdEXIT(EXIT_FAILURE);
		}
		/* now get the next number, which is the optical depth */
	   double a = (long int)p.FFmtRead();

		/* tau must be a log, and second number is energy where tau specified */
		if( lgLOG )
		{
			StopCalc.tauend = (realnum)pow(10.,a);
		}
		else
		{
			StopCalc.tauend = (realnum)a;
		}
		/* this flag says that 21cm line optical depth is the stop quantity */
		StopCalc.lgStop21cm = true;
	}
	/* stop optical depth at some energy */
	else if( p.nMatch("OPTI") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("optical depth");
		}

		/* default is for number to be log of optical depth */
		bool lgLOG = true;
		if( p.nMatch("LINE") )
		{
			/* force it to be linear not log */
			lgLOG = false;
		}

		if( a > 37. )
		{
			fprintf( ioQQQ, " optical depth too big\n" );
			cdEXIT(EXIT_FAILURE);
		}

		/* tau entered as a log unless liner specified */
		if( lgLOG )
		{
			StopCalc.tauend = (realnum)pow(10.,a);
		}
		else
		{
			StopCalc.tauend = (realnum)a;
		}

		/* energy where tau specified */
		StopCalc.taunu = (realnum)p.FFmtRead();

		if( p.lgEOL() )
		{
			if( p.nMatch("LYMA") )
			{
				/* largest Lyman limit optical depth */
				StopCalc.taunu = 1.;
			}
			else if( p.nMatch("BALM") )
			{
				/* 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" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		else
		{
			/* if second number is negative then log of energy in rydbergs */
			if( StopCalc.taunu < 0. )
			{
				StopCalc.taunu = (realnum)pow((realnum)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 );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* vary option */
		if( optimize.lgVarOn )
		{
			strcpy( optimize.chVarFmt[optimize.nparm], "STOP OPTICAL DEPTH = %f LOG AT %f RYD" );
			/* pointer to where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;
			optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.tauend);
			optimize.vparm[1][optimize.nparm] = (realnum)StopCalc.taunu;
			optimize.vincr[optimize.nparm] = 0.5;
			optimize.nvarxt[optimize.nparm] = 2;
			++optimize.nparm;
		}
	}

	/* stop optical depth at extinction in V filter */
	else if( p.nMatch(" AV ") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("optical depth in V");
		}
		/* 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) multiplied by scat opacity
		 * if keyword point occurs then for point source, otherwise for extended source */
		if( p.nMatch("EXTE" ) )
		{
			StopCalc.AV_extended = (realnum)a;
		}
		else
		{
			/* default is point, as measured in ism work */
			StopCalc.AV_point = (realnum)a;
		}
	}

	/* stop when a fraction of molecules frozen out on grain surfaces is reached */
	else if( p.nMatch("MOLE") && p.nMatch("DEPL") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("molecular depletion");
		}
		if( a <= 0. )
		{
			StopCalc.StopDepleteFrac = (realnum)pow(10.,a);
		}
		else
		{
			StopCalc.StopDepleteFrac = (realnum)a;
		}
	}

	/* stop when absolute value of flow velocity falls below this value */
	else if( p.nMatch("VELO") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("flow velocity");
		}
		/* entered in km/s but stored as cm/s */
		StopCalc.StopVelocity = (realnum)(a*1e5);
	}

	/* stop at a given computed mass */
	else if( p.nMatch("MASS") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("mass");
		}
		/* 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 mass */
		StopCalc.xMass = (realnum)a;
		/* NB 0 is sentinel for not set, if a is zero we must reset it */
		if( StopCalc.xMass == 0 )
			StopCalc.xMass = SMALLFLOAT;

		/* vary option */
		if( optimize.lgVarOn )
		{
			strcpy( optimize.chVarFmt[optimize.nparm], "STOP MASS = %f LOG" );
			/* pointer to where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;
			optimize.vparm[0][optimize.nparm] = (realnum)StopCalc.xMass;
			optimize.vincr[optimize.nparm] = 0.5;
			optimize.nvarxt[optimize.nparm] = 1;
			++optimize.nparm;
		}
	}

	/* 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( p.nMatch("THIC") || p.nMatch("DEPT") || p.nMatch("RADI") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("distance");
		}
		const double convl = p.nMatch("PARS") ? log10( PARSEC ) : 0.;
		bool lgStopRadius = p.nMatch("RADI") ? true : false ;
		const char* what = lgStopRadius ? "radius" : "thickness";

		if( p.nMatch("LINE") )
		{
			if( a > 0. )
			{
				a = log10(a) + convl;
			}
			else
			{
				fprintf(ioQQQ,"The %s is negative and linear is set - this is impossible.\n", what);
				cdEXIT(EXIT_FAILURE);
			}
		}
		else
		{
			a += convl;
		}
		if( a > 37. )
		{
			fprintf( ioQQQ, "DISASTER %s too large\n", what );
			cdEXIT(EXIT_FAILURE);
		}
		if( lgStopRadius )
			radius.StopRadius[0] = pow(10.,a);
		else
			radius.StopThickness[0] = pow(10.,a);

		/* can stop at different thickness on each iteration */
		for( j=1; j < iterations.iter_malloc; j++ )
		{
			a = p.FFmtRead();
			if( p.lgEOL() )
			{
				if( lgStopRadius )
					radius.StopRadius[j] = radius.StopRadius[j-1];
				else
					radius.StopThickness[j] = radius.StopThickness[j-1];
			}
			else
			{
				if( p.nMatch("LINE") )
				{
					if( a > 0. )
					{
						a = log10(a) + convl;
					}
					else
					{
						fprintf(ioQQQ,"The %s is negative and linear is set -"
							" this is impossible.\n", what);
						cdEXIT(EXIT_FAILURE);
					}
				}
				else
				{
					a += convl;
				}
				if( a > 37. )
				{
					fprintf( ioQQQ, "DISASTER %s too large\n", what );
					cdEXIT(EXIT_FAILURE);
				}
				if( lgStopRadius )
					radius.StopRadius[j] = pow(10.,a);
				else
					radius.StopThickness[j] = pow(10.,a);
			}
		}

		/* vary option */
		if( optimize.lgVarOn )
		{
			optimize.nvarxt[optimize.nparm] = 1;
			/* pointer to where to write */
			optimize.nvfpnt[optimize.nparm] = input.nRead;
			if( lgStopRadius )
			{
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP RADIUS %f LOG" );
				optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.StopRadius[0]);
			}
			else
			{
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP THICKNESS %f LOG" );
				optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.StopThickness[0]);
			}
			optimize.vincr[optimize.nparm] = 0.5f;
			++optimize.nparm;
		}
	}

	/* stop at a particular zone, for each iteration */
	else if( p.nMatch("ZONE") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("zone number");
		}
		/* 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;

		long int nZoneMax = geometry.nend[0];
		for( j=1; j < iterations.iter_malloc; j++ )
		{
			geometry.nend[j] = (long)p.FFmtRead();
			/* if eol on this iteration, set to previous.  In most cases
			 * all will be equal to the first */
			if( p.lgEOL() )
			{
				geometry.nend[j] = geometry.nend[j-1];
			}
			else
			{
				/* 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] );
			}
			nZoneMax = max( nZoneMax , geometry.nend[j] );
		}

		if( nZoneMax>2000 )
			fprintf(ioQQQ,"CAUTION - it will take a lot of memory to save"
			" results for %li zones.  Is this many zones really necessary?\n",
			nZoneMax );
	}

	/* stop when a prescribed continuum flux is reached */
	else if( p.nMatch("CONT") && p.nMatch("FLUX") )
	{
		/* first read the continuum energy and add this point to PredCont */
		double energy = p.FFmtRead();
		if( p.lgEOL() )
			p.NoNumb("energy");
		const char* unit = p.StandardEnergyUnit();
		long ind = t_PredCont::Inst().add( energy, unit );
		Energy E( energy, unit );

		double flux = p.FFmtRead();
		if( p.lgEOL() )
			p.NoNumb("flux");
		if( flux <= 0. || p.nMatch( " LOG") )
			flux = pow(10.,flux);
		Flux F( E, flux, p.StandardFluxUnit() );

		StopCalc.ContIndex.push_back( ind );
		StopCalc.ContNFnu.push_back( F );
	}

	/* stop at this electron fraction, relative to hydrogen */
	else if( p.nMatch("EFRA") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("electron fraction");
		}
		if( a <= 0. )
		{
			StopCalc.StopElecFrac = (realnum)pow(10.,a);
		}
		else
		{
			StopCalc.StopElecFrac = (realnum)a;
		}
	}

	/* stop at a hydrogen molecular fraction, relative to total hydrogen,
	 * this is 2H_2 / H_total*/
	else if( p.nMatch("MFRA") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("hydrogen molecular fraction");
		}
		if( a <= 0. )
		{
			StopCalc.StopH2MoleFrac = (realnum)pow(10.,a);
		}
		else
		{
			StopCalc.StopH2MoleFrac = (realnum)a;
		}
	}

	/* stop at a ionized hydrogen fraction, relative to total hydrogen,
	 * this is H+ / H_total */
	else if( p.nMatch("PFRA") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("ionized hydrogen fraction");
		}
		if( a <= 0. )
		{
			StopCalc.StopHPlusFrac = (realnum)pow(10.,a);
		}
		else
		{
			StopCalc.StopHPlusFrac = (realnum)a;
		}
	}

	/* stop at a particular column density */
	else if( p.nMatch("COLU") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("column density");
		}
		char chLabel[CHARS_SPECIES];

		/* check for linear option, if present take log since a being
		 * log column density is default */
		if( p.nMatch( "LINE" ) )
			a = log10(a);

		if( !p.GetQuote( chLabel , false ) )
		{
			/* species was given in double quotes */
			StopCalc.col_species = (realnum)pow(10.,a);
			StopCalc.lgStopSpeciesColumn = true;
			strncpy( StopCalc.chSpeciesColumn, chLabel, CHARS_SPECIES );
		}
		/*  stop at an effective column density */
		else if( p.nMatch("EFFE") )
		{
			/* actually stop at certain optical depth at 1keV */
			effcol = pow(10.,a);
			StopCalc.tauend = (realnum)(effcol*2.14e-22);
			StopCalc.taunu = (realnum)(1000./EVRYD);
			/* vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP EFFECTIVE COLUMN DENSITY %f LOG" );
				/* pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/* log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (realnum)log10(effcol);
				optimize.vincr[optimize.nparm] = 0.5f;
				++optimize.nparm;
			}
		}

		else if( p.nMatch("IONI") )
		{
			/*  this is ionized column */
			if( a > 37. )
			{
				fprintf( ioQQQ, " column too big\n" );
				cdEXIT(EXIT_FAILURE);
			}

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

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

		/*  stop at a neutral column */
		else if( p.nMatch("NEUT") )
		{
			StopCalc.colnut = (realnum)pow(10.,a);

			/*  vary option */
			if( optimize.lgVarOn )
			{
				optimize.nvarxt[optimize.nparm] = 1;
				strcpy( optimize.chVarFmt[optimize.nparm], "STOP NEUTRAL COLUMN DENSITY %f LOG");
				/*  pointer to where to write */
				optimize.nvfpnt[optimize.nparm] = input.nRead;
				/*  log of temp will be pointer */
				optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.colnut);
				optimize.vincr[optimize.nparm] = 0.5f;
				++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( p.nMatch(" H2 ") )
		{
			/* 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 */
			j = (long int)a;
			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" );
				cdEXIT(EXIT_FAILURE);
			}
			a = p.FFmtRead();
			StopCalc.col_h2 = (realnum)pow(10.,a);

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

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

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

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

		/* fall through default is total hydrogen column density */
		else
		{
			/*  both HII and HI */
			if( a > 37. )
			{
				fprintf( ioQQQ, " column too big\n" );
				cdEXIT(EXIT_FAILURE);
			}

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

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

	/* stop when electron density falls below this value, linear or log */
	else if( p.nMatch("EDEN") )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("electron density");
		}
		/* stop if electron density falls below this value
		 * LINEAR option */
		if( p.nMatch("LINE") )
		{
			StopCalc.StopElecDensity = (realnum)a;
		}
		else
		{
			StopCalc.StopElecDensity = (realnum)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( p.nMatch("LINE") )
	{
		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 */
		p.GetQuote( chLabel , true );

		/* copy first four char of label into caps, and null terminate*/
		strncpy( StopCalc.chStopLabel1[StopCalc.nstpl], chLabel , 4 );
		StopCalc.chStopLabel1[StopCalc.nstpl][4] = 0;

		// default intrinsic intensity, accept emergent
		StopCalc.nEmergent[StopCalc.nstpl] = 0;
		if( p.nMatch("EMER") )
			StopCalc.nEmergent[StopCalc.nstpl] = 1;

		/* get line wavelength */
		StopCalc.StopLineWl1[StopCalc.nstpl] = (realnum)p.getWaveOpt();

		/* get relative intensity */
		StopCalc.stpint[StopCalc.nstpl] = (realnum)p.FFmtRead();
		if( p.lgEOL() )
		{
			fprintf( ioQQQ, " There MUST be a relative intensity  entered "
				"for first line in STOP LINE command.  Sorry.\n" );
			cdEXIT(EXIT_FAILURE);
		}

		/* check for second line - use Hbeta is not specified */

		/* get label for the line - must do this first so we clear the label string before
		 * trying to read the wavelength */
		
		j = p.GetQuote( chLabel , false );

		if( j != 0 )
		{
			/* no line label, normalization line will be H beta */
			strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], "TOTL" , 4 );
			StopCalc.chStopLabel2[StopCalc.nstpl][4] = 0;
			StopCalc.StopLineWl2[StopCalc.nstpl] = 4861.f;
		}
		else
		{
			/* copy first four char of label into caps, and null terminate*/
			strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], chLabel , 4 );
			StopCalc.chStopLabel2[StopCalc.nstpl][4] = 0;

			/* wavelength of second line, may be absent and so zero -
			 * we will use Hbeta if not specified */
			StopCalc.StopLineWl2[StopCalc.nstpl] = (realnum)p.getWaveOpt();

		}
		/* increment number of stop lines commands entered */
		StopCalc.nstpl = MIN2(StopCalc.nstpl+1,MXSTPL-1);
	}

	else if( p.nMatch("NTOTALIO" ) )
	{
		double a = p.FFmtRead();
		
		if( p.lgEOL() && !p.nMatch(" OFF") )
		{
			p.NoNumb("number of calls to conv_base");
		}
		/* this counter is incremented at end of conv_base - total number of
		 * times conv base was called 
		 * this is a debugging aid - code aborts */
		StopCalc.nTotalIonizStop = (long)a;
	}

	/* oops! no keyword that we could find */
	else
	{
		fprintf( ioQQQ, " I did not recognize a keyword on this STOP line, line image follows;\n" );
		p.PrintLine(ioQQQ);
		fprintf( ioQQQ, "Sorry.\n");
		cdEXIT(EXIT_FAILURE);
	}
	return;
}