/* 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_continuum_shield_fcn computing continuum shielding due to single line */
/*conpmp local continuum pumping rate radiative transfer for all lines */
/*vfun voit function helper routine */
/*conpmp local continuum pumping rate radiative transfer for all lines */
#include "cddefines.h"
#include "physconst.h"
#include "rt.h"
/*conpmp local continuum pumping rate radiative transfer for all lines */
static double conpmp(EmLine * t);
static float PumpDamp , PumpTau;

/*rt_continuum_shield_fcn computing continuum shielding due to single line */
float RT_continuum_shield_fcn( EmLine *t )
{
	float value;

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

	value = -1.f;

	if( rt.nLineContShield == LINE_CONT_SHIELD_PESC )
	{
		/* set continuum shielding pesc - the default, 1 */
		if( t->iRedisFun == ipPRD )
		{
			value =  (float)esc_PRD_1side(t->TauCon,t->damp);
		}
		else if( t->iRedisFun == ipCRD )
		{
			value =  (float)esca0k2(t->TauCon);
		}
		else if( t->iRedisFun == ipCRDW )
		{
			value =  (float)esc_CRDwing_1side(t->TauCon,t->damp);
		}
		else if( t->iRedisFun == ipLY_A )
		{
			value = (float)esc_PRD_1side(t->TauCon,t->damp);
		}
		else
			TotalInsanity();
	}
	else if( rt.nLineContShield == LINE_CONT_SHIELD_FEDERMAN )
	{
		/* set continuum shielding federman */
		float core, wings;
		/* this option used if command "set continuum shielding federman" appears */

		/* these expressions implement the appendix of
		 * >>refer	line	shielding	Federman, S.R., Glassgold, A.E., & 
		 * >>refercon	Kwan, J. 1979, ApJ, 227, 466 */
		/* doppler core - equation A8 */
		if( t->TauCon < 2. )
		{
			core = (float)sexp( t->TauCon * 0.66666 );
		}
		else if( t->TauCon < 10. )
		{
			core = (float)(0.638 * pow(t->TauCon,-1.25f ));
		}
		else if( t->TauCon < 100. )
		{
			core = (float)(0.505 * pow(t->TauCon,-1.15f ));
		}
		else
		{
			core = (float)(0.344 * pow(t->TauCon,-1.0667f ));
		}

		/* do we add damping wings? */
		wings = 0.;
		/*if( t->TauCon>1e3 )*/
		/* >>chng 04 jun 14, evaluate if TauCon > 1, to make pump fcn smoother */
		if( t->TauCon > 1.f && t->damp>0. )
		{
			/* equation A6 */
			float t1 = (float)(3.02*pow(t->damp*1e3,-0.064 ) );
			float u1 = (float)(sqrt(t->TauCon*t->damp )/SDIV(t1));
			/* >>chng 04 jun 09, Gargi Shaw found error in following */
			/*wings = (t->damp/t->TauCon/(float)sqrt( 0.78540 + POW2(u1) ) );*/
			wings = (t->damp/SDIV(t1)/(float)sqrt( 0.78540 + POW2(u1) ) );
		}
		value = core + wings;
		/* some x-ray lines have vastly large damping constants, greater than 1.
		 * in these cases the previous wings value does not work - approximation
		 * is for small damping constant - do not let pump efficiency exceed unity
		 * in this case */
		if( t->TauCon>0. )
			value = MIN2(1.f, value );
	}
	else if( rt.nLineContShield == LINE_CONT_SHIELD_FERLAND )
	{
		/* set continuum shielding ferland */
		value = (float)conpmp( t );
	}
	else if( rt.nLineContShield == 0 )
	{
		/* set continuum shielding none */
		value = 1.f;
	}
	else
	{
		TotalInsanity();
	}

	/* the returned pump shield function must be greater than zero,
	 * and less than 1 if a maser did not occur */
	ASSERT( value>=0 && (value<=1.||t->TauCon<0.) );

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

	return value;
}

/*vfun voit function helper routine */
static double vfun(double x)
{
	double vfun_v;

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

	vfun_v = sexp(x*x) + PumpDamp/SQRTPI/(1. + x*x);

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

/*opfun  routine used to get continuum pumping of lines 
 * used in conpmp in call to qg32 */
static double opfun(double x)
{
	double opfun_v, 
	  v;

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

	v = vfun(x);
	opfun_v = sexp(PumpTau*v)*v;

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

#define	BREAK_	3.
/*conpmp local continuum pumping rate radiative transfer for all lines */
static double conpmp(EmLine * t)
{
	double a0, 
	  conpmp_v, 
	  tau, 
	  yinc1, 
	  yinc2;

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


	/* fit to results for tau less than 10 */
#define FITTED(t)	((0.98925439 + 0.084594094*(t))/(1. + (t)*(0.64794212 + (t)*0.44743976)))

	/* tau used will be optical depth in center of next zone
	 * >>chng 96 july 6, had been ipLnTauIn, did not work when sphere static set */
	tau = t->TauCon;
	/* compute pumping probability */
	if( tau <= 10. )
	{
		/* for tau<10 a does not matter, and one fit for all */
		conpmp_v = FITTED(tau);
	}
	else if( tau > 1e6 )
	{
		/* this far in winds line opacity well below electron scattering
			* so ignore for this problem */
		conpmp_v = 0.;
	}
	else
	{
		/* following two are passed on to later subs */
		PumpDamp = t->damp;
		PumpTau = (float)tau;
		a0 = 0.886227*(1. + PumpDamp);
		yinc1 = qg32(0.,BREAK_,opfun);
		yinc2 = qg32(BREAK_,100.,opfun);
		conpmp_v = (yinc1 + yinc2)/a0;
	}

	/* EscProb is escape probability, will not allow conpmp to be greater than it
	 * on second iteration with thick lines, pump prob=1 and esc=0.5
	 * conpmp = MIN( conpmp , t->t(ipLnEscP) )
	 * */

#	ifdef DEBUG_FUN
	fputs( " <->conpmp()\n", debug_fp );
#	endif
	return( conpmp_v );
#undef	FITTED
}

