/* 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_line_all do escape and destruction probs for all lines in code.  
 * Called with FALSE arg by ionize, to only redo destruction probabilites,
 * and with TRUE by cloudy to do both escape and destruction */
#include "cddefines.h"
#include "taulines.h"
#include "atomfeii.h"
#include "dense.h"
#include "converge.h"
#include "atoms.h"
#include "rfield.h"
#include "iso.h"
#include "h2.h"
#include "opacity.h"
#include "trace.h"
#include "lines_service.h"
#include "atmdat.h"
#include "hydrogenic.h"
#include "rt.h"

void RT_line_all(
	/* this is TRUE if we want to do both escape and destruction probs,
	 * and FALSE if only destruction probabilities are needed */
	int lgDoEsc ,
	/* flag saying whether to update fine opacities */
	int lgUpdateFineOpac )
{
	long int i,
		ion,
		ipISO,
		nelem;
	long ipHi , ipLo;
	double factor,
		coloi;
	int lgTOnSave;
	float SaveLyaPesc[NISO][LIMELM] , 
		SaveLyaPdest[NISO][LIMELM];

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

	if( trace.lgTrace )
		fprintf( ioQQQ, "     RT_line_all called\n" );

	/* skip line transfer if requested, and not first call in this zone */
	if( !rfield.lgDoLineTrans )
	{
		if( conv.nPres2Ioniz  )
			return;
	}

	/* this array is huge and takes significant time to zero out or update, only 
	 * do so when needed, */
	rfield.lgFine_opac_update = lgUpdateFineOpac;
	if( rfield.lgFine_opac_update )
	{
		/*fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac memset zero\n");*/
		memset(rfield.fine_opac , 0 , (unsigned long)rfield.nfine*sizeof(float) );
	}

	/* hydrogenic lyman lines special since outward optical depths always set,
	 * trick routines for Lyman lines only */
	lgTOnSave = opac.lgTauOutOn;
	opac.lgTauOutOn = TRUE;

#	if 0
	/* this code is a copy of the code within line_one that does cloaking
	 * within this zone.  it can be used to see how a particular line
	 * is being treated. */
	{
#include "doppvel.h"
		double dTau , aa;
		ipISO = 0; nelem = 0;ipLo = 0;
		ipHi = 23;
		dTau =  (float)EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc * 
			EmisLines[ipISO][nelem][ipHi][ipLo].opacity / 
			DoppVel.doppler[nelem] + (float)opac.opacity_abs[EmisLines[ipISO][nelem][ipHi][ipLo].ipCont-1];
		dTau *= radius.drad_x_fillfac_mean;
		aa = log(1. + dTau ) / SDIV(dTau);
		fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau,
			radius.drad_x_fillfac_mean,
			 aa);
	}
