/* 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 */
/*ContPump local continuum pumping rate radiative transfer for all lines */
/*con_pump_op  routine used to get continuum pumping of lines 
 * used in ContPump in call to qg32 */
/*vfun approximate form of Voit function */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "doppvel.h"
#include "radius.h"
#include "continuum.h"

/* damping constant used for pumping */
static float damp;
/* variable used for inward optical depth for pumping */
static float PumpTau;

/*vfun approximate form of Voit function */
static double vfun(double x);

/*con_pump_op  routine used to get continuum pumping of lines 
 * used in ContPump in call to qg32 */
static double con_pump_op(double x);

/*TODO	2	if used, add damp as arg since calling routine probably evaluated it */
double ContPump(EmLine * t)
{
	double a0, 
	  ContPump_v, 
	  tau, 
	  yinc1, 
	  yinc2;

#	ifdef DEBUG_FUN
	fputs( "<+>ContPump()\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)))

	if( !rfield.lgInducProcess )
	{
		/* option to turn off continuum pumping with no fluorescense */
		ContPump_v = 0.;
	}
	else
	{
		/* tau used will be optical depth in center of next zone
		 * >>chng 96 jul 06, had been ipLnTauIn, did not work when sphere static set */
		tau = t->TauIn + t->PopOpc * t->opacity / DoppVel.doppler[t->nelem-1]*radius.dRNeff;
		/* compute pumping probability */
		if( tau <= 10. )
		{
			/* for tau<10 a does not matter, and one fit for all */
			ContPump_v = FITTED(tau);
		}
		else if( tau > 1e6 )
		{
			/* this far in winds line opacity well below electron scattering
			 * so ignore for this problem */
			ContPump_v = 0.;
		}
		else
		{
			/* following two are passed on to later subs */
			if( t->iRedisFun > 0 )
			{
				damp = t->damp;
			}
			else
			{
				damp = 0.;
			}
			PumpTau = (float)tau;
			a0 = 0.886227*(1. + damp);
#			define	BREAK_	3.
			yinc1 = qg32(0.,BREAK_,con_pump_op);
			yinc2 = qg32(BREAK_,100.,con_pump_op);
			ContPump_v = (yinc1 + yinc2)/a0;
		}
	}

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

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

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

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

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


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

/*vfun approximate form of Voit function */
static double vfun(double x)
{
	double vfun_v;

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

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

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