root/branches/newmole/source/cont_pump.cpp

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*DrvContPump local continuum pumping rate radiative transfer for all lines */
/*con_pump_op  routine used to get continuum pumping of lines 
 * used in DrvContPump in call to qg32 */
#include "cddefines.h"
#include "rfield.h"
#include "doppvel.h"
#include "radius.h"
#include "continuum.h"

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

/*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 DrvContPump(transition * t)
{
        double a0, 
          ContPump_v, 
          tau, 
          yinc1, 
          yinc2;

        DEBUG_ENTRY( "DrvContPump()" );

        /* 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 fluorescence */
                ContPump_v = 0.;
        }
        else
        {
                /* tau used will be optical depth in center of next zone */
                tau = t->Emis->TauIn + t->Emis->PopOpc * t->Emis->opacity / DoppVel.doppler[t->Hi->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->Emis->iRedisFun > 0 )
                        {
                                damp = t->Emis->damp;
                        }
                        else
                        {
                                damp = 0.;
                        }
                        PumpTau = (realnum)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) )
         * */
        return( ContPump_v );
#       undef   FITTED
}

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

        DEBUG_ENTRY( "con_pump_op()" );

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

        return( opfun_v );
}