/* 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_escape do line radiative transfer,
 * evaluates escape and destruction probability */
/* NB no wind distinction - that is done where optical depths are incremented
 * with line opacity - rt_line_one_tauinc and in line width for rad pressure */
/*RT_line_fine_opacity do fine opacities for one line */
/*RT_line_electron_scatter evaluate electron scattering escape probability */
/*RT_line_pumping pumping by external and locally emitted radiation fields */
#include "cddefines.h"
#include "rfield.h"
#include "doppvel.h"
#include "dense.h"
#include "opacity.h"
#include "transition.h"
#include "conv.h"
#include "radius.h"
#include "rt.h"
#include "physconst.h"
#include "cosmology.h"
#include "thirdparty.h"
#include "hydrogenic.h"

/*RT_line_pumping pumping by external and locally emitted radiation fields */
STATIC void RT_line_pumping(
		 /* the em line we will work on  */
		 const TransitionProxy& t ,
		 /* this is option to not include line self shielding across this zone.
		 * this can cause pump to depend on zone thickness, and leads to unstable
		 * feedback in some models with the large H2 molecule, due to Solomon
		 * process depending on zone thickness and level populations. */
		 bool lgShield_this_zone,
		 realnum DopplerWidth)
{
	DEBUG_ENTRY( "RT_line_pumping()" );

	ASSERT( t.ipCont() >= 1 );

	/* pumping by incident and diffuse continua */
	/* option to kill induced processes */
	if( !rfield.lgInducProcess )
	{
		t.Emis().pump() = 0.;
	}
	else if( conv.lgFirstSweepThisZone || t.Emis().iRedisFun() == ipLY_A )
	{
		double dTau;
		double shield_continuum;

		/* continuum shielding into this function */
		shield_continuum = RT_continuum_shield_fcn( t );

		/* continuum upward pumping rate, A gu/gl abs prob occnum
		* the "no induced" command causes continuum pumping to be set to 0 
		* this includes pumping by diffuse continuum */
		t.Emis().pump() = t.Emis().Aul() * (*t.Hi()).g() / (*t.Lo()).g() * shield_continuum *(
			rfield.OccNumbIncidCont[t.ipCont()-1] + rfield.OccNumbContEmitOut[t.ipCont()-1] );
		
		dTau =  (t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth + 
			opac.opacity_abs[t.ipCont()-1])* radius.drad_x_fillfac;
	
		if( lgShield_this_zone && dTau > 1e-3 )
		{
			/* correction for line optical depth across zone */
			t.Emis().pump() *= log(1. + dTau ) / dTau;
		}
		
		/* 
		 * This is an option to account for intrinsic absorption or emission line the lyman 
		 * lines.  This is done without changing the incident coarse continuum.  
		 * Check if continuum pumping of H Lyman lines to be multiplied by scale factor
		 * hydro.xLymanPumpingScaleFactor is set with atom h-like Lyman pumping scale command 
		 * Lya pump rate is always updated - this is simplest way to finese intrinsic absorption
		 * or emission 
		 */
		if ( t.ipLo() == 0 && t.systemIs(iso_sp[ipH_LIKE][ipHYDROGEN].tr) )
		{
			t.Emis().pump() *=  hydro.xLymanPumpingScaleFactor;
		}
		/* NB NB line pumping by local diffuse emission is not included. */
	}

	return;
}

/*RT_line_electron_scatter evaluate electron scattering escape probability */
STATIC void RT_line_electron_scatter(
		 /* the em line we will work on  */
		 const TransitionProxy& t ,
		 realnum DopplerWidth)
{

	DEBUG_ENTRY( "RT_line_electron_scatter()" );

	/* escape following scattering off an electron */
	/* this is turned off with the no scattering escape command */
	if( !rt.lgElecScatEscape )
	{
		t.Emis().Pelec_esc() = 0.;
		return;
	}

	// in the transition structure PopOpc is the pop relative to the ion,
	// but has been converted to a physical population before this routine
	// was called.  No need to convert here */
	double opac_line = (*t.Lo()).Pop() * t.Emis().opacity()/DopplerWidth;

	/* the opacity in electron scattering */
	double opac_electron = dense.eden*6.65e-25;
	/* this is equation 5 of 
	 *>>refer	line	desp	Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752*/
	double opacity_ratio = opac_electron/(opac_electron+opac_line);
	/* keep total probability less than 0.1 */
	t.Emis().Pelec_esc() = (realnum)opacity_ratio * MAX2(0.f,1.f-t.Emis().Pesc()-t.Emis().Pdest());

	return;
}

