/* 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 */
/*PutLine enter local line intensity into the intensity stack for eventual printout */
/*PutExtra enter and 'extra' intensity source for some line */
/*linadd enter lines into the line storage array, called once per zone */
/*DumpLine print various information about an emission line vector, 
 * used in debugging, print to std out, ioQQQ */
/*TexcLine derive excitation temperature of line from contents of line array */
/*GetGF convert Einstein A into oscillator strength */
/*emit_frac returns fraction of populations the produce emission */
/*abscf convert gf into absorption coefficient */
/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
/*chLineLbl use information in line transfer arrays to generate a line label */
/*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
 * used to convert vacuum wavelengths to air wavelengths. */
/*eina convert a gf into an Einstein A */
/*WavlenErrorGet - find difference between two wavelengths */
/*PutCS enter a collision strength into an individual line vector */
/*lindst add local line intensity to line luminosity stack */
/*PntForLine generate pointer for forbidden line */
/*EmLineJunk set all elements of EmLine struc to dangerous values */
/*EmLineZero set all elements of EmLine struc to zero */
/*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
/*lgTauGood returns true is we have not overrun optical depth scale */
/*OccupationNumberLine - derive the photon occupation number at line center for any line */
/*outline - adds line photons to reflin and outlin */
/*MakeCS compute collision strength by g-bar approximations */
/*gbar1 compute g-bar collision strength using Mewe approximations */
/*gbar0 compute g-bar gaunt factor for neutrals */
/*totlin sum total intensity of cooling, recombination, or intensity lines */
/*FndLineHt search through line heat arrays to find the strongest heat source */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "lines.h"
#include "radius.h"
#include "elementnames.h"
#include "dense.h"
#include "rfield.h"
#include "opacity.h"
#include "grainvar.h"
#include "ipoint.h"
#include "iso.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "lines_service.h"

/*outline - adds line photons to reflin and outlin */
void outline( EmLine *t )
{
	long int ip = t->ipCont-1;
	double xInWrd = t->phots*t->FracInwd;

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

	/* the reflected part */
	rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn);

	/* inward beam that goes out since sphere set */
	rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]*
		t->ColOvTot);

	/* outward part */
	rfield.outlin[ip] += (float)(t->phots*
		(1. - t->FracInwd)*radius.BeamOutOut* t->ColOvTot);

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

	return;
}

/* define current transition probabilities and opacities for H lines */
void HLineTransOpacSet( long int nelem )
{
	long int limit , ipHi , ipLo ;
	double z4 , factor;

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

	/* will need this for some scaling laws - physical Z to the 4th power*/
	z4 = POW2(nelem+1.);
	z4 *= z4;

	/* hydrogen lines opacity is function of A's
		* above 15 just use actual A's for high density limit
		* only change opacities for levels between 4 and 15 */
	/* NB! must never set opacity for 2s-1s two-photon transition here, since
		* the "line" can become optically thick if treated as a line.  It is a 
		* continuum and the opacity is not related to the A by the usual expression */
	limit = MIN2(16,iso.numLevels[ipH_LIKE][nelem]);
	/* do Paschen and higher lines first, 
	 * do balmer below since must separate 2s and 2p */
	for( ipHi=4; ipHi < limit; ipHi++ )
	{
		for( ipLo=3; ipLo < ipHi; ipLo++ )
		{
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul = 
				(float)(hydro.HyLife[ipHi]*
				HydroBranch(ipHi,ipLo,nelem+1)*z4);
			ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul > 0.);

			/* make self-consistent opacity, convert new As back into opacities */
			EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity = 
				(float)(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
				2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
				iso.stat[ipH_LIKE][nelem][ipLo]*
				POW3(RYDLAM/(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN * WAVNRYD)));

			/* check that results are ok */
			ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity > 0.);
		}
	}

	/* the actual branching ratio from high levels down to
		* 2s and 2p depend on the density.  the code goes to
		* the correct low and high density limits - I know of
		* nothing better than can be done. */
	factor = MAX2( 0.25 , 0.32 - 0.07*dense.eden/(dense.eden + 1e7) );

	/* treat 2s = to 2p for HydroBranch, which returns total to 2s+2p */
	for( ipHi=4; ipHi < limit; ipHi++ )
	{
		/* get new effective A for this density and temperature */
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul = 
			(float)(hydro.HyLife[ipHi]*
			factor*HydroBranch(ipHi,2,nelem+1)*z4);

			/* check that results are ok */
		ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul > 0.);

		/* do 2p by scaling relative to 2s */
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul = (float)(
			EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul / factor *( 1. - factor )); 

			/* check that results are ok */
		ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul > 0.);

		/* make self-consistent opaciity for 2s, from A */
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].opacity = 
			(float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul*
			2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
			iso.stat[ipH_LIKE][nelem][ipH2s]*
			POW3(RYDLAM/(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].EnergyWN * WAVNRYD)));

		/* check that results are ok */
		ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].opacity > 0.);

		/* make self-consistent opaciity for 2p, from A */
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].opacity = 
			(float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul*
			2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
			iso.stat[ipH_LIKE][nelem][ipH2p]*
			POW3(RYDLAM/ (EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].EnergyWN * WAVNRYD)));

			/* check that results are ok */
		ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].opacity > 0.);
	}

	return;

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

}

