/*receff generate escape probability function for continua, */
#include "cddefines.h"
#include <string.h>
#include "physconst.h"
#include "e2tau.h"
#include "rfield.h"
#include "otsmin.h"
#include "phycon.h"
#include "touton.h"
#include "dtrans.h"
#include "opac.h"
#include "facexp.h"
#include "rtescprob.h"
#include "receff.h"

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

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

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

	/* dtrans.chDffTrns = "OU2" by default */
	if( ip > rfield.nflux )
	{
		receff_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->receff()\n", debug_fp );
#		endif
		return( receff_v );
	}

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

	/* inward optical depth and escape prob */
	if( strcmp(dtrans.chDffTrns,"OTS") == 0 )
	{
		tin = opac.tauabs[0][ip-1];
		if( tin < 5. )
		{
			escin = esccon(tin,hnukt);
		}
		else
		{
			escin = 1e-4;
		}

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

		if( touton.lgTauOutOn )
		{
			if( tout > 0. )
			{
				if( tout < 5. )
				{
					escout = esccon(tout,hnukt);
				}
				else
				{
					escout = 1e-4;
				}
				receff_v = 0.5*(escin + escout);
			}
			else
			{
				/* this logic added april 93 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(dtrans.chDffTrns,"OU1") == 0 )
	{
		receff_v = facexp.ExpZone[ip+1];
	}
	else if( strcmp(dtrans.chDffTrns,"OU2") == 0 )
	{
		/* this is the default rt method, asset in zero
		 * e2TauAbs is optical depth to illuminated face */
		receff_v = e2tau.e2TauAbs[ip+1];
	}
	else if( strcmp(dtrans.chDffTrns,"OU3") == 0 )
	{
		receff_v = 1.;
	}
	else if( strcmp(dtrans.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]) < 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", 
		  dtrans.chDffTrns );
		puts( "[Stop in receff]" );
		exit(1);
	}

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

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