/*RT_line_escape do line radiative transfer escape and destruction probabilities 
 * this routine sets */
STATIC void RT_line_escape(
	   /* the em line we will work on  */
	   const TransitionProxy &t,
	   /* Stark escape probability to be added to Pesc */
	   realnum pestrk,
	   realnum DopplerWidth,
	   bool lgGoodTau)
{
	int nRedis = -1;

	DEBUG_ENTRY( "RT_line_escape()" );

	/* a runaway maser */
	if( t.Emis().TauIn() < -30. )
	{
		fprintf( ioQQQ, "PROBLEM RT_line_escape called with large negative "
			"optical depth, zone %.2f, setting lgAbort true.\n",
			fnzone );
		DumpLine(t);
		/* return busted true instead of exit here, so that code will
		 * not hang under MPI */
		lgAbort = true;
		return;
	}

	if( cosmology.lgDo )
	{
		/* Sobolev escape */
		if( conv.lgFirstSweepThisZone && lgGoodTau )
		{
			realnum tau_Sobolev =  t.Emis().TauIn();

			if( tau_Sobolev < 1E-5 )
			{
				t.Emis().Pesc() = 1.;
			}
			else
			{
				t.Emis().Pesc() = ( 1.f - exp( -1.f * tau_Sobolev ) )/ tau_Sobolev;
			}

			/* inward escaping fraction */
			t.Emis().FracInwd() = rt.fracin;
		}
		fixit(); // is this correct?
		nRedis = ipDEST_K2;
	}
	/* static solution - which type of line will determine
	 * which redistribution function */
	/* iRedisFun() == 1 - alpha resonance line, partial redistribution,
	 * ipPRD == 1 */
	else if( t.Emis().iRedisFun() == ipPRD )
	{
		/* incomplete redistribution with wings */
		if( conv.lgFirstSweepThisZone && lgGoodTau )
		{
			t.Emis().Pesc() = (realnum)esc_PRD( t.Emis().TauIn(), t.Emis().TauTot(), t.Emis().damp() );

			/* >>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( pestrk > 0.f && t.Emis().Pesc() < 1.f )
				t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );

			/* inward escaping fraction */
			t.Emis().FracInwd() = rt.fracin;
		}
		nRedis = ipDEST_INCOM;
	}

	/* complete redistribution without wings - t.ipLnRedis is ipCRD == -1 */
	else if( t.Emis().iRedisFun() == ipCRD )
	{
		if( conv.lgFirstSweepThisZone && lgGoodTau )
		{
			/* >>chng 01 mar -6, escsub will call any of several esc prob routines,
			* depending of how core is set.  We always want core-only for this option,
			* so call  esca0k2(tau) directly */
			t.Emis().Pesc() = (realnum)esc_CRDcore( t.Emis().TauIn(), t.Emis().TauTot() );

			if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
				t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );

			/* inward escaping fraction */
			t.Emis().FracInwd() = rt.fracin;
		}
		nRedis = ipDEST_K2;
	}

	/* CRD with wings, = 2 */
	else if( t.Emis().iRedisFun() == ipCRDW )
	{
		/* complete redistribution with damping wings */
		if( conv.lgFirstSweepThisZone && lgGoodTau )
		{
			t.Emis().Pesc() = (realnum)esc_CRDwing( t.Emis().TauIn(), t.Emis().TauTot(), t.Emis().damp() );

			if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
				t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );

			/* inward escaping fraction */
			t.Emis().FracInwd() = rt.fracin;
		}
		nRedis = ipDEST_K2;
	}

	/* Lya is special case */
	else if( t.Emis().iRedisFun() == ipLY_A )
	{
		/* incomplete redistribution with wings, for special case of Lya
		* uses fits to Hummer & Kunasz numerical results 
		* this routine is different because escape and dest probs
		* are evaluated together, so no test of lgDoEsc */
		if( lgGoodTau )
		{
			double dest , esin;

			/* this will always evaluate escape prob, no matter what lgDoEsc is.
			* The destruction prob comes back as dest */
			t.Emis().Pesc() = (realnum)RTesc_lya( &esin, &dest, t.Emis().PopOpc(), t, DopplerWidth );

			if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
				t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );

			/* this is current destruction rate */
			t.Emis().Pdest() = (realnum)dest;

			/*  this is fraction of line which is inward escaping */
			t.Emis().FracInwd() = rt.fracin;
		}
	}
	else
	{
		fprintf( ioQQQ, " RT_line_escape called with impossible redistribution function %d\n",
				t.Emis().iRedisFun());
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* only do this if not Lya special case, since dest prob already done */
	if( lgGoodTau && t.Emis().iRedisFun() != ipLY_A && t.Emis().opacity() > 0. )
	{
		t.Emis().Pdest() = (realnum)RT_DestProb(
			/* the abundance of the species */
			t.Emis().PopOpc(),
			/* line center opacity in funny units (needs a vel) */
			t.Emis().opacity(),
			/* array index for line in continuum array,
			 * on f not c scale */
			t.ipCont(),
			/* line width in velocity units */
			DopplerWidth,
			/* escape probability */
			t.Emis().Pesc(),
			/* redistribution function */
			nRedis);
	}

	return;
}

