/* 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 */
/*cdDrive main routine to call cloudy under all circumstances) */
/*cdReasonGeo wrte why the model stopped and type of geometry on io file */
/*cdWarnings write all warnings entered into comment stack */
/*cdEms obtain the local emissivity for a line, for the last computed zone */
/*cdColm get the column density for a constituent  */
/*cdLine get the predicted line intensity, also index for line in stack */
/*cdLine_ip get the predicted line intensity, using index for line in stack */
/*cdDLine get the predicted emergent line intensity, also index for line in stack */
/*cdCautions print out all cautions after calculation, on arbitrary io unit */
/*cdGetTemp_last routine to query results and return temperature of last zone */
/*cdGetDepth_depth get depth structure from previous iteration */
/*cdTimescales returns thermal, recombination, and H2 foramtion timescales */
/*cdSurprises print out all surprises on arbitrary unit number */
/*cdNotes print stack of notes about current calculation */
/*cdGetPressure_last routine to query results and return pressure of last zone */
/*cdTalk tells the code whether to print results or be silent */
/*cdOutp redirect output to arbitrary Fortran unit number */
/*cdRead routine to read in command lines when cloudy used as subroutine */
/*cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit */
/*cdIonFrac get ionization fractions for a constituent */
/*cdTemp get mean electron temperature for any element */
/*cdGetCooling_last routine to query results and return cooling of last zone */
/*cdGetHeating_last routine to query results and return heating of last zone */
/*cdNoExec call this routine to tell code not to actually execute */
/*cdDate - puts date of code into string */
/*cdVersion produces string that gives version number of the code */
/*cdExecTime any routine can call this, and find the time since cdInit was called */
/*cdPrintCommands( FILE *) prints all input commands into file */
/*cdDrive main routine to call cloudy under all circumstances) */
/*cdNwcns get the number of cautions and warnings, to tell if calculation is ok */
/*debugLine provides a debugging hook into the main line array  */
/*cdPrtWL print line wavelengths in Angstroms in the standard format - just a wrapper */
/*cdEms_ip obtain the local emissivity for a line with known index */
/*cdGetnZone gets number of zones */
/*cdClosePunchFiles closes all the punch files that have been used */
/*cdLineListPunch create a file with a list of all emission lines punched, 
 *and their index within the emission line stack */
#include "cddefines.h"
#include "trace.h"
#include "converge.h"
#include "lines.h"
#include "pressure.h"
#include "prt.h"
#include "colden.h"
#include "radius.h"
#include "struc.h"
#include "co.h"
#include "elementnames.h"
#include "mean.h"
#include "phycon.h"
#include "called.h"
#include "parse.h"
#include "input.h"
#include "taulines.h"
#include "version.h"
#include "thermal.h"
#include "optimize.h"
#include "timesc.h"
#include "cloudy.h"
#include "warnings.h"
#include "cddrive.h"
#include "lines_service.h"

/*************************************************************************
 *
 * cdDrive - main routine to call cloudy - returns 0 if all ok, 1 if problems
 *
 ************************************************************************/

