/* 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 */
/*RT_tau_init set initial outward optical depths at start of first iteration,
 * it is only called by cloudy one time per complete calculation, just after
 * continuum set up and before start of convergence attempts. 
 * 
 * RT_tau_reset after first iteration, updates the optical depths, mirroring this
 * routine but with the previous iteration's variables */
#include "cddefines.h"
#define	TAULIM	1e8
#include "taulines.h"
#include "trace.h"
#include "doppvel.h"
#include "iso.h"
#include "h2.h"
#include "lines_service.h"
#include "rfield.h"
#include "dense.h"
#include "opacity.h"
#include "thermal.h"
#include "phycon.h"
#include "geometry.h"
#include "stopcalc.h"
#include "ipoint.h"
#include "abund.h"
#include "converge.h"
#include "tfidle.h"
#include "atomfeii.h"
#include "rt.h"

void RT_tau_init(void)
{
	long int i, 
	  nelem,
	  ipISO,
	  ipHi, 
	  ipLo;
	long lgHit; /* used for logic check */

	double AbunRatio, 
	  balc, 
	  coleff, 
	  f, 
	  /* z4 ,used for saving nelem+1**4 */
	  BalmerTauOn;

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

	/* >>chng 03 apr 11, move iso together, also move all EmLineZero to top of
	 * this routine */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				for( ipHi=1; ipHi < iso.numLevels[ipISO][nelem]; ipHi++ )
				{
					for( ipLo=0; ipLo < ipHi; ipLo++ )
					{
						/* >>chng 03 feb 14, use EmLineZero rather than explict sets */
						EmLineZero( &EmisLines[ipISO][nelem][ipHi][ipLo] );
					}
				}
				for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ )
				{
					/* >>chng 03 feb 14, use EmLineZero rather than explict sets */
					EmLineZero( &iso.ExtraLymanLines[ipISO][nelem][ipHi] );
				}
			}
		}
	}
#	if 0
	/* initial optical depth for he-like sequence,
	 * do not try to set outward looking optical depths */
	for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ipHi=1; ipHi < iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
			{
				for( ipLo=0; ipLo < ipHi; ipLo++ )
				{
					/* >>chng 03 feb 14, use EmLineZero rather than explict sets */
					EmLineZero( &EmisLines[ipHE_LIKE][nelem][ipHi][ipLo] );
				}
			}
		}
	}

	/* do optical depths in extra Lyman lines */
	for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ )
				{
					/* >>chng 03 feb 14, use EmLineZero rather than explict sets */
					EmLineZero( &iso.ExtraLymanLines[ipISO][nelem][ipHi] );
				}
			}
		}
	}

	/* >>chng 03 apr 09, add this zero out */
	/* initial optical depth for h-like sequence,
	 * below will set outward looking optical depths */
	for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ipHi=1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				for( ipLo=0; ipLo < ipHi; ipLo++ )
				{
					/* >>chng 03 feb 14, use EmLineZero rather than explict sets */
					EmLineZero( &EmisLines[ipH_LIKE][nelem][ipHi][ipLo] );
				}
			}
		}
	}
