/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*RT_line_all do escape and destruction probabilities for all lines in code. */
#include "cddefines.h"
#include "taulines.h"
#include "atomfeii.h"
#include "dense.h"
#include "conv.h"
#include "atoms.h"
#include "rfield.h"
#include "wind.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"
#include "cosmology.h"
#include "physconst.h"
#include "doppvel.h"
#include "mole.h"

/*RT_line_all do escape and destruction probabilities for all lines in code. */
void RT_line_all( void )
{
	long int i,
		ion,
		ipISO,
		nelem;
	long ipHi , ipLo;

	/* this flag says whether to update the line escape probabilities */
	bool lgPescUpdate = conv.lgFirstSweepThisZone || conv.lgIonStageTrimed;

	DEBUG_ENTRY( "RT_line_all()" );

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

	/* rfield.lgDoLineTrans - skip line transfer if requested with no line transfer
	 * conv.nPres2Ioniz is zero during search phase and on first call in this zone */
	if( !rfield.lgDoLineTrans && conv.nPres2Ioniz )
	{
		return;
	}

	/* this array is huge and takes significant time to zero out or update, 
	 * only do so when needed, */
	if( conv.lgLastSweepThisZone )
	{
		/* zero out fine opacity array */
		memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );

		if( 0 && cosmology.lgDo )
		{
			realnum dVel = (realnum)SPEEDLIGHT * ( 1.f - 1.f/POW2(1.f+cosmology.redshift_start-cosmology.redshift_current) );
			rfield.ipFineConVelShift = -(long int)( dVel/ rfield.fine_opac_velocity_width + 0.5 );
		}
		else
		{
			/* this is offset within fine opacity array caused by Doppler shift
			 * between first zone, the standard of rest, and current position 
			 * in case of acceleration away from star in outflowing wind 
			 * dWind is positive, this means that the rest frame original
			 * velocity is red shifted to lower energy */
			realnum dWind = wind.windv - wind.windv0;

			/* will add ipVelShift to index of original energy so redshift has 
			 * to be negative */
			rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5);
		}
	}

#	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 =  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() * 
			iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() / 
			GetDopplerWidth(dense.AtomicWeight[nelem]) + opac.opacity_abs[iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1];
		dTau *= radius.drad_x_fillfac;
		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,
			 aa);
	}
