/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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;

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

	/* 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.;
		
#		ifdef DEBUG_FUN
		fputs( " <->RT_recom_effic()\n", debug_fp );
#		endif
		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( opac.lgTauOutOn )
		{
			/* 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
		 * e2TauAbs is optical depth to illuminated face */
		receff_v = opac.e2TauAbs[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] < (float)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];
		}
	}
	else
	{
		fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n", 
		  rfield.chDffTrns );
		puts( "[Stop in RT_recom_effic]" );
		cdEXIT(EXIT_FAILURE);
	}

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

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