#	endif
	/* initialize heavy element line array */
	for( i=1; i <= nLevel1; i++ )
	{
		/* >>chng 03 feb 14, use EmLineZero rather than explict sets */
		EmLineZero( &TauLines[i] );
	}
	/* this is a dummy optical depth array for non-existant lines 
	 * when this goes over to struc, make sure all are set to zero here since
	 * init in grid may depend on it */
	EmLineZero( &TauDummy );
	
	/* lines in cooling function with Mewe approximate collision strengths */
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			/* inward optical depth */
			EmLineZero( &TauLine2[i] );
		}
	}
	
	/* inner shell lines */
	for( i=0; i < nUTA; i++ )
	{
		/* these are line optical depth arrays
		 * inward optical depth */
		EmLineZero( &UTALines[i] );
	}
	/* hyperfine structure lines */
	for( i=0; i < nHFLines; i++ )
	{
		EmLineZero( &HFLines[i] );
	}
	
	/* CO carbon monoxide lines */
	for( i=0; i < nCORotate; i++ )
	{
		/* >>chng 03 feb 14, use main zero function */
		EmLineZero( &C12O16Rotate[i] );
		EmLineZero( &C13O16Rotate[i] );
	}

	/* initialize optical depths in H2 */
	H2_LineZero();

	/* initialize large atom FeII arrays */
	FeII_LineZero();

	/* ==================================================================*/
	/* end setting lines to zero */

	/* >>chng 02 feb 08, moved to here from opacitycreateall, which was called later.
	 * bug inhibited convergence in some models.  Caught by PvH */
	/* this is option to stop at certain optical depth */
	if( StopCalc.taunu > 0. )
	{
		StopCalc.iptnu = ipoint(StopCalc.taunu);
		StopCalc.iptnu = MIN2(StopCalc.iptnu,rfield.nupper-1);
	}
	else
	{
		StopCalc.iptnu = rfield.nupper;
		StopCalc.tauend = 1e30f;
	}

	/* if taunu was set with a stop optical depth command, then iptnu must
	 * have also been set to a continuum index before this code is reaced */
	ASSERT( StopCalc.taunu == 0. || StopCalc.iptnu >= 0 );

	/* set initial and total optical depths in arrays
	 * TAUNU is set when lines read in, sets stopping radius */
	if( StopCalc.taunu > 0. )
	{
		/*  an optical depth has been specified */
		if( StopCalc.iptnu >= iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] )
		{
			/* at ionizing energies */
			for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
			{
				/* taumin set to 1e-20, can be reset with taumin command */
				opac.TauAbsGeo[1][i] = opac.taumin;
				opac.TauScatGeo[1][i] = opac.taumin;
				opac.TauTotalGeo[1][i] = opac.taumin;
			}

			for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
			{
				/* TauAbsGeo(i,2) = tauend * (anu(i)/anu(iptnu))**(-2.43) */
				opac.TauAbsGeo[1][i] = StopCalc.tauend;
				opac.TauScatGeo[1][i] = opac.taumin;
				opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
			}
		}

		else
		{
			/* not specified at ionizing energies */
			for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
			{
				opac.TauAbsGeo[1][i] = StopCalc.tauend;
				opac.TauScatGeo[1][i] = StopCalc.tauend;
				opac.TauTotalGeo[1][i] = StopCalc.tauend;
			}

			for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
			{
				opac.TauAbsGeo[1][i] = (float)(TAULIM*pow(rfield.anu[i],-2.43f));
				opac.TauScatGeo[1][i] = opac.taumin;
				opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
			}

		}
	}

	else
	{
		/*  ending optical depth not specified, assume 1E8 at 1 Ryd */
		for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
		{
			opac.TauAbsGeo[1][i] = opac.taumin;
			opac.TauScatGeo[1][i] = opac.taumin;
			opac.TauTotalGeo[1][i] = opac.taumin;
		}

		for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
		{
			opac.TauAbsGeo[1][i] = (float)(TAULIM*pow(rfield.anu[i],-2.43f));
			opac.TauScatGeo[1][i] = opac.taumin;
			opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
		}
	}

	/* if lgSphere then double outer, set inner to half
	 * assume will be optically thin at sub-ionizing energies */
	if( geometry.lgSphere || opac.lgCaseB )
	{
		for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
		{
			opac.TauAbsGeo[0][i] = opac.taumin;
			opac.TauAbsGeo[1][i] = opac.taumin*2.f;
			opac.TauScatGeo[0][i] = opac.taumin;
			opac.TauScatGeo[1][i] = opac.taumin*2.f;
			opac.TauTotalGeo[0][i] = 2.f*opac.taumin;
			opac.TauTotalGeo[1][i] = 4.f*opac.taumin;
		}

		for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
		{
			opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i];
			opac.TauAbsGeo[1][i] *= 2.;
			opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i];
			opac.TauScatGeo[1][i] *= 2.;
			opac.TauTotalGeo[0][i] = opac.TauTotalGeo[1][i];
			opac.TauTotalGeo[1][i] *= 2.;
		}

		if( StopCalc.taunu > 0. )
		{
			/* ending optical depth specified, and lgSphere also */
			StopCalc.tauend *= 2.;
		}
	}

	/* fix escape prob for He, metals, first set log of Te, needed by RECEFF
	 * do not do this if temperature has been set by command */
	if( !thermal.lgTSetOn )
	{
		/* this is a typical temperature for the H+ zone, and will use it is
		 * the line widths & opt depth in the following.  this is not meant to be the first
		 * temperature during the search phase. */
		phycon.te = 2e4;

		/* propagate this temperature through the code */
		/* must not call tfidle at this stage since not all vectors have been allocated */
		/*tfidle(FALSE);*/
	}

	/* set inward optical depths for hydrogenic ions to small number proportional to abundance */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* will need this for some scaling laws - physical Z to the 4th power
		z4 = POW2(nelem+1.);
		z4 *= z4;*/

		if( dense.lgElmtOn[nelem] )
		{
			HLineTransOpacSet( nelem );
#			if 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 = 0.32 - 0.07*dense.eden/(dense.eden + 1e7);
			/* upper limit for HydroBranch */
			limit = MIN2(15,iso.numLevels[ipH_LIKE][nelem]);
			/* first reset line opacities and branching ratios */
			for( ipHi=4; ipHi < limit; ipHi++ )
			{
				EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul = 
					(float)(hydro.HyLife[ipHi]*
					factor*HydroBranch(ipHi,2,nelem+1)*z4);
				ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul > 0.);

				/* do 2p in terms of 2s */
				EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul = (float)(
					EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul / factor * (1. - factor ) ); 
				ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul > 0.);

				/* make self-consistent opaciity */
				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)));
				ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].opacity > 0.);

				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)));
				ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].opacity > 0. );

				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)));
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity > 0.);
				}
			}
