/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */
/*esc_PRD_1side fundamental escape probability radiative transfer routine for incomplete redistribution */
/*RTesc_lya escape prob for hydrogen atom Lya, using Hummer and Kunasz results,
 * called by hydropesc */
/*esc_PRD escape probability radiative transfer for incomplete redistribution */
/*esc_CRDwing escape probability for CRD with wings */
/*esc_CRDcore escape probability for CRD with no wings */
/*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
/*RT_DestProb returns line destruction probability due to continuum opacity */
/*RT_LineWidth determine half width of any line with known optical depths */
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#define	SCALE	2.
#include "dense.h"
#include "converge.h"
#include "rfield.h"
#include "opacity.h"
#include "lines_service.h"
#include "taulines.h"
#include "doppvel.h"
#include "pressure.h"
#include "wind.h"
#include "rt.h"
/*
 *variables used to pass continuum optical depth and temp to
 *routine that integrates over continuum
 */
static double chnukt_ContTkt, chnukt_ctau;

/*escmase escape probability for negative (masing) optical depths,*/
static double escmase(double tau);
/*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */
static void RTesc_lya_1side(double taume, 
  double beta, 
  float *esc, 
  float *dest,
  /* position of line in frequency array on c scale */
  long ipLine );

double esc_PRD_1side(double tau, 
  double a)
{
	double atau, 
	  b, 
	  escinc_v;

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


	/* this is one of the three fundamental escape probability routines
	 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
	 * it computes esc prob for incomplete redistribution
	 * */
#	if 0
	if( strcmp(rt.chEscFunSubord,"SIMP") == 0 )
	{
		/* this set with "escape simple" command, used for debugging */
		escinc_v = 1./(1. + tau);
		
#		ifdef DEBUG_FUN
		fputs( " <->esc_PRD_1side()\n", debug_fp );
#		endif
		return( escinc_v );
	}
#	endif

	if( tau < 0. )
	{
		/* line mased */
		escinc_v = escmase(tau);
	}
	else if( tau < 10. )
	{
		/* linear part of doppler core */
		escinc_v = 1./(1. + 1.6*tau);
	}
	else
	{
		/* first find coeficient b(tau) */
		atau = a*tau;
		if( atau > 1. )
		{
			b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + atau);
		}
		else
		{
			b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + 1./sqrt(atau));
		}
		b = MIN2(6.,b);

		escinc_v = 1./(1. + b*tau);
	}

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