/*emit_frac returns fraction of populations the produce emission */
double emit_frac( EmLine *t )
{
	/* collisional deexcitation and dest by background opacitities are loss of photons
	 * without net emission */
	double deexcit_loss = t->cs * phycon.cdsqte + t->Aul*t->Pdest;
	/* this is what is observed */
	double rad_deexcit = t->Aul*(t->Pelec_esc + t->Pesc);

	return rad_deexcit/(deexcit_loss + rad_deexcit);
}

/*DumpLine print various information about an emission line vector, 
 * used in debugging, print to std out, ioQQQ */
void DumpLine(EmLine * t)
{
	char chLbl[11];

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

	/* routine to print contents of line arrays */
	strcpy( chLbl, chLineLbl(t) );

	fprintf( ioQQQ, 
		"   %10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e col%.1e het%.1e conopc%.1e albdo%.2e\n", 
	  chLbl, 
	  phycon.te, 
	  dense.eden, 
	  t->cs, 
	  t->Aul, 
	  TexcLine(t), 
	  t->cool, 
	  t->heat ,
	  opac.opacity_abs[t->ipCont-1],
	  opac.albedo[t->ipCont-1]);

	fprintf( ioQQQ, 
		"  Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n", 
	  t->TauIn, 
	  t->TauTot, 
	  t->Pesc, 
	  t->Pelec_esc, 
	  t->Pdest, 
	  t->pump, 
	  t->ots, 
	  t->PopLo, 
	  t->PopHi ,
	  t->PopOpc );

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

	return;
}


/*OccupationNumberLine - derive the photon occupation number at line center for any line */
double OccupationNumberLine( EmLine * t )
{
	double OccupationNumberLine_v;

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

	/* routine to evaluate line photon occupation number */
	if( t->PopLo > SMALLFLOAT )
	{
		/* the lower population with correction for stimulated emission */
		double PopLo_corr = t->PopLo / t->gLo - t->PopHi / t->gHi;
		OccupationNumberLine_v = ( t->PopHi / t->gHi )/SDIV( PopLo_corr ) *
			(1. - t->Pesc);
	}
	else
	{
		OccupationNumberLine_v = 0.;
	}

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

	return( OccupationNumberLine_v );
}

/*TexcLine derive excitation temperature of line from contents of line array */
double TexcLine(EmLine * t)
{
	double TexcLine_v;

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

	/* routine to evaluate line excitation temp using contents of line array
	 * */
	if( t->PopHi * t->PopLo > 0. )
	{
		TexcLine_v = ( t->PopHi / t->gHi )/( t->PopLo / t->gLo );
		TexcLine_v = log(TexcLine_v);
		/* protect against infinite temp limit */
		if( fabs(TexcLine_v) > SMALLFLOAT )
		{
			TexcLine_v = - t->EnergyK / TexcLine_v;
		}
	}
	else
	{
		TexcLine_v = 0.;
	}

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

/*eina convert a gf into an Einstein A */
double eina(double gf, 
	  double enercm, 
	  double gup)
{
	double eina_v;

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

	/* derive the transition prob, given the following
	 * call to function is gf, ener in cm^-1, g_up
	 * gf is product of g and oscilator strength
	 * eina = ( gf / 1.499e-8 ) / (wl/1e4)**2 / gup  */
	eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);

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

	return( eina_v );
}

/*GetGF convert Einstein A into oscillator strength */
double GetGF(double trans_prob, 
	  double enercm, 
	  double gup)
{
	double GetGF_v;

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

	ASSERT( enercm > 0. );
	ASSERT( trans_prob > 0. );
	ASSERT( gup > 0.);

	/* derive the transition prob, given the following
	 * call to function is gf, ener in cm^-1, g_up
	 * gf is product of g and oscilator strength
	 * trans_prob = ( GetGF/gup) / 1.499e-8 / ( 1e4/enercm )**2 */
	GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
	
#	ifdef DEBUG_FUN
	fputs( " <->GetGF()\n", debug_fp );
#	endif

	return( GetGF_v );
}