#			endif
			/* set in lines_service.h */
#			if LOWDEN_LYMAN
			/* for fix Lyman lines, in terms of alpha transition */
			for( ipHi=3; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
					float Ratio_lyman_alpha_Z1[25]={
						0.,0.,0.,
5.52E-01f,/* this mult by 1.33 since Ha split into two */
3.30E-01f,
2.96E-01f,
2.80E-01f,
2.74E-01f,
2.72E-01f,
2.72E-01f,
2.74E-01f,
2.76E-01f,
2.78E-01f,
2.81E-01f,
2.84E-01f,
2.86E-01f,
2.89E-01f,
2.92E-01f,
2.95E-01f,
2.98E-01f,
3.01E-01f,
3.05E-01f,
3.09E-01f,
3.13E-01f,
3.18E-01f};
					float Ratio_lyman_alpha_Z2[25]={
						0.,0.,0.,
4.52E-01f,
2.38E-01f,
1.98E-01f,
1.80E-01f,
1.71E-01f,
1.66E-01f,
1.64E-01f,
1.63E-01f,
1.63E-01f,
1.64E-01f,
1.65E-01f,
1.66E-01f,
1.68E-01f,
1.69E-01f,
1.71E-01f,
1.72E-01f,
1.73E-01f,
1.75E-01f,
1.76E-01f,
1.78E-01f,
1.79E-01f,
1.80E-01f};

						if( nelem==0 )
						{
							EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].Aul = EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].Aul *
							Ratio_lyman_alpha_Z1[MIN2(23,ipHi) ];
						}
						else
						{
							EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].Aul = EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].Aul *
							Ratio_lyman_alpha_Z2[MIN2(23,ipHi) ];
						}

						/* derive the abs coef, call to function is gf, wl (A), g_low */
						EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity = 
							(float)(abscf(
						  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].gf,
						  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].EnergyWN,
						  iso.stat[ipH_LIKE][nelem][ipH1s]));
			}