/*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */
double esc_CRDwing_1side(double tau, 
  double a )
{
	double esccom_v;

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

	/* this is one of the three fundamental escape probability routines
	 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
	 * it computes esc prob for complete redistribution with wings
	 * computes escape prob for complete redistribution in one direction
	 * */

	/* this is the only case that this routine computes,
	 * and is the usual case for subordinate lines, 
	 * complete redistribution with damping wings */
	esccom_v = esca0k2(tau);
	if( tau > 1e3 )
	{
		esccom_v += 0.333*sqrt(a/(SQRTPI*tau));
	}

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

/*RTesc_lya escape prob for hydrogen atom Lya, using 
 >>refer	La	escp	Hummer, D.G., & Kunasz, P.B., 1980, ApJ, 236, 609
 * called by hydropesc, return value is escape probability */
double RTesc_lya(
	/* the inward escape probability */
	float *esin, 
	/* the destruction probility */
	float *dest, 
	/* abundance of the species */
	double abund, 
	/* element number,  0 for H */
	long int nelem)
{
	double beta, 
	  conopc, 
	  escla_v;
	 float dstin, 
	  dstout;

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

	/* 
	 * this is one of the three fundamental escape probability functions
	 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
	 * evaluate esc prob for LA
	 * optical depth in outer direction always defined
	 */

	/* check charge */
	ASSERT( nelem >= 0 && nelem < LIMELM );

	if( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot - 
	  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn < 0. )
	{
		/* this is the case if we overrun the optical depth scale
		 * just leave things as they are */
		escla_v = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pesc;
		rt.fracin = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].FracInwd;
		*esin = rt.fracin;
		*dest = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest;
		
#		ifdef DEBUG_FUN
		fputs( " <->RTesc_lya()\n", debug_fp );
#		endif
		return( escla_v );
	}

	/* incomplete redistribution */
	conopc = opac.opacity_abs[EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1];
	if( abund > 0. )
	{
		/* the continuous opacity is positive, we have a valid soln */
		beta = conopc/(abund/SQRTPI*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity/
		  DoppVel.doppler[nelem] + conopc);
	}
	else
	{
		/* abundance is zero, set miniumum dest prob */
		beta = 1e-10;
	}

	/* find rt.wayin, the escape prob in inward direction */
	RTesc_lya_1side(
	  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn,
	  beta,
	  &rt.wayin,
	  &dstin , 
	  /* position of line in energy array on C scale */
	  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1);

	ASSERT( (rt.wayin <= 1.) && (rt.wayin >= 0.) && (dstin <= 1.) && (dstin >= 0.) );

	/* find rt.wayout, the escape prob in outward direction */
	RTesc_lya_1side(MAX2(EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot/100.,
	  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot-EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn),
	  beta,
	  &rt.wayout,
	  &dstout, 
	  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1);

	ASSERT( (rt.wayout <= 1.) && (rt.wayout >= 0.) && (dstout <= 1.) && (dstout >= 0.) );

	/* esc prob is mean of in and out */
	escla_v = (rt.wayin + rt.wayout)/2.;
	/* the inward escaping part of the line */
	*esin = rt.wayin;

	/* dest prob is mean of in and out */
	*dest = (dstin + dstout)/2.f;
	/* >>chng 02 oct 02, sum of escape and dest prob must be less then unity,
	 * for very thin models this forces dest prob to go to zero, 
	 * rather than the value of DEST0, caught by Jon Slavin */
	*dest = (float)MIN2( *dest , 1.-escla_v );
	/* but dest prob can't be negative */
	*dest = (float)MAX2(0., *dest );

	/* fraction of line emitted in inward direction */
	rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
	ASSERT( escla_v >=0. && *dest>=0. && *esin>=0. );

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

/*esc_PRD escape probability radiative transfer for incomplete redistribution */
double esc_PRD(double tau, 
  double tau_out, 
  double damp )
{
	double escgrd_v, 
	  tt;

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

	/* find escape prob for incomp redis, average of two 1-sided probs*/

	if( opac.lgTauOutOn )
	{
		/*  outward optical depth if defined */
		tt = tau_out - tau;
		/*  help convergence by not letting tau go to zero at back edge of
		 *  when there was a bad guess for the total optical depth
		 *  note that this test is seldom hit since RTMakeStat does check
		 *  for overrun */
		if( tt < 0. )
		{
			tt = tau/SCALE;
		}

		rt.wayin = (float)esc_PRD_1side(tau,damp);
		rt.wayout = (float)esc_PRD_1side(tt,damp);
		rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
		escgrd_v = 0.5*(rt.wayin + rt.wayout);
	}
	else
	{
		/*  outward optical depth not defined, dont estimate fraction out */
		rt.fracin = 0.;
		rt.wayout = 1.;
		escgrd_v = esc_PRD_1side(tau,damp);
		rt.wayin = (float)escgrd_v;
	}

	ASSERT( escgrd_v > 0. );

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

/*esc_CRDwing escape probability radiative transfer for CRDS in core only */
double esc_CRDwing(double tau_in, 
  double tau_out, 
  double damp)
{
	double escgrd_v, 
	  tt;

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

	/* find escape prob for CRD with damping wings, average of two 1-sided probs*/

	/* crd with wings */
	if( opac.lgTauOutOn )
	{
		/*  outward optical depth if defined */
		/* >>chng 03 jun 07, add test for masers here */
		if( tau_out <0 || tau_in < 0. )
		{
			/* we have a maser, use smallest optical depth to damp it out */
			tt = MIN2( tau_out , tau_in );
			tau_in = tt;
		}
		else
		{
			tt = tau_out - tau_in;
			/*  help convergence by not letting tau_in go to zero at back edge of
			 *  when there was a bad guess for the total optical depth
			 *  note that this test is seldom hit since RTMakeStat does check
			 *  for overrun */
			if( tt < 0. )
			{
				tt = tau_in/SCALE;
			}
		}

		rt.wayin = (float)esc_CRDwing_1side(tau_in,damp);
		rt.wayout = (float)esc_CRDwing_1side(tt,damp);
		rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
		escgrd_v = 0.5*(rt.wayin + rt.wayout);
	}
	else
	{
		/*  outward optical depth not defined, dont estimate fraction out */
		rt.fracin = 0.;
		rt.wayout = 1.;
		escgrd_v = esc_CRDwing_1side(tau_in,damp);
		rt.wayin = (float)escgrd_v;
	}

	ASSERT( escgrd_v > 0. );

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

/*esc_CRDwing escape probability radiative transfer for incomplete redistribution */
double esc_CRDcore(double tau_in, 
  double tau_out)
{
	double escgrd_v, 
	  tt;

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

	/* find escape prob for CRD with damping wings, average of two 1-sided probs*/

	/* crd with wings */
	if( opac.lgTauOutOn )
	{
		/*  outward optical depth if defined */
		/* >>chng 03 jun 07, add test for masers here */
		if( tau_out <0 || tau_in < 0. )
		{
			/* we have a maser, use smallest optical depth to damp it out */
			tt = MIN2( tau_out , tau_in );
			tau_in = tt;
		}
		else
		{
			tt = tau_out - tau_in;
			/*  help convergence by not letting tau_in go to zero at back edge of
			 *  when there was a bad guess for the total optical depth
			 *  note that this test is seldom hit since RTMakeStat does check
			 *  for overrun */
			if( tt < 0. )
			{
				tt = tau_in/SCALE;
			}
		}

		rt.wayin = (float)esca0k2(tau_in);
		rt.wayout = (float)esca0k2(tt);
		rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
		escgrd_v = 0.5*(rt.wayin + rt.wayout);
	}
	else
	{
		/*  outward optical depth not defined, dont estimate fraction out */
		rt.fracin = 0.;
		rt.wayout = 1.;
		escgrd_v = esca0k2(tau_in);
		rt.wayin = (float)escgrd_v;
	}

	ASSERT( escgrd_v > 0. );

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

/*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
double esca0k2(double taume)
{
	double arg, 
	  esca0k2_v, 
	  suma, 
	  sumb, 
	  sumc, 
	  sumd, 
	  tau;
	static double a[5]={1.00,-0.1117897,-0.1249099917,-9.136358767e-3,
	  -3.370280896e-4};
	static double b[6]={1.00,0.1566124168,9.013261660e-3,1.908481163e-4,
	  -1.547417750e-7,-6.657439727e-9};
	static double c[5]={1.000,19.15049608,100.7986843,129.5307533,-31.43372468};
	static double d[6]={1.00,19.68910391,110.2576321,169.4911399,-16.69969409,
	  -36.664480000};

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

	/* compute Hummer's K2 escape probability function for a=0
	 * using approx from 
	 * >>refer	line	escp	Hummer, D.G., xxxx, JQRST, 26, 187.
	 *
	 * convert to David's opacity */
	tau = taume*SQRTPI;

	if( tau < 0. )
	{
		/* the line mased */
		esca0k2_v = escmase(taume);

	}
	else if( tau < 0.01 )
	{
		esca0k2_v = 1. - 2.*tau;

	}
	else if( tau <= 11. )
	{
		suma = a[0] + tau*(a[1] + tau*(a[2] + tau*(a[3] + a[4]*tau)));
		sumb = b[0] + tau*(b[1] + tau*(b[2] + tau*(b[3] + tau*(b[4] + 
		  b[5]*tau))));
		esca0k2_v = tau/2.5066283*log(tau/SQRTPI) + suma/sumb;

	}
	else
	{
		/* large optical depth limit */
		arg = 1./log(tau/SQRTPI);
		sumc = c[0] + arg*(c[1] + arg*(c[2] + arg*(c[3] + c[4]*arg)));
		sumd = d[0] + arg*(d[1] + arg*(d[2] + arg*(d[3] + arg*(d[4] + 
		  d[5]*arg))));
		esca0k2_v = (sumc/sumd)/(2.*tau*sqrt(log(tau/SQRTPI)));
	}

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

/*escmase escape probability for negative (masing) optical depths */
void FindNeg( void )
{
	long int i;

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

	/* do the level 1 lines */
	for( i=1; i <= nLevel1; i++ )
	{
		/* check if a line was a strong maser */
		if( TauLines[i].TauIn < -1. )
			DumpLine(&TauLines[i]);
	}

	/* now do the level 2 lines */
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			/* check if a line was a strong maser */
			if( TauLine2[i].TauIn < -1. )
				DumpLine(&TauLine2[i]);
		}
	}

	/* now do the hyperfine structure lines */
	for( i=0; i < nHFLines; i++ )
	{
		/* check if a line was a strong maser */
		if( HFLines[i].TauIn < -1. )
			DumpLine(&HFLines[i]);
	}

	/* now do the co carbon monoxoide lines */
	for( i=0; i < nCORotate; i++ )
	{
		/* check if a line was a strong maser */
		if( C12O16Rotate[i].TauIn < -1. )
			DumpLine(&C12O16Rotate[i]);
		/* check if a line was a strong maser */
		if( C13O16Rotate[i].TauIn < -1. )
			DumpLine(&C13O16Rotate[i]);
	}

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

static double escmase(double tau)
{
	double escmase_v;

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

	/* this is the only routine that computes maser escape proba */
	ASSERT( tau <= 0. );

	if( tau > -0.1 )
	{
		escmase_v = 1. - tau*(0.5 + tau/6.);
	}
	else if( tau > -30. )
	{
		escmase_v = (1. - exp(-tau))/tau;
	}
	else
	{
		fprintf( ioQQQ, " escmase called with 2big tau%10.2e\n", 
		  tau  );
		fprintf( ioQQQ, " This is zone number%4ld\n", nzone );
		FindNeg();
		ShowMe();
		puts( "[Stop in escmase]" );
		cdEXIT(EXIT_FAILURE);
	}

	ASSERT( escmase_v >= 1. );

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

/*escConE2 one of the forms of the continuum escape probability */
/*cone2 generate e2 function needed for continuum transport */
double cone2(double t);

double escConE2(double x)
{
	double conesc_v;

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

	conesc_v = exp(-chnukt_ContTkt*(x-1.))/
		x*cone2(chnukt_ctau/x/x/x);

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

/*cone2 generate e2 function needed for continuum transport */

double cone2(double t)
{
	double cone2_v, 
	  remain, 
	  tln;

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

	if( t < 80. )
	{
		tln = exp(-t);
	}
	else
	{
		cone2_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->cone2()\n", debug_fp );
#		endif
		return( cone2_v );
	}

	/* fit of second exponential integral;
	 * T is optical depth, and TLN is EXP(-t)
	 * */
	if( t < 0.3 )
	{
		remain = (1.998069357 + t*(66.4037741 + t*107.2041376))/(1. + 
		  t*(37.4009646 + t*105.0388805));

	}
	else if( t < 20. )
	{
		remain = (1.823707708 + t*2.395042899)/(1. + t*(2.488885899 - 
		  t*0.00430538));

	}
	else
	{
		remain = 1.;
	}

	cone2_v = remain*tln/(2. + t);

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

/* a continuum escape probability */
static double conrec(double x)
{
	double conrec_v;

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

	conrec_v = exp(-chnukt_ContTkt*(x-1.))/x;

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

/*esccon continuum escape probability */
double esccon(double tau, 
  double hnukt)
{
	double dinc, 
	  escpcn_v, 
	  sumesc, 
	  sumrec;

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

	/* computes continuum escape probabilities */
	if( tau < 0.01 )
	{
		escpcn_v = 1.;
		
#		ifdef DEBUG_FUN
		fputs( " <->esccon()\n", debug_fp );
#		endif
		return( escpcn_v );
	}

	else if( hnukt > 1. && tau > 100. )
	{
		escpcn_v = 1e-20;
		
#		ifdef DEBUG_FUN
		fputs( " <->esccon()\n", debug_fp );
#		endif
		return( escpcn_v );
	}

	chnukt_ContTkt = hnukt;
	chnukt_ctau = tau;

	dinc = 10./hnukt;
	sumrec = qg32(1.,1.+dinc,conrec);
	sumesc = qg32(1.,1.+dinc,escConE2);

	if( sumrec > 0. )
	{
		escpcn_v = sumesc/sumrec;
	}
	else
	{
		escpcn_v = 0.;
	}

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

/*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */
static void RTesc_lya_1side(double taume, 
  double beta, 
  float *esc, 
  float *dest,
  /* position of line in frequency array on c scale */
  long ipLine )
{
	double esc0, 
	  fac, 
	  fac1, 
	  fac2, 
	  tau, 
	  taucon, 
	  taulog;

	/* DEST0 is the smallest destruction probability to return
	 * in high metallicity models, in rt.h
	const double DEST0=1e-8 ;*/

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

	/* fits to numerical results of Hummer and Kunasz Ap.J. 80 */
	tau = taume*SQRTPI;

	/* this is the real escape probability */
	esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau)));

	esc0 = MAX2(0.,esc0);
	esc0 = MIN2(1.,esc0);

	if( tau > 0. )
	{
		taulog = log10(MIN2(1e8,tau));
	}
	else
	{
		/* the line mased */
		hydro.xLaMase = MIN2(hydro.xLaMase,(float)tau);
		taulog = 0.;
		*dest = 0.;
		*esc = (float)esc0;
	}

	if( beta > 0. )
	{
		taucon = MIN2(2.,beta*tau);

		if( taucon > 1e-3 )
		{
			fac1 = -1.25 + 0.475*taulog;
			fac2 = -0.485 + 0.1615*taulog;
			fac = -fac1*taucon + fac2*POW2(taucon);
			fac = pow(10.,fac);
			fac = MIN2(1.,fac);
		}
		else
		{
			fac = 1.;
		}

		*esc = (float)(esc0*fac);
		/* MIN puts cat at 50 */
		*dest = (float)(beta/(0.30972 - MIN2(.28972,0.03541667*taulog)));
	}

	else
	{
		*dest = 0.;
		*esc = (float)esc0;
	}

	*dest = MIN2(*dest,1.f-*esc);
	*dest = MAX2(0.f,*dest);

	/* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering.
	 * in this case scattering is much more likely than absorption on this event */
	*dest = (float)( (1. - opac.albedo[ipLine]) * *dest + opac.albedo[ipLine]*DEST0);
	/* this is for debugging H Lya */
	{
		/*@-redef@*/
		enum {BUG=FALSE};
		/*@+redef@*/
		if( BUG )
		{
			fprintf(ioQQQ,"scatdest tau %.2e beta%.2e 1-al%.2e al%.2e dest%.2e \n",
			taume,
			beta, 
			(1. - opac.albedo[ipLine]), 
			opac.albedo[ipLine] ,
			*dest 
			);
		}
	}

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

/*RT_DestProb returns line destruction probability due to continuum opacity */
double RT_DestProb(
	  /* abundance of species */
	  double abund, 
	  /* its line absorption cross section */
	  double crsec, 
	  /* pointer to energy within continuum array, to get background opacity,
	   * this is on the f not c scale */
	  long int ipanu, 
	  /* line width */
	  double widl, 
	  /* escape probability */
	  double escp, 
	  /* type of redistribution function */
	  int nCore)
{
	/* this will be the value we shall return */
	double eovrlp_v;

	double conopc, 
	  beta;

	/* DEST0 is the smallest destruction probability to return
	 * in high metallicity models 
	 * this was set to 1e-8 until 99nov18,
	 * in cooling flow model the actual Lya ots dest prob was 1e-16,
	 * and this lower limit of 1e-8 caused energy balance problems,
	 * since dest prob was 8 orders of magnitude too great.  
	 * >>chng 99 nov 18, to 1e-20, but beware that comments indicate that
	 * this will cause problems with high metallicity clouds(?) */
	/* >>chng 00 jun 04, to 0 since very feeble ionization clouds, with almost zero opacity,
	 * this was a LARGE number */
	/*const double DEST0=1e-20 ;
	const double DEST0=0. ;*/
	
#	ifdef DEBUG_FUN
	fputs( "<+>RT_DestProb()\n", debug_fp );
#	endif

	/* computes "escape probability" due to continuum destruction of
	 *
	 * if esc prob gt 1 then line is masing - return small number for dest prob */
	/* >>>chng 99 apr 10, return min dest when scattering greater than abs */
	/* no idea of opacity whatsoever, on very first soln for this model */
	if( escp >= 1.0 || !conv.nTotalIoniz)
	{
		eovrlp_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->RT_DestProb()\n", debug_fp );
#		endif
		return( eovrlp_v );
	}

	/* find continuum opacity */
	conopc = opac.opacity_abs[ipanu-1];

	ASSERT( crsec > 0. );

	/* may be no population, cannot use above test since return 0 not DEST0 */
	if( abund <= 0. || conopc <= 0. )
	{
		/* do not set this to DEST0 since energy not then conserved */
		eovrlp_v = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->RT_DestProb()\n", debug_fp );
#		endif
		return( eovrlp_v );
	}

	/* fac of 1.7 convert to Hummer convention for line opacity */
	beta = conopc/(abund*SQRTPI*crsec/widl + conopc);
	/* >>chng 04 may 10, rm * 1-pesc)
	beta = MIN2(beta,(1.-escp)); */

	if( nCore == ipDEST_INCOM )
	{
		/*  fits to 
		 *  >>>refer	la	esc	Hummer and Kunasz 1980 Ap.J. 236,609.
		 *  the max value of 1e-3 is so that we do not go too far
		 *  beyond what Hummer and Kunasz did, discussed in
		 * >>refer	rt	esc proc	Ferland, G.J., 1999, ApJ, 512, 247 */
		/*TODO 2 this min is because there are no calculations that show what to do
		 * for beta beyound this value */
		eovrlp_v = MIN2(1e-3,8.5*beta);/**/
	}
	else if( nCore == ipDEST_K2 )
	{
		/*  Doppler core only; a=0., Hummer 68 
		eovrlp_v = RT_DestHummer(beta);*/
		eovrlp_v = MIN2(1e-3,8.5*beta);/**/
	}
	else if( nCore == ipDEST_SIMPL )
	{
		/*  this for debugging only 
		eovrlp_v = 8.5*beta;*/
		/* >>chng 04 may 13, use same min function */
		eovrlp_v = MIN2(1e-3,8.5*beta);/**/
	}
	else
	{
		fprintf( ioQQQ, " chCore of %i not understood by RT_DestProb.\n", 
		  nCore );
		puts( "[Stop in RT_DestProb]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* renorm to unity */
	eovrlp_v /= 1. + eovrlp_v;

	/* multiply by 1-escape prob, since no destruction when optically thin */
	eovrlp_v *= 1. - escp;

	/*check results in bounds */
	ASSERT( eovrlp_v >= 0.  );
	ASSERT( eovrlp_v <= 1.  );

	{
		/* debugging code for Lya problems */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			/*lint -e777 float test equality */
			if( rfield.anu[ipanu-1]>0.73 && rfield.anu[ipanu-1]<0.76 &&
				abund==EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc ) 
			/*lint +e777 float test equality */
			{
				fprintf(ioQQQ,"%li RT_DestProb\t%g\n",
					nzone, eovrlp_v  );
			}
		}
	}

	/* >>chng 04 may 10, rm min */
	/* this hack removed since no fundamental reason for it to be here,
	 * this should be added to scattering escape, if included at all */
#	if 0
	/* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering.
	 * in this case scattering is much more likely than absorption on this event
	eovrlp_v = (1. - opac.albedo[ipanu-1]) * eovrlp_v + 
		opac.albedo[ipanu-1]*DEST0; */
	/* >>chng 01 aug 11, add factor of 3 for increase in mean free path, and min on 0 */
	/*eovrlp_v = MAX2(DEST0,1. - 3.*opac.albedo[ipanu-1]) * eovrlp_v + 
		opac.albedo[ipanu-1]*DEST0;*/
	eovrlp_v = POW2(1. - opac.albedo[ipanu-1]) * eovrlp_v + 
		opac.albedo[ipanu-1]*0.;
#	endif

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

	return( eovrlp_v );
}

#if 0
/*compue linewidth for the special case of Lya - outer optical depths always defined */
/*RT_LyaWidth wrapper to call RT_LineWidth for the special case of Lya */
double RT_LyaWidth(
			 double tauin, 
			 double tauout, 
			 double a, 
			 /* answer will be in whatever units vth is in, since only generates
			  * optical depth dependent scale factor */
			 double vth)
{
	static int lgTonSav;
	double aa, 
	  atau, 
	  b, 
	  r, 
	  t, 
	  tau, 
	  therm, 
	  widlin_v;

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

	/* optical depth in outer direction is always defined for la */
	lgTonSav = opac.lgTauOutOn;
	opac.lgTauOutOn = TRUE;

	/* lgTauOutOn=.TRUE. if optical depth in outer direction is defined
	 * smaller of inner and outer optical depths is chosen for esc prob
	 * thermal broadening alone follows, lvg at 99
	 * following uses line width from Bonilha et al. Ap.J. (1979) 233 649
	 * return value is half width*(1-ESC PROB)
	 * TAU is used for width of source function
	 * T is used for escape probability
	 * this assumes incompl redis, a.tau^1/3 width
	 * */
	if( opac.lgTauOutOn )
	{
		tau = tauout/2.;
		t = MIN2(tauin,tauout-tauin);
		t = MAX2(tauin/100.,t);
	}
	else
	{
		tau = tauin;
		t = tauin;
	}

	/* max optical depth is thermalization length */
	therm = 5.3e16/dense.eden;
	if( tau > therm )
	{
		pressure.lgPradDen = TRUE;
		tau = therm;
	}

	/* >>chng 04 may 13, use thermal for subsonic cases */
	/*TODO	1	dynamics; this test assumes that neg vel are subsonic, so that sobolev length
	 * would overestimate the optical depth, since ion is at most present over computed
	 * slab, and possibly more.  */
	if( wind.windv <= 0 )
	{
		/* converted from go-to 15 june 95 */
		atau = log(MAX2(0.0001,tau));
		if( tau <= 20. )
		{
			aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
			b = 6.5*tau - atau;
		}
		else
		{
			aa = 1. + 2.*atau/pow(1. + 0.3*a*tau,0.6667) + pow(6.5*
			  a*tau,0.333);
			b = 1.6 + 1.5/(1. + 0.20*a*tau);
		}

		/* this is half width of line */
		widlin_v = vth*0.8862*aa/b*(1. - esc_PRD_1side(t,a));

	}

	else
	{
		/* WIND */
		r = a*tau/PI;
		if( r <= 1. )
		{
			widlin_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
			
			ASSERT( widlin_v >= 0. );
#			ifdef DEBUG_FUN
			fputs( " <->RT_LyaWidth()\n", debug_fp );
#			endif
			return( widlin_v );
		}

		widlin_v = 2.*fabs(wind.windv0);
		if( r*vth > fabs(widlin_v ) )
		{ 
			ASSERT( widlin_v >= 0. );
#			ifdef DEBUG_FUN
			fputs( " <->RT_LyaWidth()\n", debug_fp );
#			endif
			return( widlin_v );
		}

		widlin_v = vth*r*log(widlin_v/(r*vth));
	}

	/* reset the optical depths flag */
	opac.lgTauOutOn = lgTonSav;

	ASSERT( widlin_v >= 0. );

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

/*RT_LineWidth compute line width (cm/sec), using optical depth array information */
double RT_LineWidth(EmLine * t)
{
	double RT_LineWidth_v, 
	  aa, 
	  atau, 
	  b, 
	  r, 
	  tau, 
	  tauin, 
	  tauout, 
	  therm, 
	  vth;

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

	/* lgTauOutOn=.TRUE. if optical depth in outer direction is defined
	 * smaller of inner and outer optical depths is chosen for esc prob
	 * thermal broadening alone follows, lvg at 99
	 * following uses line width from 
	 * >>refer	esc	prob	Bonilha et al. Ap.J. (1979) 233 649
	 * return value is half width*(1-ESC PROB)
	 * TAU is used for width of source function
	 * T is used for escape probability
	 * this assumes incompl redis, damp.tau^1/3 width */
	tauin = t->TauIn;
	vth = DoppVel.doppler[ t->nelem -1 ];

	if( opac.lgTauOutOn )
	{
		/* find optical depth to nearest boundary */
		/*tau = tauout/2.;*/

		/* >>chng 04 may 22, tauout had been just TauTot, make tau to shielded face */
		tauout = t->TauTot-tauin;

		/* >>chng 99 apr 22 use smaller optical depth */
		tau = MIN2( tauin , tauout );
	}
	else
	{
		tau = tauin;
	}

	/* max optical depth is thermalization length */
	therm = 5.3e16/dense.eden;

	if( tau > therm )
	{
		pressure.lgPradDen = TRUE;
		tau = therm;
	}

	/* >>chng 01 jan 08, can be very slightly negative , 
	 * >>chng 04 may 22, which will cause 0 to be returned below
	tau = SDIV( tau );*/

	/* >>chng 01 jun 23, use wind vel instead of rt since rt deleted */
	/* >>chng 04 may 13, use thermal for subsonic cases */
	/*TODO	1	dynamics; this test assumes that neg vel are subsonic, so that sobolev length
	 * would overestimate the optical depth, since ion is at most present over computed
	 * slab, and possibly more.  */
	if( wind.windv <= 0. )
	{
		/* static geometry */
#		if 0
		/* >>chng 04 may 12 following tried to avoid max and float operation, but
		 * ignored Pdest and Pelec_esc */
		/* always use 1 - pesc */
		if( 0 && tau< 0.0001 )
		{
			/* >>chng 99 apr 22, this branch added for very optically thin lines */
			aa = 14.;
			b = 9.21;
			/* in very optically thin limt assume 1-escp is tau */
			/* >>chng 01 aug 24, subtract off taumin, the smallest it can be */
			RT_LineWidth_v = vth*0.8862*aa/b*MAX2(0.,tau-opac.taumin) ;
		}
#		endif
		/* >>chng 04 may 22, esc prob has noise if smaller than FLT_EPSILON, or is masing */
		if( (tau-opac.taumin)/100. < FLT_EPSILON )
		{
			RT_LineWidth_v = 0.;
		}
		else if ( tau <= 20. )
		{
			atau = log(MAX2(0.0001,tau));
			aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
			b = 6.5*tau - atau;
			/*RT_LineWidth_v = vth*0.8862*aa/b*(1. - t->Pesc);*/
			/* >>chng 01 feb 23, include dest prob since this could dominate the level pops */
			RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( (t->Pesc+ t->Pelec_esc + t->Pdest) , 1.) );
			/* small number roundff can domintate this process */
			if( t->Pesc < 10. * FLT_EPSILON )
				RT_LineWidth_v = 0.;
		}
		else
		{
			ASSERT( t->damp*tau >= 0.);
			atau = log(MAX2(0.0001,tau));
			aa = 1. + 2.*atau/pow(1. + 0.3*t->damp*tau,0.6667) + pow(6.5*
			  t->damp*tau,0.333);
			b = 1.6 + 1.5/(1. + 0.20*t->damp*tau);
			/*RT_LineWidth_v = vth*0.8862*aa/b*(1. - t->Pesc);*/
			/* >>chng 01 feb 23, include dest prob since this could dominate the level pops */
			RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( (t->Pesc+ t->Pelec_esc + t->Pdest) , 1.) );
		}

		/* we want full width, not half width */
		RT_LineWidth_v *= 2.;

	}
	else
	{
		/* wind */
		r = t->damp*tau/PI;
		if( r <= 1. )
		{
			RT_LineWidth_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
		}
		else
		{
			RT_LineWidth_v = 2.*fabs(wind.windv0);
			if( r*vth <= RT_LineWidth_v )
			{
				RT_LineWidth_v = vth*r*log(RT_LineWidth_v/(r*vth));
			}
		}
	}
	
	ASSERT( RT_LineWidth_v >= 0. );

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

/*RT_DestHummer evaluate Hummer's betaF(beta) function */
double RT_DestHummer(double beta) /* beta is ratio of continuum to mean line opacity,
									 * returns dest prob = beta F(beta) */
{
	double fhummr_v, 
	  x;

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

	/* evaluates Hummer's F(beta) function for case where damping
	 * constant is zero, are returns beta*F(beta)
	 * fit to Table 1, page 80, of Hummer MNRAS 138, 73-108.
	 * beta is ratio of continuum to line opacity; FUMMER is
	 * product of his F() times beta; the total destruction prob
	 * this beta is Hummer's normalization of the Voigt function */

	ASSERT( beta >= 0.) ;/* non-positive is unphysical */
	if( beta <= 0. )
	{
		fhummr_v = 0.;
	}
	else
	{
		x = log10(beta);
		if( x < -5.5 )
		{
			fhummr_v = 3.8363 - 0.56329*x;
		}
		else if( x < -3.5 )
		{
			fhummr_v = 2.79153 - 0.75325*x;
		}
		else if( x < -2. )
		{
			fhummr_v = 1.8446 - 1.0238*x;
		}
		else
		{
			fhummr_v = 0.72500 - 1.5836*x;
		}
		fhummr_v *= beta;
	}

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