#	endif

	/*this is loop over h-like species - this is special case since 2s-2p are
	 * resolved but degenerate so must bring balmer line optical depths together */
	if( lgDoEsc )
	{
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			/* evaluate hydrogenic balance if ionization reaches this high */
			if( dense.IonHigh[nelem] == (nelem + 1) )
			{
				/* update As and opacities for higher n lines in h-like series */
				HLineTransOpacSet( nelem );
			}
		}
	}

	/* find Stark escape probabilities for hydrogen itself */
	RT_stark();

	/*if( lgUpdateFineOpac )
		fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* loop over all iso-electronic sequences */
		for( nelem=ipISO; nelem<LIMELM; ++nelem )
		{
			/* this is the parent ion, for H lines, is 1, 
			 * for He is 1 for he-like and 2 for h-like */
			ion = nelem+1-ipISO;
			if( ion <=dense.IonHigh[nelem] )
			{
				/* save escape and destruction probs for each Lya since
				 * these are special cases, and quantities are not damped
				 * in routine that evaluates each, test is to avoid
				 * evaluating iso sequence that does not exist for an element,
				 * as in He-like hydrogen */
				if( ipISO<=nelem )
				{
					SaveLyaPesc[ipISO][nelem] = EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pesc;
					SaveLyaPdest[ipISO][nelem] = EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pdest;
				}

				/* convert pops to per unit vol rather than per ion */
				if( dense.xIonDense[nelem][ion] > 1e-30 )
				{
					factor = dense.xIonDense[nelem][ion];
				}
				else
				{
					/* case where almost no parent ion - this will make
					* very large line opacity, so background dest small */
					factor = 1.;
				}
				/* all lyman line optical depths are define on first iteration,
				 * but only for H-like species */
				if( ipISO==ipH_LIKE )
						opac.lgTauOutOn = TRUE;
				/* loop over all lines */
				for( ipLo=0; ipLo < (iso.numLevels[ipISO][nelem]-1); ++ipLo )
				{
					for( ipHi=ipLo+1; ipHi < iso.numLevels[ipISO][nelem]; ++ipHi )
					{
						/* negative ipCont means this is not a real line, so do not
						 * transfer it */
						if( EmisLines[ipISO][nelem][ipHi][ipLo].ipCont < 0 ) 
							continue;

						/* must temporarily make ipLnPopOpc physical */
						EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc *= factor;

						/* generate escape prob, pumping rate, destruction prob, 
						* inward outward fracs */
						RT_line_one(&EmisLines[ipISO][nelem][ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,TRUE);
						{
							/* set true to print pump rates*/
							/*@-redef@*/
							enum {DEBUG_LOC=FALSE};
							/*@+redef@*/
							if( DEBUG_LOC && nelem==1&& ipLo==0  /*&& iteration==2*/ )
							{
								fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
									ipHi ,
									fnzone,
									EmisLines[ipISO][nelem][ipHi][ipLo].Pdest);
							}
						}

						/* go back to original units so that final correction ok */
						EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc /= factor;
					}

					/* reset this flag, so only Lyman lines forced to include outward */
					opac.lgTauOutOn = lgTOnSave;
				}
				/* >>chng 02 nov 18, add this cap to reduce outward Lya */
				/* >>chng 03 may 06, did not have lgDoEsc on this, so Lya Pesc
				 * became ever smaller as model was converged in one zone */
				if( ipISO > 0 && lgDoEsc )
				{
					/* don't let too much Lya escape outward since so important */
					EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pesc *= 
						opac.ExpmTau[EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].ipCont-1];
				}

				/* now update two photon induced rates */
				atmdat_2phot_rate(ipISO , nelem);

				/* this is option to not do destruction probabilities
				 * for case b no pdest option */
				if( opac.lgCaseB_no_pdest )
				{
					ipLo = 0;
					for( ipHi=ipLo+1; ipHi < iso.numLevels[ipISO][nelem]; ++ipHi )
					{
						/* hose the previously computed destruction probability */
						EmisLines[ipISO][nelem][ipHi][ipLo].Pdest = SMALLFLOAT; 
					}
				}

				ipLo = 0;
				/* these are the extra lyman lines for the iso sequences */
				/*for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )*/
				for( ipHi=iso.quant_desig[ipISO][nelem][iso.numLevels[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ )
				{
					/* must make ipLnPopOpc physical */
					iso.ExtraLymanLines[ipISO][nelem][ipHi].PopOpc = 
						/* just use this line to get population of ground state */
						EmisLines[ipISO][nelem][3][ipLo].PopOpc *factor;

					/* actually do the work */
					RT_line_one(&iso.ExtraLymanLines[ipISO][nelem][ipHi] , lgDoEsc , lgUpdateFineOpac,TRUE);

					/* reset to funny population of ground state */
					iso.ExtraLymanLines[ipISO][nelem][ipHi].PopOpc = 
						EmisLines[ipISO][nelem][3][ipLo].PopOpc;
				}

			}/* if nelem if ion <=dense.IonHigh */
		}/* loop over nelem */
	}/* loop over ipISO */

	/* note that pesc and dest are updated no matter what escprob logic we
	 * specify, and that these are not updated if we have overrun the
	 * optical depth scale.  only update here in that case */
	if(	lgTauGood( &EmisLines[ipH_LIKE][ipHYDROGEN][iso.nLyaLevel[ipH_LIKE]][0] ) )
	{
		if( lgDoEsc )
		{
			/* these loops exclude lya, so are only evaluated with escprob are
			 * evaluated.
			 * add on Stark escape probabilities for H itself */
			for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][ipHYDROGEN] - 1); ipLo++ )
			{
				/* >>chng 02 jun 12, do not do Lya here, since done above */
				for( ipHi=MAX2((long)3,ipLo+1); ipHi < iso.numLevels[ipH_LIKE][ipHYDROGEN]; ipHi++ )
				{
					/* >>chng 03 jun 07, do not clobber esp prob when line is masing -
					 * this had effect of preventing total escape prob from getting larger than 1 */
					if( EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Pesc<1. )
					{
						EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Pesc = (float)(MIN2(1.,
							EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Pesc+
							hydro.pestrk[ipLo][ipHi]));
					}
				}
			}
		}

		/* always add on stark broadening term to Lya*/
		if( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc < 1. )
		{
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc = (float)(MIN2(1.,
				EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc+
				hydro.pestrk[ipH1s][ipH2p]));
		}

		/* do the 8446 problem */
		atom_oi_calc(&coloi);

		EmisLines[ipH_LIKE][ipHYDROGEN][3][ipH1s].Pesc = atoms.pmph31/
			EmisLines[ipH_LIKE][ipHYDROGEN][3][ipH1s].Aul;

		if( trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE] )
		{
			fprintf( ioQQQ, "       RT_line_all calls P8446 who found pmph31=%10.2e\n", 
				atoms.pmph31 );
		}

		/* add on destruction of hydrogen Lya by FeII
		* now add in FeII deexcitation via overlap,
		* but only as loss, not photoionization, source
		* dstfe2lya is Ly-alpha deexcit by FeII overlap - conversion into FeII em */
		/* find FeII overlap destruction rate, 
		* this does NOTHING when large FeII atom is turned on, 
		* in this case evaluation is done in call to FeIILevelPops */
		atoms_fe2ovr();

		/* >>chng 00 jan 06, let dest be large than 1 to desaturate the atom */
		/* >>chng 01 apr 01, add test for tout >= 0., 
		 * must not add to Pdest when it has not been refreshed here */
		EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest += hydro.dstfe2lya;
	}

	/* take mean of old and new to damp oscillations in Lya (only) */
	if( nzone > 1 )
	{
		for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
		{
			/* loop over all iso-electronic sequences */
			for( nelem=ipISO; nelem<LIMELM; ++nelem )
			{
				ion = nelem+1-ipISO;
				if( ion <=dense.IonHigh[nelem] && ipISO<=nelem )
				{
					/* now take mean of old and new escape dest probs for Lya */
					/* only ipLY_A redist did not already take average */
					if( EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].iRedisFun==ipLY_A &&
						lgTauGood( &EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) )
					{
						/* the Lya routine is different in that escape and dest
						 * probs are always updated, unless we have overrun the
						 * optical depth scale */
						/*lint -e644 SaveLyaPdest not init */
						EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pdest = 
							(SaveLyaPdest[ipISO][nelem] + 
							EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pdest) / 2.f;
						EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pesc = 
							(SaveLyaPesc[ipISO][nelem] + 
							EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Pesc) / 2.f;
						/*lint +e644 SaveLyaPdest not init */
					}
				}
			}
		}
	}
	/*5986 fprintf(ioQQQ,"DEBUG he rt_all\t%li\t%li\t%.4e\t%.4e\n",
		iteration, nzone, 
				EmisLines[1][1][10][5].TauIn,
				EmisLines[1][1][10][5].Pesc);*/

	{
		/* following should be set true to print ots contributors for he-like lines*/
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ )
		{
			/* option to dump a line  */
			DumpLine(&EmisLines[ipH_LIKE][ipHYDROGEN][2][0] );
#			if 0
			fprintf(ioQQQ,"debugggg\t%.3f\t%i\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				fnzone,
				lgDoEsc ,
				EmisLines[ipH_LIKE][ipHYDROGEN][2][0].Pesc ,
				EmisLines[ipH_LIKE][ipHYDROGEN][2][0].Pdest,
				EmisLines[ipH_LIKE][ipHYDROGEN][2][0].pump ,
				EmisLines[ipH_LIKE][ipHYDROGEN][2][0].TauCon ,
				EmisLines[ipH_LIKE][ipHYDROGEN][2][0].TauIn ,
				EmisLines[ipH_LIKE][ipHYDROGEN][2][0].TauTot,
				EmisLines[ipH_LIKE][ipHYDROGEN][3][2].Pesc ,
				EmisLines[ipH_LIKE][ipHYDROGEN][3][2].Pdest ,
				EmisLines[ipH_LIKE][ipHYDROGEN][3][2].pump );
#			endif
		}
	}

	/* level 1 lines */
	for( i=1; i <= nLevel1; i++ )
	{
		RT_line_one(&TauLines[i] , lgDoEsc , lgUpdateFineOpac,TRUE);
	}
	/* co carbon monoxide */
	for( i=0; i < nCORotate; i++ )
	{
		RT_line_one(&C12O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,TRUE);
		RT_line_one(&C13O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,TRUE);
	}
	for( i=0; i < nHFLines; i++ )
	{
		RT_line_one(&HFLines[i] , lgDoEsc , lgUpdateFineOpac,TRUE);
	}
	for( i=0; i < nUTA; i++ )
	{
		if( UTALines[i].Aul > 0. )
		{
			/* these are not defined in cooling routines */
			UTALines[i].PopOpc = dense.xIonDense[UTALines[i].nelem-1][UTALines[i].IonStg-1];
			UTALines[i].PopLo = dense.xIonDense[UTALines[i].nelem-1][UTALines[i].IonStg-1];
			UTALines[i].PopHi = 0.;
			RT_line_one(&UTALines[i] , lgDoEsc , lgUpdateFineOpac,TRUE);
		}
	}

	/* >>chng 03 aug 22, must call RT_line_one every single time, since that routine
	 * establishes the fine opacities for the level 2 lines */
	/* >>chng 01 aug 11, the level 2 lines were only updated when lgDoEsc was true,
	 * this meant that dest probs were not updated but once, causing very high-Z models
	 * to develop numerical oscialltions.  removed if, now evaluate level 2 always */
	/* only update their dest/esc probs one time per zone */
	/* lgLevel2_OTS_Imp set true in dimacool if ots rates were significant */
	/* level 2 heavy element lines in cooling with only g-bar,*/
	for( i=0; i < nWindLine; i++ )
	{
		/* >>chng 02 sug 11, do not include ions done in iso-sequences */
		/*if( TauLine2[i].IonStg <= TauLine2[i].nelem-1 )*/
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			RT_line_one(&TauLine2[i] , lgDoEsc , lgUpdateFineOpac,TRUE);
		}
	}

	/* the large H2 molecule */
	H2_RTMake( lgDoEsc , lgUpdateFineOpac);

	/* The large model FeII atom */
	FeIIRTMake( lgDoEsc , lgUpdateFineOpac);

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