int cdDrive(void )
{
	int lgBAD;

#	ifdef DEBUG_FUN
	fputs( "<+>cdDrive()\n", debug_fp );
#	endif
	/*********************************
	 * main routine to call cloudy   *
	 *********************************/

	/* this is set FALSE when code loaded, set TRUE when cdInit called,
	 * this is check that cdInit was called first */
	if( !lgcdInitCalled )
	{
		printf(" cdInit was not called first - this must be the first call.\n");
		puts( "[Stop in cdDrive]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( trace.lgTraceInput )
	{
		fprintf( ioQQQ, 
			"cdDrive: lgOptimr=%1i lgVaryOn=%1i lgNoVary=%1i input.nSave:%li\n",
			optimize.lgOptimr , optimize.lgVaryOn , optimize.lgNoVary, input.nSave );
	}

	/* should we call cloudy, or the optimization driver? */
	/* possible to have VARY on line without OPTIMIZE being set */
	if( optimize.lgOptimr && optimize.lgVaryOn )
	{
		optimize.lgVaryOn = TRUE;
	}
	else
	{
		optimize.lgVaryOn = FALSE;
	}

	/* this is option to have command line saying "no optimize"
	 * to turn off optimizer */
	if( optimize.lgNoVary )
		optimize.lgVaryOn = FALSE;

	if( optimize.lgVaryOn )
	{
		if( trace.lgTraceInput )
			fprintf( ioQQQ, "cdDrive: calling optimize_do\n");
		/* option to drive optimizer set if OPTIMIZE was in input stream */
		lgBAD = optimize_do();
	}
	else
	{
		if( trace.lgTraceInput )
			fprintf( ioQQQ, "cdDrive: calling cloudy\n");
		/* optimize did not occur, only compute one model, call cloudy */
		lgBAD = cloudy();
	}

	/* reset flag saying that cdInit has not been called */
	lgcdInitCalled = FALSE;

	if( conv.lgAbort || lgBAD )
	{
		if( trace.lgTraceInput )
			fprintf( ioQQQ, "cdDrive: returning failure during call. \n");
		/* lgAbort set true if something wrong, so return lgBAD false. */
#		ifdef DEBUG_FUN
		fputs( " <->cdDrive()\n", debug_fp );
#		endif
		return(1);
	}
	else
	{
		/* everything is ok, return 0 */
#		ifdef DEBUG_FUN
		fputs( " <->cdDrive()\n", debug_fp );
#		endif
		return(0);
	}
}

/*************************************************************************
 *
 * cdReasonGeo wrte why the model stopped and type of geometry on io file 
 *
 ************************************************************************/


/*cdReasonGeo wrte why the model stopped and type of geometry on io file */
void cdReasonGeo(FILE * ioOUT)
{

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

	/*this is the reason the calculation stopped*/
	fprintf( ioOUT, "%s", warnings.chRgcln[0] );
	fprintf( ioOUT , "\n" );
	/* this is the geometry */
	fprintf( ioOUT, "%s", warnings.chRgcln[1] );
	fprintf( ioOUT , "\n" );

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


/*************************************************************************
 *
 * cdWarnings write all warnings entered into comment stack 
 *
 ************************************************************************/

/*cdWarnings write all warnings entered into comment stack */

void cdWarnings(FILE *ioPNT )
{
	long int i;

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

	if( warnings.nwarn > 0 )
		{
			for( i=0; i < warnings.nwarn; i++ )
			{
				/* these are all warnings that were entered in comment */
				fprintf( ioPNT, "%s", warnings.chWarnln[i] );
				fprintf( ioPNT, "\n" );
			}
		}


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


/*************************************************************************
 *
 * cdCautions print out all cautions after calculation, on arbitrary io unit 
 *
 ************************************************************************/

/*cdCautions print out all cautions after calculation, on arbitrary io unit */

void cdCautions(FILE * ioOUT)
{
	long int i;

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

	if( warnings.ncaun > 0 )
	{
		for( i=0; i < warnings.ncaun; i++ )
		{
			fprintf( ioOUT, "%s", warnings.chCaunln[i] );
			fprintf( ioOUT, "\n" );
		}
	}

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

/*************************************************************************
 *
 * cdTimescales returns thermal, recombination, and H2 foramtion timescales 
 *
 ************************************************************************/

void cdTimescales(
	/* the thermal cooling timescale */
	double *TTherm , 
	/* the hydrogen recombination timescale */
	double *THRecom , 
	/* the H2 formation timescale */
	double *TH2 )
{

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

	/* these were all evaluated in AgeCheck, which was called by PrtComment */

	/* thermal or cooling timescale */
	*TTherm = timesc.ttherm;

	/* the hydrogen recombination timescale */
	*THRecom = timesc.threc;

	/* longer of the the H2 formation and destruction timescales */
	*TH2 = MAX2( timesc.time_H2_Dest_longest , timesc.time_H2_Form_longest );

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


/*************************************************************************
 *
 * cdSurprises print out all surprises on arbitrary unit number 
 *
 ************************************************************************/

/*cdSurprises print out all surprises on arbitrary unit number */

void cdSurprises(FILE * ioOUT)
{
	long int i;

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

	if( warnings.nbang > 0 )
	{
		for( i=0; i < warnings.nbang; i++ )
		{
			fprintf( ioOUT, "%s", warnings.chBangln[i] );
			fprintf( ioOUT, "\n" );
		}
	}


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


/*************************************************************************
 *
 * cdNotes print stack of notes about current calculation 
 *
 ************************************************************************/

/*cdNotes print stack of notes about current calculation */

void cdNotes(FILE * ioOUT)
{
	long int i;

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

	if( warnings.nnote > 0 )
		{
			for( i=0; i < warnings.nnote; i++ )
			{
				fprintf( ioOUT, "%s", warnings.chNoteln[i] );
				fprintf( ioOUT, "\n" );
			}
		}

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

/*************************************************************************
 *
 * cdGetCooling_last routine to query results and return cooling of last zone 
 *
 ************************************************************************/

/*cdGetCooling_last routine to query results and return cooling of last zone */
double cdGetCooling_last(void) /* return cooling for last zone */
{
	return (thermal.ctot);
}

/*************************************************************************
 *
 * cdVersion - puts version number of code into string 
 * incoming string must have at least 8 char and will become null
 * terminated string
 *
 ************************************************************************/

void cdVersion(char chString[] ) 
{
	if( !lgcdInitCalled )
	{
		printf(" cdInit was not called first - this must be the first call.\n");
		puts( "[Stop in cdVersion]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* version was set by cdInit */
	strcpy( chString , version.chVersion );
	return ;
}

/*************************************************************************
 *
 * cdDate - puts date of code into string 
 * incoming string must have at least 8 char and will become null
 * terminated string
 *
 ************************************************************************/

/* cdDate - puts date of code into string  */
void cdDate(char chString[] ) 
{
	/* date was set by cdInit */
	strcpy( chString , version.chDate );
	return ;
}


/*************************************************************************
 *
 * cdGetHeating_last routine to query results and return heating of last zone
 *
 ************************************************************************/

/*cdGetHeating_last routine to query results and return heating of last zone */

double cdGetHeating_last(void) /* return heating for last zone */
{
	return (thermal.htot);
}


/*************************************************************************
 *
 * cdNoExec call this routine to tell code not to actually execute
 *
 ************************************************************************/

/*cdNoExec call this routine to tell code not to actually execute */
#include "noexec.h"

void cdNoExec(void)
{

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

	/* this is an option to have the code read in all input quanties
	 * ot to NOT execute the actual model - to check on input parameters
	 * */
	noexec.lgNoExec = TRUE;


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


/*************************************************************************
 *
 * cdSetExecTime routine to initialize variables keeping track of time at start of calculation
 *
 ************************************************************************/

/* will be used to save initial time */
static clock_t before;
/* this must be true when cdExeTime is called, and is set true
 * when SetExecTime is called */
static int lgCalled=FALSE;

/* this routine is called by cdinit when everything is initialized,
 * so that every time cdExecTime is called the elapsed time is returned */
void cdSetExecTime(void)
{
	/* save startup time */
	before = clock();
	lgCalled = TRUE;
}

/*************************************************************************
 *
 * cdExecTime any routine can call this, and find the time since cdInit was called
 *
 ************************************************************************/

double cdExecTime(void)
{

	/* check that we were properly initialized */
	if( lgCalled )
	{
		/* this is elapsed time in seconds */
		return(  (double)(clock() - before) / (double)CLOCKS_PER_SEC )  ;
	}
	else
	{
		/* this is a big problem, we were asked for the elapsed time but
		 * the timer was not initialized by calling SetExecTime */
		fprintf( ioQQQ, " cdExecTime was called before SetExecTime, impossible.\n" );
		fprintf( ioQQQ, "Sorry.\n" );
		puts( "[Stop in cdExecTime]" );
		cdEXIT(EXIT_FAILURE);
	}
}

/*************************************************************************
 *
 * cdPrintCommands prints all input commands into file
 *
 ************************************************************************/

/* cdPrintCommands( FILE *)
 * prints all input commands into file */
void cdPrintCommands( FILE * ioOUT )
{
	long int i;
	fprintf( ioOUT, " Input commands follow:\n" );
	fprintf( ioOUT, "c ======================\n" );

	for( i=0; i <= input.nSave; i++ )
	{
		fprintf( ioOUT, "%s\n", input.chCardSav[i] );
	}
	fprintf( ioOUT, "c ======================\n" );
}


/*************************************************************************
 *
 * cdEms obtain the local emissivity for a line, for the last computed zone
 *
 ************************************************************************/

long int cdEmis(
	/* return value will be index of line within stack,
	 * negative of number of lines in the stack if the line could not be found*/
	/* 4 char null terminated string label */
	char *chLabel,
	/* line wavelength */
	float wavelength, 
	/* the vol emissivity of this line in last computed zone */
	double *emiss )
{
	/* use following to store local image of character strings */
	char chCARD[INPUT_LINE_LENGTH];
	char chCaps[5];
	long int j;
	float errorwave;

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

	/* routine returns the emissivity in the desired line
	 * only used internally by the code, to do punch lines structure */

	strcpy( chCARD, chLabel );

	/* make sure chLabel is all caps */
	caps(chCARD);/* convert to caps */

	/* get the error assocated with 4 significant figures */
	errorwave = WavlenErrorGet( wavelength );

	for( j=0; j < LineSave.nsum; j++ )
	{
		/* change chLabel to all caps to be like input chAmLab */
		cap4(chCaps , (char*)LineSv[j].chALab);

		/* check wavelength and chLabel for a match */
		/*if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength)<errorwave && 
			strcmp(chCaps,chCARD) == 0 ) */
		if( fabs(LineSv[j].wavelength-wavelength) < errorwave && strcmp(chCaps,chCARD) == 0 )
		{
			/* match, so set emiss to emissivity in line */
			*emiss = LineSv[j].emslin;
			
#			ifdef DEBUG_FUN
			fputs( " <->cdEms()\n", debug_fp );
#			endif
			/* and announce success by returning line index within stack */
			return j ;
		}
	}

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

	/* we fell through without finding the line - return false */
	return -LineSave.nsum;
}


/*************************************************************************
 *
 * cdEms_ip obtain the local emissivity for a line with known index
 *
 ************************************************************************/

void cdEmis_ip(
	/* index of the line in the stack */
	long int ipLine, 
	/* the vol emissivity of this line in last computed zone */
	double *emiss )
{

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

	/* routine returns the emissivity in the desired line
	 * used to do punch lines structure */

	/* match, so set emiss to emissivity in line */
	ASSERT( ipLine >= 0 && ipLine < LineSave.nsum );
	*emiss = LineSv[ipLine].emslin;

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

	return;
}

/*************************************************************************
 *
 * cdColm get the column density for a constituent - 0 return if ok, 1 if problems 
 *
 ************************************************************************/

int cdColm(
	/* return value is zero if all ok, 1 if errors happened */
	/* 4-char + eol string that is first 
	 * 4 char of element name as spelled by Cloudy, upper or lower case */
	char *chLabel,	

	/* integer stage of ionization, 1 for atom, 2 for A+, etc, 
	 * 0 is special flag for CO, H2, OH, or excited state */
	long int ion,

	/* the theoretical column density derived by the code */
	double *theocl )	
{
	long int nelem;
	/* use following to store local image of character strings */
	char chLABEL_CAPS[20];

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

	/* check that chLabel[4] is null - supposed to be 4 char + end */
	if( chLabel[4]!=0 )
	{
		fprintf( ioQQQ, " cdColm called with insane ion label, =%s, must be 4 character + end of string.\n", 
		  chLabel );
#		ifdef DEBUG_FUN
		fputs( " <->cdColm()\n", debug_fp );
#		endif
		return(1);
	}

	strcpy( chLABEL_CAPS, chLabel );
	/* convert element label to all caps */
	caps(chLABEL_CAPS);

	/* zero ionization stage has special meaning.  The quanties recognized are
	 * the molecules, "H2  ", "OH  ", "CO  ", etc 
	 * "CII*" excited state C+ */
	if( ion < 0 )
	{
		fprintf( ioQQQ, " cdColm called with insane ion, =%li\n", 
		  ion );
#		ifdef DEBUG_FUN
		fputs( " <->cdColm()\n", debug_fp );
#		endif
		return(1);
	}
	else if( ( strcmp( chLABEL_CAPS , "HYDR" )==0 ) && ion==3 )
	{
		/* special case, H  3 is H2 */
		*theocl = colden.colden[ipCOLH2];
	}
	else if( ion == 0 )
	{
		/* this case molecular column density */
		/* want the molecular hydrogen column density */
		if( strcmp( chLABEL_CAPS , "H2  " )==0 )
		{
			*theocl = colden.colden[ipCOLH2];
		}

		/* carbon monoxide column density */
		else if( strcmp( chLABEL_CAPS , "CO  " )==0 )
		{
			*theocl = co.hevcol[ipCO];
		}

		/* OH column density */
		else if( strcmp( chLABEL_CAPS , "OH  " )==0 )
		{
			*theocl = co.hevcol[ipOH];
		}

		/* H2O column density */
		else if( strcmp( chLABEL_CAPS , "H2O " )==0 )
		{
			*theocl = co.hevcol[ipH2O];
		}

		/* O2 column density */
		else if( strcmp( chLABEL_CAPS , "O2  " )==0 )
		{
			*theocl = co.hevcol[ipO2];
		}

		/* SiO column density */
		else if( strcmp( chLABEL_CAPS , "SIO " )==0 )
		{
			*theocl = co.hevcol[ipSIO];
		}

		/* CII^* column density, population of J=3/2 upper level of split ground term */
		else if( strcmp( chLABEL_CAPS , "CII*" )==0 )
		{
			*theocl = colden.C2Colden[1];
		}
		else if( strcmp( chLABEL_CAPS , "C11*" )==0 )
		{
			*theocl = colden.C1Colden[0];
		}
		else if( strcmp( chLABEL_CAPS , "C12*" )==0 )
		{
			*theocl = colden.C1Colden[1];
		}
		else if( strcmp( chLABEL_CAPS , "C13*" )==0 )
		{
			*theocl = colden.C1Colden[2];
		}
		else if( strcmp( chLABEL_CAPS , "O11*" )==0 )
		{
			*theocl = colden.O1Colden[0];
		}
		else if( strcmp( chLABEL_CAPS , "O12*" )==0 )
		{
			*theocl = colden.O1Colden[1];
		}
		else if( strcmp( chLABEL_CAPS , "O13*" )==0 )
		{
			*theocl = colden.O1Colden[2];
		}
		/* CIII excited states, upper level of 1909 */
		else if( strcmp( chLABEL_CAPS , "C30*" )==0 )
		{
			*theocl = colden.C3Colden[1];
		}
		else if( strcmp( chLABEL_CAPS , "C31*" )==0 )
		{
			*theocl = colden.C3Colden[2];
		}
		else if( strcmp( chLABEL_CAPS , "C32*" )==0 )
		{
			*theocl = colden.C3Colden[3];
		}
		else if( strcmp( chLABEL_CAPS , "SI2*" )==0 )
		{
			*theocl = colden.Si2Colden[1];
		}
		else if( strcmp( chLABEL_CAPS , "HE1*" )==0 )
		{
			*theocl = colden.He123S;
		}

		/* clueless as to what was meant - bail */
		else
		{
			fprintf( ioQQQ, " cdColm called with unknown element chLabel, org=\"%4.4s\" caps=\"%4.4s\"\n", 
			  chLabel , chLABEL_CAPS );
#			ifdef DEBUG_FUN
			fputs( " <->cdColm()\n", debug_fp );
#			endif
			return(1);
		}
	}
	else
	{
		/* this case, ionization stage of some element */
		/* find which element this is */
		nelem = 0;
		while( nelem < LIMELM && 
			strncmp(chLABEL_CAPS,elementnames.chElementNameShort[nelem],4) != 0 )
		{
			++nelem;
		}

		/* this is true if we have one of the first 30 elements in the label,
		 * nelem is on C scale */
		if( nelem < LIMELM )
		{

			/* sanity check - does this ionization stage exist?
			 * max2 is to pick up H2 as H 3 */
			if( ion > MAX2(3,nelem + 2) )
			{
				fprintf( ioQQQ, 
				  " cdColm asked to return ionization stage %ld for element %s but this is not physical.\n", 
				  ion, chLabel );
				return(1);
			}

			/* the column density, ion is on physics scale, but means are on C scale */
			*theocl = mean.xIonMeans[0][nelem][ion-1];
		}
		else
		{
			fprintf( ioQQQ, 
			  " cdColm did not understand this combinatino of ion %4ld and element %4.4s.\n", 
			  ion, chLabel );
			return(1);
		}
	}

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


/*************************************************************************
 *
 *cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit 
 *
 ************************************************************************/

void cdErrors(FILE *ioOUT)
{
	long int nc, 
	  nn, 
	  npe, 
	  ns, 
	  nte, 
	  nw ,
	  nIone,
	  nEdene;
	int lgAbort;

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

	/* first check for number of warnings, cautions, etc */
	cdNwcns(&lgAbort,&nw,&nc,&nn,&ns,&nte,&npe, &nIone, &nEdene );

	/* only say something is one of these problems is nonzero */
	if( nw || nc || nte || npe ||	nIone || nEdene || lgAbort )
	{
		/* say the title of the model */
		fprintf( ioOUT, "%75.75s\n", input.chTitle );

		if( lgAbort )
			fprintf(ioOUT," Calculation ended with abort!\n");

		/* print warnings on the io unit */
		if( nw != 0 )
		{
			cdWarnings(ioOUT);
		}

		/* print cautions on the io unit */
		if( nc != 0 )
		{
			cdCautions(ioOUT);
		}

		if( nte != 0 )
		{
			fprintf( ioOUT , "Te failures=%4ld\n", nte );
		}

		if( npe != 0 )
		{
			fprintf( ioOUT , "Pressure failures=%4ld\n", npe );
		}

		if( nIone != 0 )
		{
			fprintf( ioOUT , "Ionization failures=%4ld\n", nte );
		}

		if( nEdene != 0 )
		{
			fprintf( ioOUT , "Electron density failures=%4ld\n", npe );
		}
	}

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

/*************************************************************************
 *
 *cdGetDepth_depth get depth structure from previous iteration 
 *
 ************************************************************************/
void cdGetDepth_depth( double cdDepth[] )
{
	long int nz;

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

	for( nz = 0; nz<nzone; ++nz )
	{
		cdDepth[nz] = struc.depth[nz];
	}

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

/*************************************************************************
 *
 *cdGetPressure_depth routine to query results and return pressure of last iteration 
 *
 ************************************************************************/

/*
 * cdGetPressure_depth
 * This returns the pressure and its constituents for the last iteration. 
 * space was allocated in the calling routine for the vectors - 
 * before calling this, cdGetnZone should have been called to get the number of
 * zones, then space allocated for the arrays */
void cdGetPressure_depth(
	/* total pressure, all forms*/
	double TotalPressure[],			
	/* gas pressure */
	double GasPressure[],				
	/* radiation pressure */
	double RadiationPressure[])
{
	long int nz;

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

	for( nz = 0; nz<nzone; ++nz )
	{
		TotalPressure[nz] = struc.pressure[nz];
		GasPressure[nz] = struc.GasPressure[nz];
		RadiationPressure[nz] = struc.PresRadCurr[nz];
	}

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

/*************************************************************************
 *
 *cdGetPressure_last routine to query results and return pressure of last zone 
 *
 ************************************************************************/

void cdGetPressure_last(
		double *PresTotal,  /* total pressure, all forms, for the last computed zone*/
		double *PresGas,    /* gas pressure */
		double *PresRad)    /* radiation pressure */
{

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

	*PresGas = pressure.PresGasCurr;
	*PresRad = pressure.PresRadCurr;
	*PresTotal = pressure.PresGasCurr + pressure.PresRadCurr;

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

/*************************************************************************
 *
 *cdGetnZone gets number of zones
 *
 ************************************************************************/

/* returns number of zones */
long int cdGetnZone( void ) 
{
	return nzone;
}

/*************************************************************************
 *
 *cdGetTemp_last routine to query results and return temperature of last zone 
 *
 ************************************************************************/


double cdGetTemp_last(void)
{

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


#	ifdef DEBUG_FUN
	fputs( " <->cdGetTemp_last()\n", debug_fp );
#	endif
	return( phycon.te);
}

/*************************************************************************
 *
 *cdIonFrac get ionization fractions for a constituent
 *
 ************************************************************************/

int cdIonFrac(
	/* four char string, null terminzed, giving the element name */
	char *chLabel, 
	/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
	 * 0 says special case */
	long int IonStage, 
	/* will be fractional ionization */
	double *fracin, 
	/* how to weight the average, must be "VOLUME" or "RADIUS" */
	char *chWeight ,
	/* if true then weighting also has electron density, if false then only volume or radius */
	int lgDensity ) 
	/* return value is 0 if element was found, non-zero if failed */
{

	int lgVol;
	long int ip, 
		ion, /* used as index within aaa vector*/
		nelem;
	float aaa[LIMELM + 1];
	/* use following to store local image of character strings */
	char chCARD[INPUT_LINE_LENGTH];

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

	strcpy( chCARD, chWeight );
	/* make sure chWeight is all caps */
	caps(chCARD);/* convert to caps */

	/*caps(chWeight);*/

	if( strcmp(chCARD,"RADIUS") == 0 )
	{
		lgVol = FALSE;
	}

	else if( strcmp(chCARD,"VOLUME") == 0 )
	{
		lgVol = TRUE;
	}

	else
	{
		fprintf( ioQQQ, " cdIonFrac: chWeight=%6.6s makes no sense to me, valid options are RADIUS and VOLUME\n", 
		  chWeight );
		*fracin = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(1);
	}

	/* first ensure that chLabel is all caps */
	strcpy( chCARD, chLabel );
	/* make sure chLabel is all caps */
	caps(chCARD);/* convert to caps */

	if( IonStage==0 )
	{
		/* special case */
		if( strcmp(chCARD,"H2  " ) == 0 )
		{
			/* this will be request for H2, on c scale, hydro is 0 */
			nelem = 0;
			IonStage = 3;
		}
		else
		{
			fprintf( ioQQQ, " cdIonFrac: ion stage of zero and element %s makes no sense to me\n", 
			chCARD );
			*fracin = 0.;
			
#			ifdef DEBUG_FUN
			fputs( " <->cdIonFrac()\n", debug_fp );
#			endif
			return(1);
		}
	}

	else 
	{
		/* find which element this is, nelem is on physical scale, H is 1 */
		nelem = 0;
		while( nelem < LIMELM && 
			strcmp(chCARD,elementnames.chElementNameShort[nelem]) != 0 )
		{
			++nelem;
		}
		/* after this loop nelem is on c scale, H is 0 */
	}

	/* if element not recognized and above loop falls through, nelem is LIMELM 
	 * nelem counter is on physical scale, H = 1 Zn = 30 */
	if( nelem >= LIMELM )
	{
		fprintf( ioQQQ, " cdIonFrac called with unknown element chLabel, =%4.4s\n", 
		  chLabel );
		
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(1);
	}

	/* sanity check - does this ionization stage exist? 
	 * both IonStage and nelem are on physical scales, H atoms are 1 1 */

	/* IonStage came in on physics scale,
	 * ion will be used as pointer within the aaa array that contains average values,
	 * convert to C scale */
	ion = IonStage - 1;

	/*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/
	if( (ion > nelem+1 || ion < 0 ) && !(nelem==ipHYDROGEN&&ion==2))
	{
		fprintf( ioQQQ, " cdIonFrac asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
		  IonStage, chLabel );
		*fracin = -1.;
#		ifdef DEBUG_FUN
		fputs( " <->cdIonFrac()\n", debug_fp );
#		endif
		return(1);
	}

	/* get either volume or radius average, aaa is filled in from 0 */
	if( lgVol )
	{
		/* aaa is dim'd LIMELM+1 so largest arg is LIMELM
		 * 'i' means ionization, not temperature 
		 * nelem is on the physical scale */
		/* last arg says whether to include electron density */
		/* volmean uses c scale for nelem */
		MeanIonVolume('i', nelem,&ip,aaa,lgDensity);
		*fracin = pow(10.f,aaa[ion]);
	}
	else
	{
		/* last arg says whether to include electron density */
		/* volmean uses c scale for nelem */
		MeanIonRadius('i', nelem,&ip,aaa,lgDensity);
		*fracin = pow(10.f,aaa[ion]);
	}

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

	/* we succeeded - say so */
	return(0);
}

/*************************************************************************
 *
 * debugLine provides a debugging hook into the main line array 
 *
 ************************************************************************/

 /* following routine provides a debugging hook into the main line array 
 * loops over whole array and finds every line that matches length,
 * the wavelength, the argument to the function
 * put breakpoint inside if test */
long debugLine( float wavelength )
{
	long j, kount;
	float errorwave;

	kount = 0;

	/* get the error assocated with 4 significant figures */
	errorwave = WavlenErrorGet( wavelength );

	for( j=0; j < LineSave.nsum; j++ )
	{
		/* check wavelength and chLabel for a match */
		/* if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave ) */
		if( fabs(LineSv[j].wavelength-wavelength) < errorwave )
		{
			printf("%s\n", LineSv[j].chALab);
			++kount;
		}
	}
	printf(" hits = %li\n", kount );
	return(kount);
}

/*************************************************************************
 *
 *cdLineListPunch create a file with a list of all emission lines punched, 
 *and their index within the emission line stack
 *
 ************************************************************************/

/* returns total number of lines in the list */
long int cdLineListPunch( 
	/* a file handle pointing to a file that is read for writing -
	 * the calling routine must close it */
	 FILE* io )
{

	long int j;

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

	for( j=1; j < LineSave.nsum; j++ )
	{

		fprintf(io, "%li\t%s\t",
			j,
			LineSv[j].chALab);
		prt_wl( io , LineSv[j].wavelength );
		fprintf(io, "\n");
	}


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

	/* return total number of lines as debugging aid */
	return LineSave.nsum;
}

/*************************************************************************
 *
 *cdLine get the predicted line intensity, also index for line in stack 
 *
 ************************************************************************/

 /* returns array index for line in array stack if we found the line, 
  * or FALSE==0 if we did not find the line */
long int cdLine(char *chLabel, 
		/* wavelength of line in angstroms, not format printed by code */
	  float wavelength, 
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{
	char chCaps[5], 
	  chFind[5];
	long int ipobs, 
	  j;
	float errorwave;

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

	/* this is zero when cdLine called with cdNoExec called too */
	if( LineSave.nsum == 0 )
	{
		*relint =  0.;
		*absint = 0.;
		return 0;
	}
	ASSERT( LineSave.ipNormWavL >= 0);
	ASSERT( LineSave.nsum > 0);

	/* check that chLabel[4] is null - supposed to be 4 char + end */
	if( chLabel[4]!=0 )
	{
		fprintf( ioQQQ, " cdLine called with insane line label, =%s, must be 4 character + end of string.\n", 
		  chLabel );
#		ifdef DEBUG_FUN
		fputs( " <->cdColm()\n", debug_fp );
#		endif
		return(1);
	}

	/* change chLabel to all caps */
	cap4(chFind,chLabel);

	/* get the error assocated with 4 significant figures */
	errorwave = WavlenErrorGet( wavelength );

	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && fabs(wavelength-1000.)<0.01 )
		{
			fprintf(ioQQQ,"cddrive wl %.4e error %.3e\n",
				wavelength, errorwave );
		}
	}

	/* now go through entire line stack, do not do 0, which is H-beta and
	 * in stack further down - this is to stop query for H-beta from returning
	 * 0, the flag for line not found */
	for( j=1; j < LineSave.nsum; j++ )
	{

		/* check wavelength and chLabel for a match */
		/* DELTA since often want wavelength of zero */
		/*if( fabs(LineSv[j].wavelength-wavelength)/MAX2(DELTA,wavelength) < errorwave )*/
		if( fabs(LineSv[j].wavelength-wavelength) < errorwave )
		{

			/* change chLabel to all caps to be like input chALab */
			cap4(chCaps , (char*)LineSv[j].chALab);

			/* now see if labels agree */
			if( strcmp(chCaps,chFind) == 0 )
			{
				/* match, so set pointer */
				ipobs = j;

				/* does the normalization line have a positive intensity*/
				if( LineSv[LineSave.ipNormWavL].sumlin > 0. )
				{
					*relint = LineSv[ipobs].sumlin/LineSv[LineSave.ipNormWavL].sumlin*
					  LineSave.ScaleNormLine;
				}
				else
				{
					*relint = 0.;
				}

				/* return log of current line intensity if it is positive */
				if( LineSv[ipobs].sumlin > 0. )
				{
					*absint = log10(LineSv[ipobs].sumlin) + radius.Conv2PrtInten;
				}
				else
				{
					/* line intensity is actually zero, return small number */
					*absint = -37.;
				}
				
#				ifdef DEBUG_FUN
				fputs( " <->cdLine()\n", debug_fp );
#				endif
				/* we found the line, return pointer to its location */
				return ipobs;
			}
		}
	}

	*absint = 0.;
	*relint = 0.;

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

	/* if we fell down to here we did not find the line */
	/* >>chng 00 sep 02, return negative of total number of lines as debugging aid */
	return -LineSave.nsum;
}

/*************************************************************************
 *
 *cdLine_ip get the predicted line intensity, using index for line in stack 
 *
 ************************************************************************/

void cdLine_ip(long int ipLine, 
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{

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

	/* this is zero when cdLine called with cdNoExec called too */
	if( LineSave.nsum == 0 )
	{
		*relint =  0.;
		*absint = 0.;
		return;
	}
	ASSERT( LineSave.ipNormWavL >= 0);
	ASSERT( LineSave.nsum > 0);

	/* does the normalization line have a positive intensity*/
	if( LineSv[LineSave.ipNormWavL].sumlin > 0. )
	{
		*relint = LineSv[ipLine].sumlin/LineSv[LineSave.ipNormWavL].sumlin*
		  LineSave.ScaleNormLine;
	}
	else
	{
		*relint = 0.;
	}

	/* return log of current line intensity if it is positive */
	if( LineSv[ipLine].sumlin > 0. )
	{
		*absint = log10(LineSv[ipLine].sumlin) + radius.Conv2PrtInten;
	}
	else
	{
		/* line intensity is actually zero, return small number */
		*absint = -37.;
	}

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

	return;
}

/*************************************************************************
 *
 *cdDLine get the predicted emergent line intensity, also index for line in stack 
 *
 ************************************************************************/

 /* returns array index for line in array stack if we found the line, 
 * or FALSE==0 if we did not find the line */
long int cdDLine(char *chLabel, 
		/* wavelength of line as printed by code*/
	  float wavelength, 
	  /* linear intensity relative to normalization line*/
	  double *relint, 
	  /* log of luminosity or intensity of line */
	  double *absint )
{
	char chCaps[5], 
	  chFind[5];
	long int ipobs, 
	  j;
	float errorwave;

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

	/* this is zero when cdDLine called with cdNoExec called too */
	if( LineSave.ndsum == 0 )
	{
		*relint =  0.;
		*absint = 0.;
		return 0;
	}
	ASSERT( LineSave.ipNormWavL >= 0);
	ASSERT( LineSave.nsum > 0);

	/* change chLabel to all caps */
	cap4(chFind,chLabel);

	/* get the error assocated with 4 significant figures */
	errorwave = WavlenErrorGet( wavelength );

	/* now go through entire line stack, do not do 0, which is H-beta and
	 * in stack further down - this is to stop query for H-beta from returning
	 * 0, the flag for line not found */
	for( j=1; j < LineSave.ndsum; j++ )
	{

		/* check wavelength and chLabel for a match */
		/* DELTA since often want wavelength of zero */
		/* if( fabs(LineDSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave ) */
		if( fabs(LineSv[j].wavelength-wavelength) < errorwave )
		{

			/* change chLabel to all caps to be like input chALab */
			cap4(chCaps , (char*)LineDSv[j].chSMDLab);

			/* now see if labels agree */
			if( strcmp(chCaps,chFind) == 0 )
			{
				/* match, so set array index */
				ipobs = j;

				/* does the normalization line have a positive intensity*/
				if( LineDSv[LineSave.ipEmerNormWavL].smdlin > 0. )
				{
					*relint = LineDSv[ipobs].smdlin/LineDSv[LineSave.ipEmerNormWavL].smdlin*
					  LineSave.ScaleNormLine;
				}
				else
				{
					*relint = 0.;
				}

				/* return log of current line intensity if it is positive */
				if( LineDSv[ipobs].smdlin > 0. )
				{
					*absint = log10(LineDSv[ipobs].smdlin) + radius.Conv2PrtInten;
				}
				else
				{
					/* line intensity is actually zero, return small number */
					*absint = -37.;
				}
				
#				ifdef DEBUG_FUN
				fputs( " <->cdDLine()\n", debug_fp );
#				endif
				/* we found the line, return pointer to its location */
				return ipobs;
			}
		}
	}

	*absint = 0.;
	*relint = 0.;

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

	/* if we fell down to here we did not find the line */
	/* >>chng 00 sep 02, return negative of total number of lines as debugging aid */
	return -LineSave.nsum;
}

/*************************************************************************
 *
 *cdNwcns get the number of cautions and warnings, to tell if calculation is ok
 *
 ************************************************************************/

void cdNwcns(
  /* abort status, this better be false */
  int *lgAbort ,
  /* the number of warnings, cautions, notes, and surprises */
  long int *NumberWarnings, 
  long int *NumberCautions, 
  long int *NumberNotes, 
  long int *NumberSurprises, 
  /* the number of temperature convergence failures */
  long int *NumberTempFailures, 
  /* the number of pressure convergence failures */
  long int *NumberPresFailures,
  /* the number of ionzation convergence failures */
  long int *NumberIonFailures, 
  /* the number of electron density convergence failures */
  long int *NumberNeFailures )
{

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

	/* this would be set true if code ended with abort - very very bad */
	*lgAbort = conv.lgAbort;
	/* this sub is called after comment, to find the number of various comments */
	*NumberWarnings = warnings.nwarn;
	*NumberCautions = warnings.ncaun;
	*NumberNotes = warnings.nnote;
	*NumberSurprises = warnings.nbang;

	/* these are counters that were incremented during convergence failures */
	*NumberTempFailures = conv.nTeFail;
	*NumberPresFailures = conv.nPreFail;
	*NumberIonFailures = conv.nIonFail;
	*NumberNeFailures = conv.nNeFail;

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

/*************************************************************************
 *
 * cdOutp redirect output to arbitrary openned file
 *
 ************************************************************************/

void cdOutp( FILE *ioOut)
{

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

	/* ioQQQ is pointer to output fiile */
	ioQQQ = ioOut;

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

/*************************************************************************
 *
 * cdInp redirect line input from arbitrary openned file
 *
 ************************************************************************/

void cdInp( FILE *ioInp)
{

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

	/* ioQQQ is pointer to output fiile */
	ioStdin = ioInp;

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

/*************************************************************************
 *
 *cdTalk tells the code whether to print results or be silent
 *
 ************************************************************************/

void cdTalk(int lgTOn)
{

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

	called.lgTalk = lgTOn;

	if( lgTOn )
	{
		/* means talk has not been forced off */
		called.lgTalkForcedOff = FALSE;
	}
	else
	{
		/* means talk is forced off */
		called.lgTalkForcedOff = TRUE;
	}

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

/*************************************************************************
 *
 * cdTemp get mean electron temperature for any element 
 *
 ************************************************************************/

/* This routine finds the mean electron temperature for any ionization stage 
 * It returns 0 if it could find the species, 1 if it could not find the species.
 * The first argument is a null terminated 4 char string that gives the element
 * name as spelled by Cloudy.  
 * The second argument is ion stage, 1 for atom, 2 for first ion, etc
 * This third argument will be returned as result,
 * Last parameter is either "VOLUME" or "RADIUS" to give weighting 
 *
 * if the ion stage is zero then the element label will have a special meaning.
 * The string "21CM" is will return the 21 cm temperature.
 * The string "H2  " will return the temperature weighted wrt the H2 density */
/*TODO	2	this should have a last argument like cdIonFrac for whether or not weighting
 * is wrt electron density */

/* return value is o if things are ok and element was found, 
 * non-zero if element not found or there are problems */
int cdTemp(
	/* four char string, null terminzed, giving the element name */
	char *chLabel, 
	/* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
	 * 0 means that chLabel is a special case */
	long int IonStage, 
	/* will be temperature */
	double *TeMean, 
	/* how to weight the average, must be "VOLUME" or "RADIUS" */
	char *chWeight ) 
{

	int lgVol;
	long int ip, 
		ion, /* used as pointer within aaa vector*/
		nelem;
	float aaa[LIMELM + 1];
	/* use following to store local image of character strings */
	char chWGHT[INPUT_LINE_LENGTH] , chELEM[INPUT_LINE_LENGTH];

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

	/* make sure chWeight is all caps */
	strcpy( chWGHT, chWeight );
	caps(chWGHT);/* convert to caps */

	/* ensure that chLabel is all caps */
	strcpy( chELEM, chLabel );
	caps(chELEM);/* convert to caps */

	/* now see if volume or radius weighting */
	if( strcmp(chWGHT,"RADIUS") == 0 )
	{
		lgVol = FALSE;
	}
	else if( strcmp(chWGHT,"VOLUME") == 0 )
	{
		lgVol = TRUE;
	}
	else
	{
		fprintf( ioQQQ, " cdTemp: chWeight=%6.6s makes no sense to me, the options are RADIUS and VOLUME.\n", 
		  chWeight );
		*TeMean = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(1);
	}

	if( IonStage == 0 )
	{
		/* return atomic hydrogen weighted harmonic mean of gas kinetic temperature */
		if( strcmp(chELEM,"21CM") == 0 )
		{
			if( lgVol )
			{
				/* this is mean of inverse temperature weighted over volume */
				if( mean.HarMeanTempVolume[1] > SMALLFLOAT )
				{
					*TeMean = mean.HarMeanTempVolume[0]/mean.HarMeanTempVolume[1];
				}
				else
				{
					*TeMean = 0.;
				}
			}
			else
			{
				/* this is mean of inverse temperature weighted over radius */
				if( mean.HarMeanTempRadius[1] > SMALLFLOAT )
				{
					*TeMean = mean.HarMeanTempRadius[0]/mean.HarMeanTempRadius[1];
				}
				else
				{
					*TeMean = 0.;
				}
			}
		}
		/* return atomic hydrogen weighted harmonic mean 21 cm spin temperature,
		 * using actual level populations with 1s of H0 */
		else if( strcmp(chELEM,"SPIN") == 0 )
		{
			*TeMean = mean.H_21cm_spin_mean_radius[0] / 
				SDIV(mean.H_21cm_spin_mean_radius[1]);
		}
		/* return temperature deduced from ratio of 21 cm and Lya optical depths */
		else if( strcmp(chELEM,"OPTI") == 0 )
		{
			/* this is the temperature derived from Lya - 21 cm optical depths */
			*TeMean = 
				3.84e-7 * EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon /
				SDIV( HFLines[0].TauCon );
		}
		/* special case, mean wrt H_2 */
		else if( strcmp(chELEM,"H2  ") == 0 )
		{
			if( lgVol )
			{
				/* mean temp with H2 over volume */
				if( mean.H2MeanTempVolume[1] > SMALLFLOAT )
				{
					*TeMean = mean.H2MeanTempVolume[0] / mean.H2MeanTempVolume[1];
				}
				else
				{
					*TeMean = 0.;
				}
			}
			else
			{
				/* mean temp with H2 over radius */
				if( mean.H2MeanTempRadius[1] > SMALLFLOAT )
				{
					*TeMean = mean.H2MeanTempRadius[0] / mean.H2MeanTempRadius[1];
				}
				else
				{
					*TeMean = 0.;
				}
			}
		}
		else
		{
			fprintf( ioQQQ, " cdTemp called with ion=0 and unknown quantity, =%4.4s\n", 
			  chLabel );
			fprintf( ioQQQ, " I know about 21CM and H2__ \n");
			
#			ifdef DEBUG_FUN
			fputs( " <->cdTemp()\n", debug_fp );
#			endif
			return(1);
		}

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

		/* say things are ok */
		return(0);
	}

	/* find which element this is */
	nelem = 0;
	while( nelem < LIMELM && 
		strcmp(chELEM,elementnames.chElementNameShort[nelem]) != 0 )
	{
		++nelem;
	}
	/* after this loop nelem is atomic number of element, H is 1 */

	/* if element not recognized and above loop falls through, nelem is LIMELM+1 
	 * nelem counter is on physical scale, H = 1 Zn = 30 */
	if( nelem >= LIMELM )
	{
		fprintf( ioQQQ, " cdTemp called with unknown element chLabel, =%4.4s\n", 
		  chLabel );
		
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(1);
	}

	/* sanity check - does this ionization stage exist? 
	 * both IonStage on physical scale, nelem on c */

	/* ion will be used as pointer within the aaa array that contains average values,
	 * done this way to prevent lint from falsing in acess to aaa array */
	ion = IonStage - 1;

	/*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/
	if( ion > nelem+1 || ion < 0 || ion > LIMELM )
	{
		fprintf( ioQQQ, " cdTemp asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 
		  IonStage, chLabel );
#		ifdef DEBUG_FUN
		fputs( " <->cdTemp()\n", debug_fp );
#		endif
		return(1);
	}

	/* get either volume or radius average, aaa is filled in from 0 */
	if( lgVol )
	{
		/* aaa is dim'd LIMELM+1 so largest arg is LIMELM */
		/* MeanIonVolume uses C scale for nelem */
		MeanIonVolume('t', nelem,&ip,aaa,FALSE);
		*TeMean = pow(10.f,aaa[ion]);
	}
	else
	{
		MeanIonRadius('t', nelem,&ip,aaa,FALSE);
		*TeMean = pow(10.f,aaa[ion]);
	}

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

/*************************************************************************
 *
 * cdRead routine to read in command lines 
 * called by user when cloudy used as subroutine 
 * called by maincl when used as a routine
 *
 ************************************************************************/

/* returns the number of lines available in command stack 
 * this is limit to how many more commands can be read */
int cdRead(
	/* the string containing the commands */
	char *chInputLine )	
{
	char *chEOL , /* will be used to seach for end of line symbols */
		chCARD[INPUT_LINE_LENGTH],
		chLocal[INPUT_LINE_LENGTH];

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

	/* this is set FALSE when code loaded, set TRUE when cdInit called,
	 * this is check that cdInit was called first */
	if( !lgcdInitCalled )
	{
		printf(" cdInit was not called first - this must be the first call.\n");
		puts( "[Stop in cdRead]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* totally ignore any card starting with a #, *, space, //, or % */
	if( ( chInputLine[0]=='#') || 
		 (chInputLine[0]=='%') ||
		 (strncmp(chInputLine,"//", 2 )==0 ) || 
		 (chInputLine[0]=='*') || 
		 /* these two are end-of-input-stream sentinals */
		 (chInputLine[0]=='\n')||
		 (chInputLine[0]==' ') )
	{
		/* return value is number of lines that can still be stuffed in */
		return(NKRD - input.nSave);
	}

	/***************************************************************************
	* validate a location to store this line image, then store the version     *
	* that has been truncated from special end of line characters              *
	* stored image will have <=80 valid characters                             *
	***************************************************************************/

	/* this will now point to the next free slot in the line image save array 
	 * this is where we will stuff this line image */
	++input.nSave;

	if( input.nSave >= NKRD )
	{
		/* too many input commands were entered - bail */
		fprintf( ioQQQ, 
			" Too many line images entered to cdRead.  The limit is %d\n", 
		  NKRD );
		fprintf( ioQQQ, 
			" The limit to the number of allowed input lines is %d.  This limit was exceeded.  Sorry.\n", 
		  NKRD );
		fprintf( ioQQQ, 
			" This limit is set by the variable NKRD which appears in input.h \n" );
		fprintf( ioQQQ, 
			" Increase it everywhere it appears.\n" );
		puts( "[Stop in cdRead]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now kill any part of line image after special end of line character,
	 * this stops info user wants ignored from entering after here */
	if( (chEOL = strchr(chInputLine , '\n' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	if( (chEOL = strchr(chInputLine , '%' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	/* >>chng 02 apr 10, add this char */
	if( (chEOL = strchr(chInputLine , '#' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	if( (chEOL = strchr(chInputLine , ';' ) ) !=NULL )
	{
		*chEOL = '\0';
	}
	if( (chEOL = strstr(chInputLine , "//" ) ) !=NULL )
	{
		*chEOL = '\0';
	}

	/* check that there are less than INPUT_LINE_LENGTH characters to the null character 
	 * first find the end of string mark */
	chEOL = strchr(chInputLine , '\0' );

	/* do something if no end of string, or if string longer than 
	 * INPUT_LINE_LENGTH characters */
	if( (chEOL==NULL) || (chEOL - chInputLine)>INPUT_LINE_LENGTH )
	{
		/* remember, in C [INPUT_LINE_LENGTH] is the n+1th character */
		chInputLine[INPUT_LINE_LENGTH-1] = '\0' ;
		fprintf(ioQQQ," cdRead, while parsing the following input line:\n %s\n",
			chInputLine);
		fprintf(ioQQQ," found that the line is longer than %i characters, the longest possible line.\n",
			INPUT_LINE_LENGTH);
		fprintf(ioQQQ," Please make the command line shorter than this limit.\n");
	}

	/* now do it again, since we now want to make sure that there is a trailing space
	 * if the line is shorter than 80 char, test on null is to keep lint happy */
	if( (chEOL = strchr(chInputLine , '\0' )) == NULL )
		TotalInsanity();

	/* pad with a space if short enough,
	 * if not short enough for this to be done, then up to user to create correct input */
	strcpy( chLocal , chInputLine );
	if( chEOL-chInputLine <INPUT_LINE_LENGTH-4 )
	{
		strcat( chLocal , "  " );
	}

	/* save string in master array for later use in readr */
	strcpy( input.chCardSav[input.nSave], chLocal );

	/* copy string over the chCARD, then convert this to caps to check for optimize*/
	strcpy( chCARD, chLocal );
	caps(chCARD);/* convert to caps */

	/* check whether this is a trace command, turn on printout if so */
	if( strncmp(chCARD,"TRACE",5) == 0 )
	{
		trace.lgTraceInput = TRUE;
	}

	/* print input lines if trace specified */
	if( trace.lgTraceInput )
	{
		fprintf( ioQQQ,"cdRead=%s=\n",input.chCardSav[input.nSave] );
	}

	/* now check whether VARY is specified */
	if( lgMatch("VARY",chCARD) )
	{
		/* a optimize flag was on the line image */
		optimize.lgVaryOn = TRUE;
	}

	/* now check whether line is "no optimize" command */
	if( strncmp(chCARD,"NO VARY",7) == 0 )
	{
		optimize.lgNoVary = TRUE;
	}

	if( strncmp(chCARD,"OPTI",4) == 0 )
	{
		/* optimize command read in */
		optimize.lgOptimr = TRUE;
	}

#	ifdef DEBUG_FUN
	fputs( " <->cdRead()\n", debug_fp );
#	endif
	return( NKRD - input.nSave );
}

/* print line wavelengths in Angstroms in the standard format - just a wrapper */
void cdPrtWL( FILE *io , float wl )
{

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

	prt_wl( io , wl );

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

/* wrapper to close all punch files */
void cdClosePunchFiles( void )
{

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

	ClosePunchFiles();

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

#if 0
/*cdTemp21cmLya this is the temperature derived from Lya - 21 cm optical depths */
double cdTemp21cmLya( void )
{
	double t;

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

	/*fprintf(ioQQQ,"\n bug debug 21cm %.3e %.3e %.3e\n", 
		EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon ,
		 HFLines[0].TauCon , 
		 3.84e-7 * EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon /
		 SDIV( HFLines[0].TauCon ) );*/

	/* this is the temperature derived from Lya - 21 cm optical depths */
	t = 
		3.84e-7 * EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauCon /
		SDIV( HFLines[0].TauCon );

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