#	endif

	/* find Stark escape probabilities for hydrogen itself */
	if( lgPescUpdate )
		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 )
		{
			/* parent ion stage, for H is 1, for He is 1 for He-like and 
			 * 2 for H-like */
			ion = nelem+1-ipISO;

			/* element turned off */
			if( !dense.lgElmtOn[nelem] )
				continue;
			/* need we consider this ion? */
			if( ion <= dense.IonHigh[nelem] )
			{
				/* loop over all lines */
				for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ++ipHi )
				{
					for( ipLo=0; ipLo < ipHi; ++ipLo )
					{
						/* negative ipCont means this is not a real line, so do not
						 * transfer it */
						if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() < 0 ) 
							continue;

						/* generate escape prob, pumping rate, destruction prob, 
						 * inward outward fracs */
						fixit();  // should this use pestrk_up or pestrk?
						RT_line_one( iso_sp[ipISO][nelem].trans(ipHi,ipLo),
							     true,(realnum)iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up,
								 GetDopplerWidth(dense.AtomicWeight[nelem]));

						/* set true to print pump rates*/
						enum {DEBUG_LOC=false};
						if( DEBUG_LOC && nelem==1&& ipLo==0  /*&& iteration==2*/ )
						{
							fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
								ipHi ,
								fnzone,
								iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest());
						}
					}
				}

				/* this is option to not do destruction probabilities
				 * for case b no pdest option */
				if( opac.lgCaseB_no_pdest )
				{
					ipLo = 0;
					/* okay to let this go to numLevels_max. */
					for( ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
					{
						if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
							continue;

						/* hose the previously computed destruction probability */
						iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest() = SMALLFLOAT; 
					}
				}

				ipLo = 0;
				/* these are the extra Lyman lines for the iso sequences */
				/*for( ipHi=2; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
				/* only update if significant abundance and need to update fine opac */
				if( dense.xIonDense[nelem][ion] > 1e-30 && ( conv.lgFirstSweepThisZone || conv.lgLastSweepThisZone ) )
				{
					for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
					{
						TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
						/* we just want the population of the ground state */
						(*tr).Emis().PopOpc() = iso_sp[ipISO][nelem].st[0].Pop();
						(*(*tr).Lo()).Pop() =
							iso_sp[ipISO][nelem].st[ipLo].Pop();

						/* actually do the work */
						RT_line_one( *tr, true, 0.f,
							GetDopplerWidth(dense.AtomicWeight[nelem]));
					}
				}

			}/* 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( iso_sp[ipH_LIKE][ipHYDROGEN].trans(iso_ctrl.nLyaLevel[ipH_LIKE],0) ) )
	{
		/* 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 */
		t_fe2ovr_la::Inst().atoms_fe2ovr();
		/*fprintf(ioQQQ,"DEBUG fe2 %.2e %.2e\n", hydro.dstfe2lya ,
			hydro.HLineWidth);*/

		/* >>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 */
		iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest() += hydro.dstfe2lya;
	}

	/* is continuum pumping of H Lyman lines included?  yes, but turned off
	 * with atom h-like Lyman pumping off command */
	if( !hydro.lgLymanPumping )
	{
		ipISO = ipH_LIKE;
		nelem = ipHYDROGEN;
		ipLo = 0;
		for( ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
		{
			iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump() = 0.;
		}
	}

	{
		/* following should be set true to print ots contributors for he-like lines*/
		enum {DEBUG_LOC=false};
		if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ )
		{
			/* option to dump a line  */
			DumpLine(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0) );
		}
	}

	/* level 1 lines */
	for( i=1; i <= nLevel1; i++ )
	{
		RT_line_one( TauLines[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
	}

	{
		// 09 mar 28, this transition, between the highest two levels in the
		// very simple version of the model Fe II atom, is troublesome due 
		// to Lya pumping, which can invert the highest two level populations.  
		// The electron scattering escape probability shows noise as a result.  
		// The sim nova_photos will throw an assert if this is removed
		// This model is unphysical since highest two levels are close together
		// and Lya is populating the highest one.  The problem does not occur
		// when the large FeII atom is turned on - that has realistic level
		// energies.  
		// evaluate electron scattering escape probability one time per zone
		static realnum P_elec_esc_ipTFe56;
		static long nSave;
		if( nzone <= 1 || nzone!=nSave )
		{
			P_elec_esc_ipTFe56 = TauLines[ipTFe56].Emis().Pelec_esc();
			nSave = nzone;
		}
		else
			TauLines[ipTFe56].Emis().Pelec_esc() = P_elec_esc_ipTFe56;
	}

	/* external database lines */
	for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
	{
		if( dBaseSpecies[ipSpecies].lgActive )
		{
			realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
			for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin(); 
				  tr != dBaseTrans[ipSpecies].end(); ++tr)
			{	
				int ipHi = (*tr).ipHi();
				if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
					continue;
				if( (*tr).ipCont() > 0 )
				{
					RT_line_one( *tr, true, 0.f, DopplerWidth );
				}
			}
		}
	}

	/* this is a major time sink for this routine - only evaluate on last
	 * sweep when fine opacities are updated since only effect of UTAs is
	 * to pump inner shell lines and add to total opacity */
	if( conv.lgFirstSweepThisZone || conv.lgLastSweepThisZone )
	{
		for( i=0; i < nUTA; i++ )
		{
			/* these are not defined in cooling routines so must be set here */
			UTALines[i].Emis().PopOpc() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
			(*UTALines[i].Lo()).Pop() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
			(*UTALines[i].Hi()).Pop() = 0.;
			RT_line_one( UTALines[i], true,0.f, 
				GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
		}
	}

	/* do satellite lines for iso sequences gt H-like
	 * H-like ions only have one electron, no satellite lines. */
	for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO )
	{
		/* loop over all iso-electronic sequences */
		for( nelem=ipISO; nelem<LIMELM; ++nelem )
		{
			if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
			{
				for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ++ipLo )
				{
					RT_line_one( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]], true,0.f, 
						GetDopplerWidth(dense.AtomicWeight[nelem]) );
				}
			}
		}
	}

	/* the large H2 molecule */
	for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
		(*diatom)->H2_RTMake( );

	return;
}