/*abscf convert gf into absorption coefficient */
double abscf(double gf, 
	  double enercm, 
	  double gl)
{
	double abscf_v;

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

	ASSERT(gl > 0. && enercm > 0. && gf > 0. );

	/* derive line absorption coefficient, given the following:
	 * gf, enercm, g_low
	 * gf is product of g and oscillator strength */
	abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);

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

	return( abscf_v );
}

/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
void chIonLbl(char *chIonLbl_v , EmLine * t )
{
	/*static char chIonLbl_v[5];*/

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

	/* function to use information within the line array
	 * to generate an ion label, giving element and
	 * ionization stage
	 * */
	if( t->nelem < 0 )
	{
		/* this line is to be ignored */
		strcpy( chIonLbl_v, "Dumy" );
	}
	else if( t->nelem-1 >= LIMELM )
	{
		/* this is one of the molecules, either 12CO or 13CO */

		/* >>chng 02 may 15, from to chElementNameShort to go mole right */
		strcpy( chIonLbl_v , elementnames.chElementNameShort[t->nelem-1] );

		/* chIonStage is four char null terminated string, starting with "_1__" 
		strcat( chIonLbl_v, "CO");*/
	}

	else
	{
		/* ElmntSym.chElementSym is null terminated, 2 ch + null, string giving very
		 * short form of element name */
		strcpy( chIonLbl_v , elementnames.chElementSym[t->nelem -1] );

		/* chIonStage is four char null terminated string, starting with "_1__" */
		strcat( chIonLbl_v, elementnames.chIonStage[t->IonStg-1]);
	}

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

	/* chIonLbl is four char null terminated string */
	return/*( chIonLbl_v )*/;
}

/*chLineLbl use information in line transfer arrays to generate a line label */
/* this label is null terminated */
/* ContCreatePointers has test this with full range of wavelengths */
char* chLineLbl(EmLine * t )
{
	static char chLineLbl_v[11];

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


	/* function to use information within the line array
	 * to generate a line label, giving element and
	 * ionization stage
	 * rhs are set in large block data */

	/* NB this funciton is profoundly slow due to sprintf statement
	 * also - it cannot be evaluated within a write statement itself*/

	if( t->WLAng > (float)INT_MAX )
	{
		sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
			elementnames.chElementSym[t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
		   (int)(t->WLAng/1e8), 'c' );
	}
	else if( t->WLAng > 9999999. )
	{
		/* wavelength is very large, convert to centimeters */
		sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
			elementnames.chElementSym[t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
			t->WLAng/1e8, 'c' );
	}
	else if( t->WLAng > 999999. )
	{
		/* wavelength is very large, convert to microns */
		sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
			elementnames.chElementSym[t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
			(int)(t->WLAng/1e4), 'm' );
	}
	else if( t->WLAng > 99999. )
	{
		/* wavelength is very large, convert to microns */
		sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c", 
			elementnames.chElementSym[t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
			t->WLAng/1e4, 'm' );
	}
	else if( t->WLAng > 9999. )
	{
		sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
			elementnames.chElementSym[ t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
		   t->WLAng/1e4, 'm' );
	}
	else if( t->WLAng >= 100. )
	{
		sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
			elementnames.chElementSym[ t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
		   (int)t->WLAng, 'A' );
	}
	/* the following two formats should be changed for the
	 * wavelength to get more precision */
	else if( t->WLAng >= 10. )
	{
		sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c", 
			elementnames.chElementSym[ t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
		   t->WLAng, 'A' );
	}
	else
	{
		sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
			elementnames.chElementSym[ t->nelem -1], 
			elementnames.chIonStage[t->IonStg-1], 
		   t->WLAng, 'A' );
	}
	/* make sure that string ends with null character - this should
	 * be redundant */
	ASSERT( chLineLbl_v[10]==0 );

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

	return( chLineLbl_v );
}

/*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
 * used to convert vacuum wavelengths to air wavelengths. */
double RefIndex(EmLine * t)
{
	double RefIndex_v, 
	  WaveMic, 
	  xl, 
	  xn;

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

	/* the wavelength in microns */
	WaveMic = 1.e4/t->EnergyWN;

	/* only do index of refraction if longward of 2000A */
	if( WaveMic > 0.2 )
	{
		/* longward of 2000A
		 * xl is 1/WaveMic^2 */
		xl = 1.0/WaveMic/WaveMic;
		/* use a formula from Allen's 3rd Edition AQ, p.124 */
		xn = 255.4/(41. - xl);
		xn += 29498.1/(146. - xl);
		xn += 64.328;
		RefIndex_v = xn/1.e6 + 1;
	}
	else
	{
		RefIndex_v = 1.;
	}

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

	return( RefIndex_v );
}

/*PutCS enter a collision strength into an individual line vector */
void PutCS(double cs, 
  /*float tarray[])*/
  EmLine * t)
{

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

	/* collision strength must not be negative, had been test for being positive,
	 * but called with zero - did not check why 98 jul 5 */
	ASSERT( cs >= 0. );

	t->cs = (float)cs;

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

	return;
}

