/*RTMakeStat do line radiative transfer for static geometry,
 * evaluates escape and destruction probability, 
 * called by HydroPEsc, RTMake  */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "twoway.h"
#include "showme.h"
#include "rfield.h"
#include "doppvel.h"
#include "fluors.h"
#include "taumin.h"
#include "rtescprob.h"
#include "radius.h"
#include "dumpline.h"
#include "touton.h"
#include "eovrlp.h"
#include "rt.h"

void RTMakeStat(EmLine * t , int lgDoEsc )
{
	int lgGoodTau;
	long int nelem;
	float damp, 
	  velocity;
	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.  */

	/* iplntauin is 1, so -1 is 0 */
	if( t->TauIn < -30. )
	{
		fprintf( ioQQQ, " RTMakeStat called with large negative optical depth.\n" );
		DumpLine(t);
		ShowMe();
		puts( "[Stop in RTMakeStat]" );
		exit(1);
	}
	/*end sanity check */

	/* this checks if we have overrun the optical depth scale
	 * possible for line to mase, reason for sec test
	 * >>chng 97 jan 8, 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 != tauminCom.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( !touton.lgTauOutOn )
	{
		lgGoodTau = TRUE;
	}

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

	/* thermal plus turbulent velocity for this element */
	velocity = DoppVel.doppler[nelem-1];

	/* damping constant for the line, save it */
	damp = t->damprel/velocity;
	t->damp = damp;

	/* static solution - which redistribution function? */
	if( t->iRedisFun > 0 )
	{
		/* incomplete redistribution  */
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				t->Pesc = (float)(escgrd(t->TauIn,t->TauTot,
				  damp));

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

			/* update pumping even if over run optical depth scale */
			if( fluors.lgFluor )
			{
				/* the inward-looking escape probability */
				twoway.wayin = (float)escinc(t->TauCon,damp);

				/* continuum 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 * twoway.wayin *
				  rfield.occnum[t->ipCont-1];

				/* >>chng 99 Dec 2, add correction for line optical depth across zone */
				/* correction for attenuation within line across zone */
				if( t->dTau * radius.dReff > 0. )
				{
					/* during very first search dTau is -1e38*/
					t->pump *= (float)(log(1. + t->dTau * radius.dReff ) / (t->dTau * radius.dReff) );
				}
			}
			else
			{
				t->pump = 0.;
			}
		}

		strcpy( chRedis , "INCO" );
	}

	else if( t->iRedisFun < 0 )
	{

		/* complete redistribution - t(ipLnRedis) is negative */
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				t->Pesc = (float)escsub(t->TauIn,t->TauTot,
				  damp);

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

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

				t->pump = t->Aul * t->gHi / t->gLo * twoway.wayin*
				  rfield.occnum[t->ipCont-1];
				/* >>chng 99 Dec 2, add correction for line optical depth across zone */
				if( t->dTau * radius.dReff > 0. )
				{
					/* during very first search dTau is -1e38*/
					t->pump *= (float)(log(1. + t->dTau * radius.dReff ) / (t->dTau * radius.dReff) );
				}
			}
			else
			{
				t->pump = 0.;
			}
		}

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

	/* >>>chng 99 dec 18, did not have the test for lgGoodTau, and so clobbered
	 * dest prob when overrun */
	if( lgGoodTau )
	{
		t->Pdest = (float)(eovrlp(
			/* the abundance of the species */
			t->PopOpc,
			/* line center opacity in funny units (needs a vel) */
			t->opacity,
			/* pointer to line in continuum array */
			t->ipCont,
			/* line width in velocity units */
			velocity,
			/* escape probability */
			t->Pesc,
			/* redistribution function */
			chRedis));
	}

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