/*esccom fundamental escape probability radiative transfer routine, for complete redistribution */
/*escinc fundamental escape probability radiative transfer routine for incomplete redistribution */
/*escla escape prob for hydrogen atom Lya, using Hummer and Kunasz results,
 * called by hydropesc */
/*escsub escape probability radiative transfer for subordinate lines */
/*escgrd escape probability radiative transfer for incomplete redistribution */
/*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
/*escmase escape probability for negative (masing) optical depths,
 * only called by routines here */
double escmase(double tau);

#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "escdst.h"
#include "zonecnt.h"
#include "showme.h"
#include "twoway.h"
#include "conrec.h"
#define	SCALE	2.
#include "chnukt.h"
#include "touton.h"
#include "taulines.h"
#include "rtescprob.h"
#include "qg32.h"

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

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


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

	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( " <->escinc()\n", debug_fp );
#	endif
	return( escinc_v );
}

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

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

	/* this is one of the three fundamental escape probability routines
	 * the three are esccom, escinc, and escla
	 * it computes esc prob for complete redistribution
	 * computes escape prob for complete redistribution in one direction
	 * */
	if( tau < 0. )
	{
		/* line mased */
		esccom_v = escmase(tau);
		
#		ifdef DEBUG_FUN
		fputs( " <->esccom()\n", debug_fp );
#		endif
		return( esccom_v );
	}

	if( strcmp(escdst.chEscFun,"SIMP") == 0 )
	{
		/* this set with "escape simple" command, used for debugging */
		esccom_v = 1./(1. + tau);
		
#		ifdef DEBUG_FUN
		fputs( " <->esccom()\n", debug_fp );
#		endif
		return( esccom_v );
	}

	/* chEscFun is the escape prob function, usually k2
	 * */
	if( strcmp(escdst.chEscFun,"INCO") == 0 )
	{
		if( tau < 10. )
		{
			esccom_v = 1./(1. + 1.6*tau);
			
#			ifdef DEBUG_FUN
			fputs( " <->esccom()\n", debug_fp );
#			endif
			return( esccom_v );
		}

		/* 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);

		esccom_v = 1./(1. + b*tau) + 0.25*sqrt(a/tau);

	}
	else if( strcmp(escdst.chEscFun,"  K2") == 0 )
	{
		/* this is the usual case, Hummer's k2 function
		 * error! what happens if negative tau? */
		esccom_v = esca0k2(tau);
	}
	else
	{
		fprintf( ioQQQ, " Unknown escape probability function=%4.4s\n", 
		  escdst.chEscFun );
		puts( "[Stop in esccom]" );
		exit(1);
	}

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

/*escla escape prob for hydrogen atom Lya, using Hummer and Kunasz results,
 * called by hydropesc */
#include "doppvel.h"
#include "opac.h"
#include "redis.h"
#include "humla.h"