/*WavlenErrorGet - given the real wavelength in A for a line
 * routine will find the error expected between the real 
 * wavelength and the wavelength printed in the output, with 4 sig figs,
 * function returns difference between exact and 4 sig fig wl, so 
 * we have found correct line is fabs(d wl) < return */
float WavlenErrorGet( float wavelength )
{
	double a;
	float errorwave;

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

	if( wavelength > 0. )
	{
		/* normal case, positive (non zero) wavelength */
		a = log10( wavelength+FLT_EPSILON);
		a = floor(a);

		errorwave = 5.f * (float)pow( 10., a - 4. );
	}
	else
	{
		/* might be called with wl of zero, this is that case */
		errorwave = 1e-4f;
		a = 0.;
	}

	return errorwave;

}

/*linadd enter lines into the line storage array, called once per zone for each line*/
void linadd(
  double xInten,	/* xInten - local emissivity per unit vol, no fill fac */
  float wavelength,	/* float wavelength */
  char *chLab,		/* string label for ion */
  char chInfo 		/* character type of entry for line - given below */
					/* 'c' cooling, 'h' heating, 'i' info only, 'r' recom line */)
{

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

	/* main routine to actually enter lines into the line storage array
	 * called at top level within routine lines
	 * called series of times in routine PutLine for lines transferred
	 */

	/* three values, -1 is just counting, 0 if init, 1 for calculation */
	if( LineSave.ipass > 0 )
	{
		/* not first pass, sum lines only
		 * total emission from vol */
		LineSv[LineSave.nsum].sumlin += xInten*radius.dVeff;
		/* local emissivity in line */
		LineSv[LineSave.nsum].emslin = xInten;
	}

	else if( LineSave.ipass == 0 )
	{
		/* first call to suff lines in array, confirm that label is one of
		 * the four correct ones */
		/* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
		ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );

		/* then save it into array */
		LineSv[LineSave.nsum].chSumTyp = (char)chInfo;

		/* number of lines ok, set parameters for first pass */
		LineSv[LineSave.nsum].sumlin = 0.;
		LineSv[LineSave.nsum].wavelength = wavelength;
		strcpy( LineSv[LineSave.nsum].chALab, chLab );
	}

	/* increment the line counter */
	++LineSave.nsum;

	/* routince can be called with negative LineSave.ipass, in this case
	 * we are just counting number of lines for current setup */

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

	return;
}

/* this is energy in ryd as set by call to PntForLine */
static double EnergyRyd;
/* this is flag saying that EnergyRyd has been set */
static int lgEnergyRydSet=FALSE;

void lindst(double xInten, 
  float wavelength, 
  char *chLab, 
  long int ipnt, 
  char chInfo, 
  int lgOutToo)
{
	double e2fn;

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

	if( LineSave.ipass > 0 )
	{
		/* not first pass, sum lines only */
		ASSERT( xInten >= 0.);
		LineSv[LineSave.nsum].sumlin += xInten*radius.dVeff;
		LineSv[LineSave.nsum].emslin = xInten;

		/* add line to outward beam */
		/* there are lots of lines that are sums of other lines, or
		 * just for info of some sort.  These have this flag false.
		 * Note that the EnergyRyd variable only has a rational
		 * value if PntForLine was called just before this - in all
		 * cases where this did not happen the flag is false. */
		if( lgOutToo )
		{

			/* >>chng 97 sep 02, add tmn */
			/*rfield.outlin[ipnt-1] += (float)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
			  radius.dVolOutwrd*opac.tmn[ipnt-1]);*/
			/* >>chng 00 jan 01, from tmn to ExpZone, since now this is thrown into
			 * outward beam AFTER attenuation of field across zone.  This attenuation
			 * is still needed */
			rfield.outlin[ipnt-1] += (float)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
			  radius.dVolOutwrd*opac.ExpZone[ipnt-1]);
			rfield.reflin[ipnt-1] += (float)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
			  radius.dVolReflec);
		}
	}
	else if( LineSave.ipass == 0 )
	{
		/* >>chng 01 apr 16, remove assert, don't understand why it was here */
		/* confirm that either the energy in rydbergs was previously set,
		 * or we don't care 
		ASSERT( lgEnergyRydSet || !lgOutToo );*/
		/* if energy was not previously set in call to PntForLine then use
		 * array index ipnt that was sent down here, for energy */
		/* >>chng 01 apr 16, if energy not previously set use array index */
		if( !lgEnergyRydSet )
			EnergyRyd = rfield.anu[ipnt-1];

		LineSv[LineSave.nsum].xLineEnergy = (float)EnergyRyd;
		/* first call to suff lines in array, confirm that label is one of
		 * the four correct ones */
		/* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
		ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
		LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
		/* number of lines ok, set parameters for first pass */
		LineSv[LineSave.nsum].sumlin = 0.;
		LineSv[LineSave.nsum].emslin = 0.;
		LineSv[LineSave.nsum].wavelength = wavelength;

		strcpy( LineSv[LineSave.nsum].chALab, chLab );
	}

	/* this is dust part */
	if( gv.lgDustOn )
	{
		if( LineSave.ipass > 0 )
		{
			/* main production pass, sum lines only
			 * evaluate 2nd exponential integral function
			 * >>chng 97 nov 05, some variables not defined above nflux */
			if( ipnt <= rfield.nflux )
			{
				/* >>chng 01 aug 23, replace e2 with stored vlue */
				/*e2fn = e2(opac.TauAbsGeo[0][ipnt-1],opac.ExpmTau[ipnt-1]);*/
				e2fn = opac.e2TauAbs[ipnt-1];
				/* DSTOTH was unit vol energy lost onto dust by cooling lines
				 * dstoth = dstoth + xInten * (1.-e2fn) */
				LineDSv[LineSave.ndsum].smdlin += (float)(xInten*radius.dVeff/
				  2.*(1. + opac.albedo[ipnt-1])*e2fn);
			}
		}
		else if( LineSave.ipass == 0 )
		{
			/* set parameters for first pass */
			LineDSv[LineSave.ndsum].smdlin = 0.;
			LineDSv[LineSave.ndsum].wavelength = wavelength;
			strcpy( LineDSv[LineSave.ndsum].chSMDLab, chLab );
		}
		/* there is also case where ipass is negative, in this one we just add counter
		 * since space for array has not been malloced yet */
		/* increment the counter to the number of lines */
		++LineSave.ndsum;
	}

	/* increment the line counter */
	++LineSave.nsum;

	/* say that energy is not set */
	lgEnergyRydSet = FALSE;
	EnergyRyd = 0.;

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

	return;
}

