/* 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_recom_effic generate escape probability function for continua, */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "phycon.h"
#include "opacity.h"
#include "rt.h"

double RT_recom_effic(long int ip)
{
	long int i;
	double dEner, 
	  denom, 
	  escin, 
	  escout, 
	  hnukt, 
	  receff_v, 
	  sum, 
	  tin, 
	  tout;

	DEBUG_ENTRY( "RT_recom_effic()" );

	/* escape probability function for continua,
	 * formally correct for photoelectric absorption only */

	ASSERT( ip > 0 && ip <= rfield.nupper );

	if( ip > rfield.nflux )
	{
		/* >>chng 01 dec 18, return had been zero, but this did not
		 * work for case where gas much hotter than continuum, as in a
		 * coronal plasma.  change to return of unity */
		receff_v = 1.;
		return( receff_v );
	}

	/* bug in following statement unocvered June 93 S. Schaefer */
	hnukt = TE1RYD*rfield.anu[ip-1]/phycon.te;

	/* rfield.chDffTrns = "OU2" by default */
	/* inward optical depth and escape prob */
	if( strcmp(rfield.chDffTrns,"OSS") == 0 )
	{
		/* which side of Lyman limit is this? */
		if( rfield.anu[ip] > 0.99 )
		{
			/* this is a simple OTS, turned on with DIFFUSE OTS SIMPLE */
			receff_v = SMALLFLOAT;
		}
		else
		{
			receff_v = 1.;
		}
	}
	else if( strcmp(rfield.chDffTrns,"OTS") == 0 )
	{
		tin = opac.TauAbsGeo[0][ip-1];
		if( tin < 5. )
		{
			escin = esccon(tin,hnukt);
		}
		else
		{
			escin = 1e-4;
		}

		/* outward optical depth */
		tout = opac.TauAbsGeo[1][ip-1] - tin;

		if( iteration > 1 )
		{
			/* check whether we have overrun the optical depth scale */
			if( tout > 0. )
			{
				/* good optical depths in both directions, take mean */
				if( tout < 5. )
				{
					escout = esccon(tout,hnukt);
				}
				else
				{
					escout = 1e-4;
				}
				receff_v = 0.5*(escin + escout);
			}
			else
			{
				/* >>chng 91 apr add logic to prevent big change in
				 * esc prob, resulting in terminal oscillations, when optical
				 * depth scale overrun
				 * tau was negative, use 5% of inward optical depth */
				escout = esccon(tin*0.05,hnukt);
				receff_v = 0.5*(escin + escout);
			}
		}
		else
		{
			receff_v = escin;
		}
	}
	else if( strcmp(rfield.chDffTrns,"OU1") == 0 )
	{
		receff_v = opac.ExpZone[ip+1];
	}
	else if( strcmp(rfield.chDffTrns,"OU2") == 0 )
	{
		/* this is the default rt method, as set in zero
		 * E2TauAbsFace is optical depth to illuminated face */
		receff_v = opac.E2TauAbsFace[ip+1];
	}
	else if( strcmp(rfield.chDffTrns,"OU3") == 0 )
	{
		receff_v = 1.;
	}
	else if( strcmp(rfield.chDffTrns,"OU4") == 0 )
	{
		/* this cannot happen, was the former outward treat
		 * optical depth for this zone */
		if( rfield.ContBoltz[ip-1] > 0. )
		{
			i = ip;
			dEner = phycon.te/TE1RYD*0.5;
			sum = 0.;
			denom = 0.;
			while( rfield.ContBoltz[i-1] > 0. &&
			   rfield.anu[i-1]-rfield.anu[ip-1] < (realnum)dEner &&
			   i <= rfield.nflux )
			{
				sum += rfield.ContBoltz[i-1]*opac.tmn[i-1];
				denom += rfield.ContBoltz[i-1];
				i += 1;
			}
			receff_v = sum/denom;
		}
		else
		{
			receff_v = opac.tmn[ip-1];
		}
	}
#if	0	
	else if( strcmp(rfield.chDffTrns,"SOB") == 0 )
	{
		long int ipRecombEdgeFine = rfield.ipnt_coarse_2_fine[ip];
		double OpacityEffective, EffectiveThickness;
		realnum tau;

		/* find line center opacity - use fine opacity if array indices are OK */
		if( ipRecombEdgeFine>=0 && ipRecombEdgeFine<rfield.nfine && rfield.lgOpacityFine )
		{
			/* use fine opacities fine grid fine mesh to get optical depth 
			 * to continuum source */
			/* total line center optical depth, all overlapping lines included */
			OpacityEffective = rfield.fine_opac_zone[ipRecombEdgeFine];
		}
		else
		{
			OpacityEffective =  opac.opacity_abs[ip];
		}

		if( cosmology.lgDo )
		{
			/* dv/dr (s-1), equal to dv/dt / v */
			/* in this case, dv/dr is just the Hubble factor */
			realnum dvdr = GetHubbleFactor(cosmology.redshift_current);
			/* here, we assume that recombining electrons have energy kT
			 * and that photons resulting from recombination must be redshifted
			 * by kT before they are too weak to photoionize the same lower level again */
			realnum width_to_shift = phycon.te_ryd/(rfield.anu[ip]+phycon.te_ryd) * SPEEDLIGHT;
			EffectiveThickness = width_to_shift / dvdr;
			tau = (realnum)(OpacityEffective * EffectiveThickness);
		}
		else
			tau = opac.taumin;

		tau = MAX2((double)opac.taumin,tau);

		ASSERT( tau >= 0. );

		if( tau < 1e-5 )
			receff_v = (1. - tau/2.);
		else
			receff_v = (1. - sexp( tau ) )/ tau;
		ASSERT( receff_v >= 0.f );
		ASSERT( receff_v <= 1.f );
	}
#endif	
	else
	{
		fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n", 
		  rfield.chDffTrns );
		cdEXIT(EXIT_FAILURE);
	}

	receff_v = MAX2((double)opac.otsmin,receff_v);
	/* can get epsilon above unity on cray */
	receff_v = MIN2(1.,receff_v);
	return( receff_v );
}