/*RT_line_fine_opacity do fine opacities for one line */
STATIC void RT_line_fine_opacity(
	/* the em line we will work on  */
	const TransitionProxy& t ,
	realnum DopplerWidth)

{
	DEBUG_ENTRY( "RT_line_fine_opacity()" );

	/* this is line center frequency, including bulk motion of gas */
	long int ipLineCenter = t.Emis().ipFine() + rfield.ipFineConVelShift;

	/* define fine opacity fine grid fine mesh */
	/* rfield.lgOpacityFine flag set false with no fine opacities command */
	if( !conv.lgLastSweepThisZone || ipLineCenter < 0 || t.Emis().PopOpc() < SMALLFLOAT ||
		ipLineCenter>rfield.nfine || !rfield.lgOpacityFine )
	{
		return;
	}

	/* number of fine opacity cells corresponding to one doppler width for current
	 * temperature, velocity field, and nuclear mass,
	 * rfield.fine_opac_velocity_width is width per cell, cm/s */
	realnum cells_wide_1x = DopplerWidth/rfield.fine_opac_velocity_width;

	/* line center opacity - type realnum since will add to fine opacity array,
	 * which is realnum */
	realnum opac_line =  (realnum)t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth;

	// this is effective optical depth to this point. Do not do line if 
	// this product is less than SMALLFLOAT
	double dTauEffec = opac_line*radius.depth_x_fillfac;
	if( dTauEffec < SMALLFLOAT )
		return;

	/* core width of optically thick line, do 4x with exponential Doppler core,
	 * must be at least one cell, but profile is symmetric */
	const bool doDamp = dTauEffec*t.Emis().damp()/9. > 0.1;
	long int nCells_core = (long)(cells_wide_1x*4.f + 1.5f);
	/* core is symmetric - make sure both upper and lower bounds
	 * are within continuum energy bin */
	if( ipLineCenter - nCells_core < 1 )
		nCells_core = ipLineCenter - 1;
	if( ipLineCenter + nCells_core > rfield.nfine )
		nCells_core = ipLineCenter -rfield.nfine - 1;

	/* want core to be at least one cell wide */
	nCells_core = MAX2( 1 , nCells_core );

	long int nCells_damp;
	/* include damping wings if optical depth is large */
	if( doDamp )
	{
		// find number of cells to extend the damping wings - cells_wide_1x is one dop width
		// tests with th85orion, stopping at half -h2 point, showed that 0.01 was
		// needed for outer edge, given the definition of dTauEffec
		realnum x = (realnum)sqrt( dTauEffec * t.Emis().damp()*100./SQRTPI ) * cells_wide_1x;
		// test on size of x, which can exceed
		// limits of long in extreme optical depths */
		if( x<LONG_MAX )
		{
			nCells_damp = (long)x;

			if( ipLineCenter-nCells_damp < 1 )
				nCells_damp = ipLineCenter-1;

			if( ipLineCenter+nCells_damp > rfield.nfine )
				nCells_damp = rfield.nfine - ipLineCenter-1;
		}
		else
		{
			/* x was too big, just set to extreme range, which is
			* number of cells to nearest boundary */
			nCells_damp = MIN2( rfield.nfine-ipLineCenter , ipLineCenter )-1;
		}
	}
	else
	{
		nCells_damp = nCells_core;
	}

	static vector<realnum> xprofile, profile;
	xprofile.resize(nCells_damp);
	profile.resize(nCells_damp);

	for( long int i=0; i<nCells_damp; ++i )
	{
		/* distance from line center in units of doppler width */
		xprofile[i] = (realnum) i/cells_wide_1x;
	}
	
	VoigtH(t.Emis().damp(), &xprofile[0], &profile[0], nCells_damp);

	/* line center itself, must not double count here */
	rfield.fine_opac_zone[ipLineCenter] += profile[0]*opac_line;
	for( long int i=1; i<nCells_damp; ++i )
	{
		rfield.fine_opac_zone[ipLineCenter+i] += profile[i]*opac_line;
		rfield.fine_opac_zone[ipLineCenter-i] += profile[i]*opac_line;
	}

	return;
}