#			endif

			/* now get actual optical depths */
			AbunRatio = abund.solar[nelem]/abund.solar[0];
			for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					/* set all inward optical depths to taumin, regardless of abundance
					 * this is a very small number, 1e-20 */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn = opac.taumin;
				}
			}

			/* La may be case B, tlamin set to 1e9 by default with case b command */
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn = opac.tlamin;

			/* scale factor so that alll other lyman lines are appropriate for this Lya optical depth*/
			f = opac.tlamin/EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity;
			for( ipHi=3; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn = 
					(float)(f*EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity);
					ASSERT( EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn > 0. );
			}

			/* after this set of if's the total lya optical depth will be known,
			 * common code later on will set rest of lyman lines
			 * if case b then set total lyman to twice inner */
			if( opac.lgCaseB )
			{
				/* force outer optical depth to twice inner if case B */
				EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot = 
					(float)(2.*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn);
				/* force off Balmer et al optical depths */
				BalmerTauOn = 0.;
			}

			else
			{
				/* usual case for H LyA; try to guess outer optical depth */
				if( StopCalc.colnut < 6e29 )
				{
					/* neutral column is defined */
					EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot = (float)(StopCalc.colnut*
					  rt.DoubleTau*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity/
					  DoppVel.doppler[nelem]*AbunRatio);
					ASSERT( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot > 0. );
				}
				/* has optical depth at some energy where we want to stop been specified?
				 * taunu is energy where
				 * this has been specified - this checks on Lyman continuum, taunu = 1 */
				else if( StopCalc.taunu < 3. && StopCalc.taunu >= 0.99 )
				{
					/* Lyman continuum optical depth defined */
					EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot = (float)(StopCalc.tauend*
					  1.2e4*1.28e6/DoppVel.doppler[nelem]*rt.DoubleTau*AbunRatio);
					ASSERT( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot > 0. );
				}
				else if( StopCalc.HColStop < 6e29 )
				{
					/* neutral column not defined, guess from total col and U */
					coleff = StopCalc.HColStop - 
						MIN2(MIN2(rfield.qhtot/dense.eden,1e24)/2.6e-13,0.6*StopCalc.HColStop);

					EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot = (float)(coleff*
					  7.6e-14*AbunRatio);
					ASSERT( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot > 0. );
				}
				else
				{
					/* no way to estimate 912 optical depth, set to large number */
					EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot = (float)(1e20*
					  AbunRatio);
					ASSERT( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot > 0. );
				}
				/* allow Balmer et al. optical depths */
				BalmerTauOn = 1.;
			}

			/* Lya total optical depth now known, is it optically thick?*/
			if( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot > 1. )
			{
				rt.TAddHLya = (float)MIN2(1.,EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot/
				  1e4);
				EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn += rt.TAddHLya;
			}

			else
			{
				rt.TAddHLya = opac.tlamin;
			}

			/* this scale factor is to set other lyman lines, given the Lya optical depth */
			f = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot/
				EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity;

			for( ipHi=3; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				/* set total optical depth for higher lyman lines */
				EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauTot = 
					(float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity* f);
				ASSERT( EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauTot > 0.);

				/* increment inward optical depths by rt all lyman lines, inward
				 * optical depth was set above, this adds to it.  this was originally
				 * set some place above */
				EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn += rt.TAddHLya*
				  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity/
				  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity;
				ASSERT( EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn > 0.);
			}

			/* try to guess what Balmer cont optical guess,
			 * first branch is case where we will stop at a balmer continuum optical
			 * depth - taunu is energy where tauend was set */
			if( StopCalc.taunu > 0.24 && StopCalc.taunu < 0.7 )
			{
				EmisLines[ipH_LIKE][nelem][3][2].TauTot = (float)(StopCalc.tauend*
				  3.7e4*BalmerTauOn*AbunRatio + 1e-20);
			}

			else
			{
				/* this is a guess based on Ferland&Netzer 1979, but it gets very large */
				balc = rfield.qhtot*2.1e-19*BalmerTauOn*AbunRatio + 1e-20;

				EmisLines[ipH_LIKE][nelem][3][2].TauTot = 
					(float)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
				ASSERT( EmisLines[ipH_LIKE][nelem][3][2].TauTot > 0.);
			}

			/* 2s - 2p strongly forbidden */
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].TauTot = 2.f*opac.taumin;
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].TauIn = opac.taumin;

			/* 2s - 1s strongly forbidden */
			EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].TauTot = 2.f*opac.taumin;
			EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].TauIn = opac.taumin;

			/* fill in rest of the alpha transitions */
			if( iso.numLevels[ipH_LIKE][nelem] > 4 )
				EmisLines[ipH_LIKE][nelem][4][3].TauTot = 
				(float)(EmisLines[ipH_LIKE][nelem][3][2].TauTot*0.01*BalmerTauOn + opac.taumin);

			if( iso.numLevels[ipH_LIKE][nelem] > 5 )
			EmisLines[ipH_LIKE][nelem][5][4].TauTot = 
				(float)(EmisLines[ipH_LIKE][nelem][3][2].TauTot* 0.0005*BalmerTauOn + opac.taumin);

			if( iso.numLevels[ipH_LIKE][nelem] > 6 )
			EmisLines[ipH_LIKE][nelem][6][5].TauTot = 
				(float)(EmisLines[ipH_LIKE][nelem][3][2].TauTot*7.7e-5*BalmerTauOn + opac.taumin);

			for( ipHi=7; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].TauTot = 
					(float)(EmisLines[ipH_LIKE][nelem][6][5].TauTot/ipHi*BalmerTauOn + opac.taumin);
			}

			/* transitions down to 2s are special since 'alpha' (2s-2p) has
			 * small optical depth, so use 3 to 2p instead */
			f = EmisLines[ipH_LIKE][nelem][3][ipH2p].TauTot/
				EmisLines[ipH_LIKE][nelem][3][ipH2p].opacity* BalmerTauOn;

			ipLo = ipH2s;
			for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot = (float)(opac.taumin +
					f* EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity);
				ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot > 0.);
			}
