/*RTMakeStat do line radiative transfer for static geometry,
 * evaluates escape and destruction probability, 
 * called by HydroPEsc, RTMake  */
#include "cddefines.h"
#include "called.h"
#include "rfield.h"
#include "doppvel.h"
#include "ionfracs.h"
#include "opacity.h"
#include "rtescprob.h"
#include "converge.h"
#include "radius.h"
#include "rt.h"

void RTMakeStat(EmLine * t , int lgDoEsc )
{
	int lgGoodTau;
	long int nelem;
	char chRedis[5];

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

	/* this is the routine that computes much of the radiative transfer
	 * physics for lines for static case.  */

	assert( t->IonStg < LIMELM+2 );
	assert( t->nelem-1 < LIMELM+2 );

	/* >>chng 01 aug 18, do not evaluate if no population */
	/* molecules are evaluated here, these must be valid */
	/* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
	 * that rates are evaluated at start of a new iteration */
	if( xIonFracs[ t->nelem-1][t->IonStg] == 0.  && !lgDoEsc )
	{
		/* zero population, return after setting everything with side effects */
		t->Pesc = 1.f;

		/* inward escaping fraction */
		t->FracInwd = 0.5;

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

		/* destruction probability */
		t->Pdest = 0.;

#		ifdef DEBUG_FUN
		fputs( " <->RTMakeStat()\n", debug_fp );
#		endif
		return;
	}
	/* check that we don't have a runaway maser */
	else if( t->TauIn < -30. )
	{
		fprintf( ioQQQ, " RTMakeStat called with large negative optical depth.\n" );
		DumpLine(t);
		/* >>>chng 00 apr 23, return busted true instead of exit here, so that code will
		 * not hang under MPI */
		called.lgBusted = TRUE;
		return;
	}

	/* this checks if we have overrun the optical depth scale
	 * possible for line to mase, reason for sec test
	 * >>chng 97 jan 08, on iterations where the first iteration had NO species,
	 * but it does exist on later iterations, added test for existence of ion */
	if( (t->TauTot*0.9 - t->TauIn < 0. && t->TauIn > 0.) && 
		  t->TauTot != opac.taumin )
	{
		/* do not do anything if we have overrun the scale, */
		lgGoodTau = FALSE;
	}
	else
	{
		lgGoodTau = TRUE;
	}

	/* say tau ok if this is first iteration and outward optical depths turned off */
	if( !opac.lgTauOutOn )
	{
		lgGoodTau = TRUE;
	}

	/* atomic number of this element */
	nelem = t->nelem;

	/* thermal plus turbulent velocity for this element */
	assert( (nelem-1) < (LIMELM+2) );

	/* damping constant for the line, save it */
	t->damp = t->damprel / DoppVel.doppler[nelem-1];

	/* static solution - which type of line will determine
	 * which redistribution function */
	/* iRedisFun == 1 - alpha resonance line, partial redistribution */
	if( t->iRedisFun == ipPRD )
	{
		/* incomplete redistribution with wings */
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				t->Pesc = (float)(esc_PRD(t->TauIn,t->TauTot,
				  t->damp));

				/* inward escaping fraction */
				t->FracInwd = rt.fracin;
			}

			/* update pumping even if over run optical depth scale */
			if( rfield.lgInducProcess )
			{
				/* the inward-looking escape probability */
				rt.wayin = (float)esc_PRD_1side(t->TauCon,t->damp);

				/* continuum upward pumping rate, A gu/gl abs prob occnum
				 * the "no induced" command causes conpmp to set 0 */
				t->pump = t->Aul * t->gHi / t->gLo * rt.wayin *
				  rfield.OccNumbIncidCont[t->ipCont-1];

				/* >>chng 99 dec 02, add correction for line optical depth across zone */
				/* correction for attenuation within line across zone */
				/* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow, 
				 * this caused an overestimate of pumping rates */
				/* >>chng 00 dec 19, nova.in oscillated since correction factor was
				 * very large across zones that were already very optically thick, dr changed,
				 * causing feedback between ionization and thickness
				 * the intention of this correction was to only apply when on the 
				 * verge of optically thick,
				 * change lower bound on product t->dTau * radius.dReff to 1e-3 to avoid
				 * constantly getting unity, also check on relative change optical depth/optical depth*/
				/*if( t->dTau * radius.dReff > 1e-8 )*/
				if( t->dTau * radius.dReff > 1e-3 && t->dTau * radius.dReff < 10. )
				{
					/* during very first search dTau is -1e38*/
					t->pump *= (float)(log(1. + t->dTau * radius.dReff ) / (t->dTau * radius.dReff) );
				}

				/* >>chng 00 Oct 03, add pumping by diffuse continuum,
				 * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in  parsedont.c
				 * with no diffuse line pumping command */
				if( t->dTau*rfield.DiffPumpOn > SMALLFLOAT )
				{
					/* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so
					 * that smallest and largest number are multiplied first to prevent overflow, PvH */
					float scale;
					if( iteration > 1 )
					{
						/* inward looking depth, t->dTau is opacity in line, cm^-1,
						 * so multiplying by this scale is same as dividing by opacity in 
						 * source function */
						scale = MIN2( 1.f/ t->dTau , (float)radius.depth );
						t->pump += t->Aul * t->gHi / t->gLo * (0.5f*scale*
							rfield.OccNumbDiffCont[t->ipCont-1]);
						/* outward looking depth */
						scale = MIN2( 1.f/ t->dTau , (float)radius.Depth2Go );
						t->pump += t->Aul * t->gHi / t->gLo * (0.5f*MAX2(0.f,scale)*
							rfield.OccNumbDiffCont[t->ipCont-1]);
					}
					else 
					{
						/* first iteration, only inward looking depths really known */
						scale = 1.f/ t->dTau;
						scale = MIN2( scale , (float)radius.depth );
						t->pump += t->Aul * t->gHi / t->gLo * (scale*
							rfield.OccNumbDiffCont[t->ipCont-1]);
					}
				}
			}
			else
			{
				t->pump = 0.;
			}
		}

		strcpy( chRedis , "INCO" );
	}

	/* complete redistribution without wings - t(ipLnRedis) is -1 */
	else if( t->iRedisFun == ipCRD )
	{

		if( lgDoEsc )
		{
			if( 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->Pesc = (float)esc_CRDcore(t->TauIn,t->TauTot);

				/* inward escaping fraction */
				t->FracInwd = rt.fracin;
			}

			/* continuum pumping rate
			 * the "no induced" command causes conpmp to set 0 */
			if( rfield.lgInducProcess )
			{
				/* wayin = esc_CRDwing_1side( t(ipLnTauCon),damp )
				 * this is to make it more like C90 was
				 * complete redistribution with doppler core only */
				rt.wayin = (float)esca0k2(t->TauCon);

				t->pump = t->Aul * t->gHi / t->gLo * rt.wayin*
				  rfield.OccNumbIncidCont[t->ipCont-1];
				/* >>chng 99 dec 02, add correction for line optical depth across zone,
				 * dTau is actually differential opacity */
				/* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow, 
				 * this caused an overestimate of pumping rates */
				/* >>chng 00 dec 19, nova.in oscillated since correction factor was
				 * very large across zones that were already very optically thick, dr changed,
				 * causing feedback between ionization and thickness
				 * the intention of this correction was to only apply when on the 
				 * verge of optically thick,
				 * change lower bound on product t->dTau * radius.dReff to 1e-3 to avoid
				 * constantly getting unity, also check on relative change optical depth/optical depth*/
				/*if( t->dTau * radius.dReff > 1e-8 )*/
				if( t->dTau * radius.dReff > 1e-3 && t->dTau * radius.dReff < 10. )
				{
					/* during very first search dTau is -1e38*/
					t->pump *= (float)(log(1. + t->dTau * radius.dReff ) / (t->dTau * radius.dReff) );
				}

				/* >>chng 00 oct 03, add pumping by diffuse continuum,
				 * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in parsedont.ce
				 * with no diffuse line pumping command */
				if( t->dTau*rfield.DiffPumpOn> SMALLFLOAT )
				{
					/* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so
					 * that smallest and largest number are multiplied first to prevent overflow, PvH */
					float scale;
					if( iteration > 1 )
					{
						/* inward looking depth */
						scale = MIN2( 1.f/ t->dTau , (float)radius.depth );
						t->pump += t->Aul * t->gHi / t->gLo * (0.5f*scale*
							rfield.OccNumbDiffCont[t->ipCont-1]);

						/* outward looking depth */
						scale = MIN2( 1.f/ t->dTau , (float)radius.Depth2Go );
						t->pump += t->Aul * t->gHi / t->gLo * (0.5f*MAX2(0.f,scale)*
							rfield.OccNumbDiffCont[t->ipCont-1]);
					}
					else 
					{
						/* first iteration, only inward looking depths really known */
						scale = 1.f/ t->dTau;
						scale = MIN2( scale , (float)radius.depth );
						t->pump += t->Aul * t->gHi / t->gLo * (scale*
							rfield.OccNumbDiffCont[t->ipCont-1]);
					}
				}
			}
			else
			{
				t->pump = 0.;
			}
		}

		strcpy( chRedis , "  K2" );
	}

	/* chng 01 feb 27, add CRD with wings, = 2 */
	else if( t->iRedisFun == ipCRDW )
	{
		/* complete redistribution with damping wings */
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				t->Pesc = (float)esc_CRDwing(t->TauIn,t->TauTot,t->damp);

				/* inward escaping fraction */
				t->FracInwd = rt.fracin;
			}

			/* continuum pumping rate
			 * the "no induced" command causes conpmp to set 0 */
			if( rfield.lgInducProcess )
			{
				/* wayin = esc_CRDwing_1side( t(ipLnTauCon),damp )
				 * this is to make it more like C90 was
				 * complete redistribution with doppler core only */
				/* >>chng 01 may 29, had been esc_PRD_1side, as per above comment,
				 * but this should include damping wings.  */
				/*rt.wayin = (float)esc_PRD_1side(t->TauCon,damp);*/
				rt.wayin = (float)esc_CRDwing_1side(t->TauCon,t->damp);

				t->pump = t->Aul * t->gHi / t->gLo * rt.wayin*
				  rfield.OccNumbIncidCont[t->ipCont-1];
				/* >>chng 99 dec 02, add correction for line optical depth across zone,
				 * dTau is actually differential opacity */
				/* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow, 
				 * this caused an overestimate of pumping rates */
				/* >>chng 00 dec 19, nova.in oscillated since correction factor was
				 * very large across zones that were already very optically thick, dr changed,
				 * causing feedback between ionization and thickness
				 * the intention of this correction was to only apply when on the 
				 * verge of optically thick,
				 * change lower bound on product t->dTau * radius.dReff to 1e-3 to avoid
				 * constantly getting unity, also check on relative change optical depth/optical depth*/
				/*if( t->dTau * radius.dReff > 1e-8 )*/
				if( t->dTau * radius.dReff > 1e-3 && t->dTau * radius.dReff < 10. )
				{
					/* during very first search dTau is -1e38*/
					t->pump *= (float)(log(1. + t->dTau * radius.dReff ) / (t->dTau * radius.dReff) );
				}

				/* >>chng 00 oct 03, add pumping by diffuse continuum,
				 * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in parsedont.ce
				 * with no diffuse line pumping command */
				if( t->dTau*rfield.DiffPumpOn> SMALLFLOAT )
				{
					/* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so
					 * that smallest and largest number are multiplied first to prevent overflow, PvH */
					float scale;
					if( iteration > 1 )
					{
						/* inward looking depth */
						scale = MIN2( 1.f/ t->dTau , (float)radius.depth );
						t->pump += t->Aul * t->gHi / t->gLo * (0.5f*scale*
							rfield.OccNumbDiffCont[t->ipCont-1]);
						/* outward looking depth */
						scale = MIN2( 1.f/ t->dTau , (float)radius.Depth2Go );
						t->pump += t->Aul * t->gHi / t->gLo * (0.5f*MAX2(0.f,scale)*
							rfield.OccNumbDiffCont[t->ipCont-1]);
					}
					else 
					{
						/* first iteration, only inward looking depths really known */
						scale = 1.f/ t->dTau;
						scale = MIN2( scale , (float)radius.depth );
						t->pump += t->Aul * t->gHi / t->gLo * (scale*
							rfield.OccNumbDiffCont[t->ipCont-1]);
					}
				}
			}
			else
			{
				t->pump = 0.;
			}
		}

		strcpy( chRedis , "  K2" );
	}
	else
	{
		fprintf( ioQQQ, " RTMakeStat called with redistribution function zero\n" );
		ShowMe();
		puts( "[Stop in RTMakeStat]" );
		cdEXIT(1);
	}

	/* >>>chng 99 dec 18, did not have the test for lgGoodTau, and so clobbered
	 * dest prob when overrun */
	if( lgGoodTau )
	{
		float DestSave = t->Pdest;
		t->Pdest = (float)(RT_DestProb(
			/* the abundance of the species */
			t->PopOpc,
			/* line center opacity in funny units (needs a vel) */
			t->opacity,
			/* array index for line in continuum array,
			 * on f not c scale */
			t->ipCont,
			/* line width in velocity units */
			DoppVel.doppler[nelem-1],
			/* escape probability */
			t->Pesc,
			/* redistribution function */
			chRedis));

		if( !conv.lgSearch )
			t->Pdest = 0.75f*t->Pdest + 0.25f * DestSave;
	}

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