/*PntForLine generate pointer for forbidden line */
void PntForLine(double wavelength, 
  char *chLabel,
  /* this is array index on the f, not c scale,
   * for the continuum cell holding the line */
  long int *ipnt)
{
	/* 
	 * maximum number of forbidden lines - this is a good bet since
	 * new lines do not go into this group, and lines are slowly 
	 * moving to level 1 
	 */
#	define	MAXFORLIN	1000
	static long int ipForLin[MAXFORLIN]={0};

	/* number of forbidden lines entered into continuum array */
	static long int nForLin;

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

	/* must be 0 or greater */
	ASSERT( wavelength >= 0. );

	if( wavelength == 0. )
	{
		/* zero is special flag to initialize */
		nForLin = 0;
		/* ipLineEnergy will only put in line label if nothing already there */
		EnergyRyd = 0.;
		/* say that energy is not set */
		lgEnergyRydSet = FALSE;
	}
	else
	{

		if( LineSave.ipass > 0 )
		{
			/* not first pass, sum lines only */
			*ipnt = ipForLin[nForLin];
		}
		else if( LineSave.ipass == 0 )
		{
			/* check if number of lines in arrays exceeded */
			if( nForLin >= MAXFORLIN )
			{
				fprintf( ioQQQ, "%5ld lines is too many for PntForLine.\n", 
				  nForLin );
				fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
				puts( "[Stop in pntforline]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* ipLineEnergy will only put in line label if nothing already there */
			EnergyRyd = RYDLAM/wavelength;
			/* say that energy is not set */
			lgEnergyRydSet = TRUE;
			ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
			*ipnt = ipForLin[nForLin];
		}
		else
		{
			/* this is case where we are only counting lines */
			*ipnt = 0;
		}
		++nForLin;
	}

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

	return;
}

static float ExtraInten;

/*PutLine enter local line intensity into the intensity stack for eventual printout */
void PutLine(EmLine * t)
{
	char chLabel[5];
	float wl;
	double xIntensity,
		other;

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

	/* routine to use line array data to generate input
	 * for emission line array */

	/* if ipass=0 then we must generate label info since first pass
	 * gt.0 then only need line intensity data */
	if( LineSave.ipass == 0 )
	{
		/* these variables not used by linadd if ipass ne 0 */
		/* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" */
		chIonLbl(chLabel, t);
		/* iWavLen converts wavelength into truncated form used in printout, ie,
		 * FIR lines are small integers, two-photon emission is 0 wavelength
		ASSERT( t->WLAng >= 0. );
		ASSERT( t->WLAng < 1e20 ); */

		/* the following converts wavelength into format for convenient printing,
		 * also returns units and possible shift of decimal point
		i = iWavLen(t, &chUnits , &chShift); */
		wl = t->WLAng;

		/* following is only non-zero for lines that are added into continuum */
		LineSv[LineSave.nsum].xLineEnergy = t->EnergyWN * (float)WAVNRYD;
		xIntensity = 0.;
	}
	else
	{
		/* both the counting and integrating modes comes here */
		strcpy( chLabel, "    " );
		wl = 0.;
		/* total line intensity or luminosity */
		xIntensity = t->xIntensity + ExtraInten;
	}
	/* initial counting case, where ipass == -1, just ignored above, call linadd below */

	/* ExtraInten is option that allows extra intensity (i.e., recomb)
	 * to be added to this line  with Call PutExtra( exta )
	 * in main lines this extra
	 * contribution must be identified explicitly */
	/* note that xLineEnergy is non-zero since defined above, for this call to linadd only,
	 * so this will be added to continuum.  Following calls to linadd will not be added to
	 * continuum since xLineEnergy will be left zero */
	ExtraInten = 0.;
	linadd(xIntensity,wl,chLabel,'i');

	/* inward part of line - do not move this away from previous lines
	 * since xIntensity is used here */
	xIntensity = xIntensity*t->FracInwd;
	linadd(xIntensity,wl,"Inwd",'i');

	/* cooling part of line */
	other = t->cool;
	linadd(other,wl,"Coll",'i');

	/* fluorescent excited part of line */
	other = t->PopOpc * t->pump * (1.-t->ColOvTot) * t->EnergyErg;
	linadd(other,wl,"Pump",'i');

	/* heating part of line */
	other = t->heat;
	linadd(other,wl,"Heat",'i');

	ASSERT( t->ipCont > 0. );
	/* this is dust part */
	if( gv.lgDustOn )
	{
		if( LineSave.ipass > 0 )
		{
			/* main production pass, sum lines only
			 * evaluate 2nd exponential integral function
			 * >>chng 97 nov 05, some variables not defined above nflux */
			if( t->ipCont <= rfield.nflux )
			{
				/* >>chng 01 aug 23, replace e2 with stored vlue */
				/*e2fn = e2(opac.TauAbsGeo[0][ipnt-1],opac.ExpmTau[ipnt-1]);*/
				/*float e2fn = opac.e2TauAbs[t->ipCont-1];*/
				/* DSTOTH was unit vol energy lost onto dust by cooling lines
				 * dstoth = dstoth + xInten * (1.-e2fn) */
				/* >>chng 01 aug 23, add to dust array, (had not been
				 * in before) also use inward fraction */
				LineDSv[LineSave.ndsum].smdlin += (float)(xIntensity*radius.dVeff*
				  (1. + opac.albedo[t->ipCont-1])*opac.e2TauAbs[t->ipCont-1]);
			}
		}
		else if( LineSave.ipass == 0 )
		{
			/* set parameters for first pass */
			LineDSv[LineSave.ndsum].smdlin = 0.;
			LineDSv[LineSave.ndsum].wavelength = wl;
			strcpy( LineDSv[LineSave.ndsum].chSMDLab, chLabel );
		}
		/* there is also case where ipass is negative, in this one we just add counter
		 * since space for array has not been malloced yet */
		/* increment the counter to the number of lines */
		++LineSave.ndsum;
	}

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

	return;
}

/* ==================================================================== */
/*PutExtra enter and 'extra' intensity source for some line */
void PutExtra(double Extra)
{

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

	ExtraInten = (float)Extra;

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

	return;
}

/*EmLineJunk set all elements of EmLine struc to dangerous values */
/*EmLineZero zeros out the emission line structure,
 * declarations are in taulines.h */
void EmLineJunk( EmLine * t )
{

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

	/* optical depth in continuum to ill face */
	t->TauCon = -FLT_MAX;

	/* inward and total line optical depths */
	t->TauIn = -FLT_MAX;
	t->TauTot = -FLT_MAX ;

	/* type of redistribution function, */
	t->iRedisFun = INT_MIN ;

	/* array offset for line within continuum array */
	t->ipCont = -10000;

	/* array offset for line within fine opacity */
	t->ipFine = -10000;

	/* inward fraction */
	t->FracInwd = -FLT_MAX;

	/* continuum pumping rate */
	t->pump = -FLT_MAX ;

	/* line intensity */
	t->xIntensity = -FLT_MAX ;

	/* number of photons emitted per sec in the line */
	t->phots = -FLT_MAX;

	/* gf value */
	t->gf = -FLT_MAX ;

	/* escape and destruction probs */
	t->Pesc = -FLT_MAX;
	t->Pdest = -FLT_MAX;
	t->Pelec_esc = -FLT_MAX;

	/* damping constant, and number related to it */
	t->dampXvel = -FLT_MAX ;
	t->damp = -FLT_MAX ;

	/* cooling and heating due to collisional excitation */
	t->cool = -FLT_MAX;
	t->heat = -FLT_MAX ;

	/* ratio of collisional to radiative excitation*/
	t->ColOvTot = -FLT_MAX ;

	/* collision strength for transition */
	t->cs = -FLT_MAX;

	 /* ion stage of element, 1 for atom, 2 ion, etc */
	t->IonStg = -10000;

	 /* atomic number of element, 1 for H, 2 for He, etc */
	t->nelem = -10000;

	/* wavelentgh, usually in A, used for printout */
	t->WLAng = -FLT_MAX;

	/* transition energy in degrees kelvin*/
	t->EnergyK = -FLT_MAX;
	/* transition energy in ergs */
	t->EnergyErg = -FLT_MAX;
	/* transition energy in wavenumbers */
	t->EnergyWN = -FLT_MAX;

	/* line opacity */
	t->opacity = -FLT_MAX;

	/* lower and upper statistical weights */
	t->gLo = -FLT_MAX;
	t->gHi = -FLT_MAX ;

	/* populations of lower and upper levels, and pop that enters net opacity */
	t->PopLo = -FLT_MAX ;
	t->PopHi = -FLT_MAX ;
	t->PopOpc = -FLT_MAX ;

	/* transition prob, Einstein A upper to lower */
	t->Aul = -FLT_MAX ;

	/* ots rate */
	t->ots = -FLT_MAX;

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

	return;
}

/*EmLineZero zeros out the emission line structure */
/*>>chng 03 feb 14, rm sets that damage the array - now only
 * zero quantities that need to be set to 0 at start of calc */
void EmLineZero( EmLine * t )
{

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

	/* total optical depth in all overlapping lines to illuminated face,
	 * used for pumping */
	t->TauCon = opac.taumin;

	/* inward and total line optical depths */
	/* >>chng 03 feb 14, from 0 to opac.taumin */
	t->TauIn = opac.taumin;

	/* total optical depths */
	t->TauTot = 1e20f;

	/* inward fraction */
	/* >>chng 03 feb 14, from 0 to 1 */
	t->FracInwd = 1.;

	/* continuum pumping rate */
	t->pump = 0. ;

	/* line intensity */
	t->xIntensity = 0. ;

	/* number of photons emitted per sec in the line */
	t->phots = 0.;

	/* escape and destruction probs */
	/* >>chng 03 feb 14, change from 0 to 1 */
	t->Pesc = 1.;
	t->Pdest = 0.;
	t->Pelec_esc = 0.;

	/* cooling and heating due to collisional excitation */
	t->cool = 0.;
	t->heat = 0. ;

	/* ratio of collisional to radiative excitation*/
	t->ColOvTot = 0. ;

	/* populations of lower and upper levels, and pop that enters net opacity */
	t->PopLo = 0. ;
	t->PopHi = 0. ;
	t->PopOpc = 0. ;

	/* ots rate */
	t->ots = 0.;

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

	return;
}

/*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
void LineConvRate2CS( EmLine * t , float rate )
{

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

	/* return is collision strength, convert from collision rate from 
	 * upper to lower, this assumes pure electron collisions, but that will
	 * also be assumed by anything that uses cs, for self-consistency */
	t->cs = rate * t->gHi / phycon.cdsqte;

	/* change assert to non-negative - there can be cases (Iin H2) where cs has
	 * underflowed to 0 on some platforms */
	ASSERT( t->cs >= 0. );

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

	return;
}

/* returns true is we have not overrun optical depth scale */
int lgTauGood( EmLine * t)
{
	int lgGoodTau;

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

	/*lint -e777 float test equality */
	if( (t->TauTot*0.9 - t->TauIn < 0. && t->TauIn > 0.) && 
		  t->TauTot != opac.taumin )
	/*lint +e777 float test equality */
	{
		/* do not do anything if we have overrun the scale, */
		lgGoodTau = FALSE;
	}
	else
	{
		lgGoodTau = TRUE;
	}

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

	return lgGoodTau;

}

/*gbar0 compute g-bar gaunt factor for neutrals */
static void gbar0(double ex, 
	  float *g)
{
	double a, 
	  b, 
	  c, 
	  d, 
	  y;

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

	/* written by Dima Verner
	 *
	 * Calculation of the effective Gaunt-factor by use of 
	 * >>refer	gbar	cs	Van Regemorter, H., 1962, ApJ 136, 906
	 * fits for neutrals
	 *  Input parameters: 
	 * ex - energy ryd - now K
	 * t  - temperature in K
	 *  Output parameter:
	 * g  - effective Gaunt factor
	 * */

	/* y = ex*157813.7/te */
	y = ex/phycon.te;
	if( y < 0.01 )
	{
		*g = (float)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
	}
	else
	{
		if( y > 10.0 )
		{
			*g = (float)(0.066/sqrt(y));
		}
		else
		{
			a = 1.5819068e-02;
			b = 1.3018207e00;
			c = 2.6896230e-03;
			d = 2.5486007e00;
			d = log(y/c)/d;
			*g = (float)(a + b*exp(-0.5*POW2(d)));
		}
	}

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

	return;
}

/*gbar1 compute g-bar collision strength using Mewe approximations */
static void gbar1(double ex, 
	  float *g)
{
	double y;

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

	/*	*written by Dima Verner
	 *
	 * Calculation of the effective Gaunt-factor by use of 
	 * >>refer	gbar	cs	Mewe,R., 1972, A&A 20, 215
	 * fits for permitted transitions in ions MgII, CaII, FeII (delta n = 0)
	 * Input parameters: 
	 * ex - excitation energy in Ryd - now K
	 * te  - temperature in K
	 * Output parameter:
	 * g  - effective Gaunt factor
	 */

	/* y = ex*157813.7/te */
	y = ex/phycon.te;
	*g = (float)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));

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

	return;
}