/*RT_line_one do rt for emission line structure - calls RT_line_escape or RT_line_wind */
void RT_line_one(
	/* the em line we will work on  */
	const TransitionProxy &t,
	/* this is option to not include line self shielding across this zone.
	 * this can cause pump to depend on zone thickness, and leads to unstable
	 * feedback in some models with the large H2 molecule, due to Solomon
	 * process depending on zone thickness and level populations. */
	bool lgShield_this_zone,
	/* Stark escape probability to be added to Pesc */
	realnum pestrk,
	realnum DopplerWidth )
{
	DEBUG_ENTRY( "RT_line_one()" );

	// do nothing is population and this is not the very first call 
	// skip line transfer if requested with 'no line transfer' command, but never skip Lya
	if( !rfield.lgDoLineTrans && (t.Emis().iRedisFun() != ipLY_A) )
	{
		return;
	}

	/* line damping constant at current temperature  */
	t.Emis().damp() = t.Emis().dampXvel() / DopplerWidth;
	ASSERT( t.Emis().damp() > 0. );

	// do not evaluate if no population 
	if( (*t.Lo()).Pop()<=SMALLFLOAT )
	{
		/* zero population, return after setting everything with side effects */
		t.Emis().Pesc() = 1.f;

		/* inward escaping fraction */
		t.Emis().FracInwd() = 0.5;

		/* pumping rate */
		t.Emis().pump() = 0.;

		/* destruction probability */
		t.Emis().Pdest() = 0.;
		t.Emis().Pelec_esc() = 0.;

		return;
	}

	/* option to keep track of population values during calls,
	 * print out data to make histogram */
	enum {DEBUG_LOC=false};
	if( DEBUG_LOC )
	{
		static long int nTau[100];
		long n;

		if( nzone==0 )
		{
			for(n=0; n<100; ++n )
				nTau[n] = 0;
		}
		if( (*t.Lo()).Pop()<=SMALLFLOAT )
			n = 0;
		else
			n = (long)log10( (*t.Lo()).Pop() )+37;
		n = MIN2( n , 99 );
		n = MAX2( n , 0 );
		++nTau[n];
		if( nzone > 183 )
		{
			for(n=0; n<100; ++n )
				fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] );
			cdEXIT(EXIT_SUCCESS);
		}
	}
	
	// transition is below plasma frequency - photons not emitted
	if( t.EnergyErg() / EN1RYD <= rfield.plsfrq )
	{
		t.Emis().Pesc() = SMALLFLOAT;
		t.Emis().Pdest() = SMALLFLOAT;
		t.Emis().Pelec_esc() = SMALLFLOAT;
		t.Emis().pump() = SMALLFLOAT;
	}
	else
	{

		/* this checks if we have overrun the optical depth scale,
		* in which case the inward optical depth is greater than the
		* previous iteration's total optical depth.
		* We do not reevaluate escape probabilities if the optical depth
		* scale has been overrun due to huge bogus change in solution
		* that would result */
		bool lgGoodTau = lgTauGood( t );

		// the last sweep through this zone is to do the fine opacities
		// the populations are not updated on the last sweep so the 
		// line transfer details also should not be updated
		if( conv.lgLastSweepThisZone )
			RT_line_fine_opacity( t , DopplerWidth );
		else
		{
			RT_line_escape( t, pestrk, DopplerWidth , lgGoodTau);
			RT_line_electron_scatter( t , DopplerWidth );
			RT_line_pumping( t , lgShield_this_zone , DopplerWidth );	
		}
	}

	return;
}