double escla(float *esin, 
  float *dest, 
  double abund, 
  long int ipZ)/* 0 for H, */
{
	double beta, 
	  conopc, 
	  damp, 
	  escla_v;
	 float dstin, 
	  dstout;

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

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

	/* check charge */
	assert( ipZ >= 0 );
	assert( ipZ < LIMELM );

	if( HydroLines[ipZ][IP2P][IP1S].TauTot - 
	  HydroLines[ipZ][IP2P][IP1S].TauIn < 0. )
	{
		/* this is the case if we overrun the optical depth scale
		 * just leave things as they are */
		escla_v = HydroLines[ipZ][IP2P][IP1S].Pesc;
		twoway.fracin = HydroLines[ipZ][IP2P][IP1S].FracInwd;
		*esin = twoway.fracin;
		*dest = HydroLines[ipZ][IP2P][IP1S].Pdest;
		
#		ifdef DEBUG_FUN
		fputs( " <->escla()\n", debug_fp );
#		endif
		return( escla_v );
	}

	if( strcmp(redis.chRedis,"INCO") == 0 )
	{
		/* incomplete redistribution */
		conopc = opac.opac[HydroLines[ipZ][IP2P][IP1S].ipCont-1];
		if( abund > 0. )
		{
			/* the continuous opacity is positive, we have a valid soln */
			beta = conopc/(abund/SQRTPI*HydroLines[ipZ][IP2P][IP1S].opacity/
			  DoppVel.doppler[ipZ] + conopc);
		}
		else
		{
			/* abundance is zero, set miniumum dest prob */
			beta = 1e-10;
		}

		/* WAYIN will be escape prob in inward direction */
		humla(HydroLines[ipZ][IP2P][IP1S].TauIn,
		  beta,&twoway.wayin,&dstin , 
		  HydroLines[ipZ][IP2P][IP1S].ipCont-1);

		if( (twoway.wayin > 1.) || (twoway.wayin < 0.) || (dstin > 1.) || (dstin < 0.) )
		{
			fprintf( ioQQQ, 
				" HUMLA returns inwd prob out of bounds;  esc=%10.2e dst="
				"%10.2e beta=%10.2e tau=%10.2e\n", 
			  twoway.wayin, dstin, beta, HydroLines[ipZ][IP2P][IP1S].TauIn );
			puts( "[Stop in escla]" );
			exit(1);
		}

		humla(MAX2(HydroLines[ipZ][IP2P][IP1S].TauTot/
		  100.,HydroLines[ipZ][IP2P][IP1S].TauTot-
		  HydroLines[ipZ][IP2P][IP1S].TauIn),
		  beta,&twoway.wayout,&dstout, 
		  HydroLines[ipZ][IP2P][IP1S].ipCont-1);

		if( (twoway.wayout > 1.) || (twoway.wayout < 0.) || (dstout > 1.) || (dstout < 0.) )
		{
			fprintf( ioQQQ, 
				" HUMLA returns outwd prob out of bounds;  esc=%10.2e dst=%10.2e beta=%10.2e tau=%10.2e\n", 
			  twoway.wayout, dstout, beta, 
			  MAX2(HydroLines[ipZ][IP2P][IP1S].TauTot/
			  100.,HydroLines[ipZ][IP2P][IP1S].TauTot-
			  HydroLines[ipZ][IP2P][IP1S].TauIn) );
			puts( "[Stop in escla]" );
			exit(1);
		}

		escla_v = (twoway.wayin + twoway.wayout)/2.;
		*esin = twoway.wayin;
		*dest = (dstin + dstout)/2.f;
		twoway.fracin = twoway.wayin/(twoway.wayin + twoway.wayout);
	}

	else if( (strcmp(redis.chRedis,"COMP") == 0 ) ||
		     (strcmp(redis.chRedis,"WIND") == 0 ) )
	{
		/* wind or complete redis */
		/*>>>chng 98 dec 18, set damp to zero before call, had not
		 * been initialized */
		damp = 0.;
		escla_v = escsub(HydroLines[ipZ][IP2P][IP1S].TauIn,
		  HydroLines[ipZ][IP2P][IP1S].TauTot,
		  damp);
		*esin = twoway.wayin;
		twoway.fracin = 0.5;
		*dest = 0.;
		/* not used, but keep lint quiet */
		beta = 1e-10;
	}

	else
	{
		/* not used, but keep lint quiet */
		beta = 1e-10;
		fprintf( ioQQQ,"escla called with insane redistribution function ==%s==\n",
			redis.chRedis );
		puts( "[Stop in escla]" );
		exit(1);
	}

	assert( escla_v >= 0. );

	/* this is for debugging H Lya */
	{
		enum {BUG=FALSE};
		if( BUG )
		{
			if( ipZ == 0 )
			{
				fprintf(ioQQQ," %.2e %.2e %.2e %.2e %.2e \n",
				HydroLines[ipZ][IP2P][IP1S].TauIn,
				beta, 
				opac.albedo[HydroLines[ipZ][IP2P][IP1S].ipCont-1] ,
				escla_v , 
				*dest 
				);
			}
		}
	}

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

/*escsub escape probability radiative transfer for subordinate lines */

double escsub(double tau, 
  double tout, 
  double damp)
{
	double escsub_v, 
	  t;

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

	if( touton.lgTauOutOn )
	{
		/* outward optical depths are defined, but don't let them go negat */
		t = tout - 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( t < 0. )
			t = tau/SCALE;
		twoway.wayin = (float)esccom(tau,damp);
		twoway.wayout = (float)esccom(t,damp);
		twoway.fracin = twoway.wayin/(twoway.wayin + twoway.wayout);
		escsub_v = 0.5*(twoway.wayin + twoway.wayout);
	}
	else
	{
		/* first iteration, outward optical depths not defined
		 * dont estimate fraction in */
		twoway.fracin = 0.;
		twoway.wayout = 1.;
		escsub_v = esccom(tau,damp);
		twoway.wayin = (float)escsub_v;
	}


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

/*escgrd escape probability radiative transfer for incomplete redistribution */

double escgrd(double tau, 
  double tout, 
  double damp )
{
	double escgrd_v, 
	  tt;

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

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

	if( touton.lgTauOutOn )
	{
		/*  outward optical depth if defined */
		tt = tout - 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;
		}

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

	assert( escgrd_v > 0. );

#	ifdef DEBUG_FUN
	fputs( " <->escgrd()\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 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 */
#include "dumpline.h"

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 the major heat agent */
		if( TauLines[i].TauIn < -1. )
			DumpLine(&TauLines[i]);
	}

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

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

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", ZoneCnt.nzone );
		FindNeg();
		ShowMe();
		puts( "[Stop in escmase]" );
		exit(1);
	}

	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 );
}

/*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 = (float)hnukt;
	chnukt.ctau = (float)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 );
}