/*MakeCS compute collision strength by g-bar approximations */
void MakeCS(EmLine * t)
{
	long int ion;
	double Abun, 
	  cs;
	float
	  gbar;

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

	/* routine to get cs from various approximations */

	/* check if abundance greater than 0 */
	ion = t->IonStg;

	Abun = dense.xIonDense[ t->nelem -1 ][ ion-1 ];
	if( Abun <= 0. )
	{
		gbar = 1.;
	}
	else
	{
		/* check if neutral or ion */
		if( ion == 1 )
		{
			/* neutral - compute gbar for evantual cs */
			gbar0(t->EnergyK,&gbar);
		}
		else
		{
			/* ion - compute gbar for evantual cs */
			gbar1(t->EnergyK,&gbar);
		}
	}

	/* above was g-bar, convert to cs */
	cs = gbar*(14.5104/WAVNRYD)*t->gf/t->EnergyWN;

	/* stuff the cs in place */
	t->cs = (float)cs;

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

	return;
}

/*totlin sum total intensity of cooling, recombination, or intensity lines */
double totlin(
	/* chInfor is 1 char, 
	'i' information, 
	'r' recombination or 
	'c' collision */
	int chInfo)
{
	long int i;
	double totlin_v;

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

	/* routine goes through set of entered line
	 * intensities and picks out those which have
	 * types agreeing with chInfo.  Valid types are
	 * 'c', 'r', and 'i'
	 *begin sanity check */
	if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
	{
		fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n", 
		  chInfo );
		puts( "[Stop in totlin]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*end sanity check */

	/* now find sum of lines of type chInfo */
	totlin_v = 0.;
	for( i=0; i < LineSave.nsum; i++ )
	{
		if( LineSv[i].chSumTyp == chInfo )
		{
			totlin_v += LineSv[i].sumlin;
		}
	}

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

	return( totlin_v );
}


/*FndLineHt search through line heat arrays to find the strongest heat source */
void FndLineHt(long int *level, 
  /* this is the index of the strongest line in the array on the c scale */
  long int *ipStrong, 
  double *Strong)
{
	long int i; 

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

	*Strong = 0.;
	*level = 0;

	/* do the level 1 lines, 0 is dummy line, <=nLevel1 is correct for c scale */
	for( i=1; i <= nLevel1; i++ )
	{
		/* check if a line was the major heat agent */
		if( TauLines[i].heat > *Strong )
		{
			*ipStrong = i;
			*level = 1;
			*Strong = TauLines[i].heat;
		}
	}

	/* now do the level 2 lines */
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			/* check if a line was the major heat agent */
			if( TauLine2[i].heat > *Strong )
			{
				*ipStrong = i;
				*level = 2;
				*Strong = TauLine2[i].heat;
			}
		}
	}

	/* now do co carbon monoxide lines */
	for( i=0; i < nCORotate; i++ )
	{
		/* check if a line was the major heat agent */
		if( C12O16Rotate[i].heat > *Strong )
		{
			*ipStrong = i;
			*level = 3;
			*Strong = C12O16Rotate[i].heat;
		}
	}
	for( i=0; i < nCORotate; i++ )
	{
		/* check if a line was the major heat agent */
		if( C13O16Rotate[i].heat > *Strong )
		{
			*ipStrong = i;
			*level = 4;
			*Strong = C13O16Rotate[i].heat;
		}
	}

	/* now do the hyperfine structure lines */
	for( i=0; i < nHFLines; i++ )
	{
		/* check if a line was the major heat agent */
		if( HFLines[i].heat > *Strong )
		{
			*ipStrong = i;
			*level = 5;
			*Strong = HFLines[i].heat;
		}
	}

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

	return;
}