#			if 0
			if( nelem==0 )
				fprintf(ioQQQ," h opacityyy %.5e %.5e %.5e %.5e %.5e %.5e\n",
				EmisLines[ipH_LIKE][nelem][3][ipH2p].TauTot ,
				EmisLines[ipH_LIKE][nelem][3][ipH2p].opacity, 
				BalmerTauOn ,
				EmisLines[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1][1].opacity,
				opac.taumin,
				EmisLines[ipH_LIKE][nelem][iso.numLevels[ipH_LIKE][nelem]-1][1].TauTot);
#			endif

			/* this is for rest of lines, scale from opacity */
			for( ipLo=ipH2p; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
			{
				f = EmisLines[ipH_LIKE][nelem][ipLo+1][ipLo].TauTot/
				  EmisLines[ipH_LIKE][nelem][ipLo+1][ipLo].opacity*
				  BalmerTauOn;

				for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot = (float)(opac.taumin +
						f* EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity);
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot > 0.);
				}
			}

			/* this loop is over all possible H lines, do some final cleanup */
			for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					/* TauCon is line optical depth to inner face used for continuum pumping rate,
					 * not equal to TauIn in case of static sphere since TauCon will never
					 * count the far side line optical depth */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauCon = 
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn;

					/* make sure inward optical depth is not larger than half total */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn = 
						MIN2(	EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn ,
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot/2.f );
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn > 0.);

					/* this is fraction of line that goes inward, not known until second iteration*/
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].FracInwd = 0.5;
				}
			}
		}
	}

	/* now do helium like sequence if case b */
	if( opac.lgCaseB )
	{
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				float Aprev,
					ratio;
				/* La may be case B, tlamin set to 1e9 by default with case b command */
				EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn = opac.tlamin;

				EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauCon = EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;

				EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauTot = 
				  2.f*EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;

				ratio = opac.tlamin/EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].opacity;

				/* this will be the trans prob of the previous lyman line, will use this to 
				 * find the next one up in the series */
				Aprev = EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Aul;

				i = ipHe2p1P+1;
				/* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess
				 * which are which - this will work for any number of levels */
				while( i < iso.numLevels[ipHE_LIKE][nelem] )
				/*while( i < N_HE_LYMAN && ipHeLyman[i] < iso.numLevels[ipHE_LIKE][nelem] )*/
				{
					/* >>chng 02 mar 15, add test for collapsed levels - As are very much
					 * smaller at boundary between collapsed/resolved so this test is needed*/
					if( EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul> Aprev/10. ||
						iso.quant_desig[ipHE_LIKE][nelem][i].s < 0 )
					{
						Aprev = EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul;

						EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn = 
							ratio*EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].opacity;
						/* reset line optical depth to continuum source */
						EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauCon = EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn;
						EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauTot = 
						  2.f*EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;
						/*fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",nelem, i, 
							EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul, EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn );*/
					}
					++ i;
				}
			}
		}
	}

	/* begin sanity check, total greater than 1/0.9 time inward */
	lgHit = FALSE;
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
				{
					if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot < 
						0.9f*EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn )
					{
						fprintf(ioQQQ,
							"RT_tau_init insanity for h line, Z=%li lo=%li hi=%li tot=%g in=%g \n",
							nelem , ipLo, ipHi , EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot , 
							EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn );
						lgHit = TRUE;
					}
				}
			}
		}
	}
	if( lgHit )
	{
		fprintf( ioQQQ," stopping due to insanity in RT_tau_init\n");
		ShowMe();
#		ifdef DEBUG_FUN
		fputs( " <->RT_tau_init()\n", debug_fp );
#		endif
		cdEXIT(EXIT_FAILURE);
	}
	/*end sanity check */

	/* fix offset for effective column density optical depth */
	rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];

	/* this is flag detected by dest prob routines to see whether ots rates are
	 * oscaillating - damp them out if so */
	conv.lgOscilOTS = FALSE;

	/* now that optical depths have been incremented, do escape prob, this
	 * is located here instead on in cloudy.c (where it belongs) because
	 * information generated by RT_line_all is needed for the following printout */

	RT_line_all(TRUE , FALSE );

	/* rest of routine is printout in case of trace */
	if( trace.lgTrace )
	{
		if( trace.lgHBug && trace.lgIsoTraceFull[ipH_LIKE] )
		{
			fprintf( ioQQQ, "\n\n   up EmisLines[ipH_LIKE][hi][lo].TauTot array as set in RT_tau_init ipZTrace (nelem)= %ld\n", 
			  trace.ipIsoTrace[ipH_LIKE] );
			for( ipHi=2; ipHi < iso.numLevels[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
			{
				fprintf( ioQQQ, " %3ld", ipHi );
				for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e",
						EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].TauTot ));
				}
				fprintf( ioQQQ, "\n" );
			}

			fprintf( ioQQQ, "\n\n TauIn array\n" );
			for( ipHi=1; ipHi < iso.numLevels[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
			{
				fprintf( ioQQQ, " %3ld", ipHi );
				for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e", 
						EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].TauIn ));
				}
				fprintf( ioQQQ, "\n" );
			}

			fprintf( ioQQQ, "\n\n Aul As array\n" );
			for( ipHi=1; ipHi < iso.numLevels[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
			{
				fprintf( ioQQQ, " %3ld", ipHi );
				for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e", 
						EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Aul) );
				}
				fprintf( ioQQQ, "\n" );
			}

			fprintf( ioQQQ, "\n\n Aul*esc array\n" );
			for( ipHi=1; ipHi < iso.numLevels[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
			{
				fprintf( ioQQQ, " %3ld", ipHi );
				for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e", 
						EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Aul*
					  (EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Pdest + 
					  EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Pelec_esc +
					  EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Pesc) ));
				}
				fprintf( ioQQQ, "\n" );
			}

			fprintf( ioQQQ, "\n\n H opac array\n" );
			for( ipHi=1; ipHi < iso.numLevels[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
			{
				fprintf( ioQQQ, " %3ld", ipHi );
				for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e", 
						EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].opacity ));
				}
				fprintf( ioQQQ, "\n" );
			}
		}

		else
		{
			fprintf( ioQQQ, " RT_tau_init called.\n" );
		}
	}
	ASSERT( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn> 0. );
	ASSERT( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc>= 0. );

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

