/home66/gary/public_html/cloudy/c08_branch/source/rt_line_one.cpp

Go to the documentation of this file.

/* 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 */
/*RT_line_static do line radiative transfer,
 * evaluates escape and destruction probability, 
 * called by HydroPEsc, RT_line_all  */
/* ND wind distinction is not here - that is done where optical depths are incremented
 * with line opacity - rt_line_one_tauinc and in line width for rad pressure, in 
 * rt escprob */
#include "cddefines.h"
#include "rfield.h"
#include "doppvel.h"
#include "dense.h"
#include "opacity.h"
#include "lines_service.h"
#include "conv.h"
#include "radius.h"
#include "rt.h"
#include "physconst.h"

/*RT_line_static do line radiative transfer,
 * called by RT_line_all */
STATIC void RT_line_static(
        /* the em line we will work on  */
        transition * t , 
        /* when this is true do both esc and dest probs,
         * when falst, only dest probs */
        bool lgDoEsc ,
        /* this is option to not include line self shielding across this zone.
         * this can cause pump to depend on zone thickness, and leads to unstable
         * feedback in some models with the large H2 molecule, due to Solomon
         * process depending on zone thickness and level populations. */
        bool lgShield_this_zone );

/*RT_line_one do rt for emission line structure - calls RT_line_static or RT_line_wind */
void RT_line_one(
        /* the em line we will work on  */
        transition *t , 
        /* when this is true do both esc and dest probs,
         * when falst, only dest probs */
        bool lgDoEsc ,
        /* if true then we want to update the fine opacity scale */
        bool lgDoFine_opac_update ,
        /* this is option to not include line self shielding across this zone.
         * this can cause pump to depend on zone thickness, and leads to unstable
         * feedback in some models with the large H2 molecule, due to Solomon
         * process depending on zone thickness and level populations. */
        bool lgShield_this_zone )
{
                static long int nTau[100],
                        n;

        DEBUG_ENTRY( "RT_line_one()" );

        /*>>chng 06 apr 05,  effects of simply returning upon call with no pop */
        /* skip line transfer if requested with 'no line transfer' command, but never skip Lya */
        if( (t->Lo->Pop<=SMALLFLOAT && conv.nTotalIoniz) || 
                (!rfield.lgDoLineTrans && (t->Emis->iRedisFun != ipLY_A) ) )
        {
                return;
        }

        {
                /* option to keep track of population values during calls,
                 * print out data to make histogram */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        if( nzone==0 )
                        {
                                for(n=0; n<100; ++n )
                                        nTau[n] = 0;
                        }
                        if( t->Lo->Pop<=SMALLFLOAT )
                                n = 0;
                        else
                                n = (long)log10( t->Lo->Pop )+37;
                        n = MIN2( n , 99 );
                        n = MAX2( n , 0 );
                        /*if( n==37 && nzone > 5 )
                        {
                                fprintf(ioQQQ,"HIT\n");
                        }*/
                        ++nTau[n];
                        if( nzone > 183 )
                        {
                                for(n=0; n<100; ++n )
                                        fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] );
                                cdEXIT(EXIT_SUCCESS);
                        }
                }
        }

        /* always call single routine */
        RT_line_static(
                /* the line we will work on  */
                t , 
                /* when this is true do both esc and dest probs,
                        * when falst, only dest probs */
                lgDoEsc ,
                lgShield_this_zone);

        /* this is line center frequency, including bulk motion of gas */
        long int ipLineCenter = t->Emis->ipFine + rfield.ipFineConVelShift;

        /* define fine opacity fine grid fine mesh */
        /* rfield.lgOpacityFine flag set flast with no fine opacities command */
        if( lgDoFine_opac_update && ipLineCenter >= 0 && t->Emis->PopOpc > 0. &&
                ipLineCenter<rfield.nfine && rfield.lgOpacityFine )
        {
                long int i;

                /* number of fine opacity cells corresponding to one doppler width for current
                 * temperature, velocity field, and nuclear mass,
                 * rfield.fine_opac_velocity_width is width per cell, cm/s */
                realnum cells_wide_1x = DoppVel.doppler[t->Hi->nelem-1]/rfield.fine_opac_velocity_width;

                /* line center opacity - want realnum since will add to fine opacity array,
                 * which is realnum */
                realnum opac_line =  (realnum)t->Emis->PopOpc * t->Emis->opacity / DoppVel.doppler[t->Hi->nelem-1];

                /* this is effective optical depth to this point. Do not work with line if this product
                 * is substantially less than unity */
                double dTauEffec = opac_line*radius.depth_x_fillfac;

                /* this is line center for current velocity, which is measured
                * relative to illuminated face */

                /* use dTauEffec to choose whether line is important */

                if( dTauEffec > SMALLFLOAT )
                {
                        /* width of optically thick line, do 4x with exponential core,
                         * must be at least one cell, but profile is symmetric */
                        long int nCells_core = (long)(cells_wide_1x*4.f+1.5f);
                        realnum opacity;

                        /* core is symmetric - make sure both upper and lower bounds
                         * are within continuum energy bin */
                        if( ipLineCenter - nCells_core < 1 )
                                nCells_core = ipLineCenter - 1;
                        if( ipLineCenter + nCells_core > rfield.nfine )
                                nCells_core = ipLineCenter -rfield.nfine - 1;

                        /* want core to be at least one cell wide */
                        nCells_core = MAX2( 1 , nCells_core );

                        /* line center itself, must not double count here */
                        rfield.fine_opac_zone[ipLineCenter] += opac_line;
                        for( i=1; i<nCells_core; ++i )
                        {
                                /* square of distance from line center in units of doppler width */
                                realnum xsq = POW2(i/cells_wide_1x);

                                /* the line opacity at that point, includes doppler core and damping wings */
                                opacity = opac_line*( sexp(xsq) + t->Emis->damp / MAX2(1.f,(1.772f*xsq)) );

                                rfield.fine_opac_zone[ipLineCenter+i] += opacity;
                                rfield.fine_opac_zone[ipLineCenter-i] += opacity;
                        }

                        /* include damping wings if optical depth is large */
                        /* >>chng 04 mar 29, rewrite to use symmetry in two halves of line */
                        if( dTauEffec*t->Emis->damp/9. > 0.1 )
                        {
                                /* do damping wings if very optically thick
                                 * this is number of cells to extend the damping wings - cells_wide is one dop width 
                                 * 56.419 is 1 / (0.01 * sqrt pi) ) */
                                /* tests with th85orion, stopping at half -h2 point, showed that 0.01 was
                                 * needed for outer edge, given the definition of dTauEffec */
                                realnum x = (realnum)sqrt( dTauEffec * t->Emis->damp *56.419 ) * cells_wide_1x;
                                long int nCells_damp;
                                /* >>chng 04 mar 24, add test on size of x, which can exceed
                                * limits of long in extreme optical depths */
                                if( x<LONG_MAX )
                                {
                                        nCells_damp = (long)x;

                                        if( ipLineCenter-nCells_damp < 1 )
                                                nCells_damp = ipLineCenter-1;

                                        if( ipLineCenter+nCells_damp > rfield.nfine )
                                                nCells_damp = rfield.nfine - ipLineCenter-1;
                                }
                                else
                                {
                                        /* x was too big, just set to extreme range, which is
                                         * number of cells to nearest boundary */
                                        nCells_damp = MIN2( rfield.nfine-ipLineCenter , ipLineCenter )-1;
                                }

                                ASSERT( nCells_core>0 );
                                for( i=nCells_core; i<nCells_damp; ++i )
                                {
                                        /* distance from line center in doppler units */
                                        realnum xsq = POW2(i/cells_wide_1x);
                                        /* term in denominator is x^2 */
                                        opacity = opac_line*t->Emis->damp/(1.772f*xsq );
                                        rfield.fine_opac_zone[ipLineCenter+i] += opacity;
                                        rfield.fine_opac_zone[ipLineCenter-i] += opacity;
                                }
                        }
                }
        }
        return;
}

/*RT_line_static do line radiative transfer for static geometry */
STATIC void RT_line_static(
        /* the em line we will work on  */
        transition * t , 
        /* when this is true do both esc and dest probs,
         * when false, only dest probs */
        bool lgDoEsc ,
        /* this is option to not include line self shielding across this zone.
         * this can cause pump to depend on zone thickness, and leads to unstable
         * feedback in some models with the large H2 molecule, due to Solomon
         * process depending on zone thickness and level populations. */
        bool lgShield_this_zone )
{
        bool lgGoodTau;
        long int nelem;
        int nRedis = -1;

        DEBUG_ENTRY( "RT_line_static()" );

        /* this is the routine that computes much of the radiative transfer
         * physics for lines for static case.  */

        /*ASSERT( t->IonStg < LIMELM+4 );*/
        /*ASSERT( t->nelem-1 < LIMELM+4 );*/

        /* >>chng 01 aug 18, do not evaluate if no population */
        /* molecules are evaluated here, these must be valid */
        /* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
         * that rates are evaluated at start of a new iteration */
        if( dense.xIonDense[ t->Hi->nelem-1][t->Hi->IonStg-1] == 0.  && !lgDoEsc )
        {
                /* zero population, return after setting everything with side effects */
                t->Emis->Pesc = 1.f;

                /* inward escaping fraction */
                t->Emis->FracInwd = 0.5;

                /* pumping rate */
                t->Emis->pump = 0.;

                /* destruction probability */
                t->Emis->Pdest = 0.;
                t->Emis->Pelec_esc = 0.;

                return;
        }
        /* check that we don't have a runaway maser */
        else if( t->Emis->TauIn < -30. )
        {
                fprintf( ioQQQ, "PROBLEM RT_line_static called with large negative optical depth, zone %.2f, setting lgAbort true.\n",
                        fnzone );
                DumpLine(t);
                /* >>>chng 00 apr 23, return busted true instead of exit here, so that code will
                 * not hang under MPI */
                lgAbort = true;
                return;
        }

        /* this checks if we have overrun the optical depth scale,
         * in which case the inward optical depth is greater than the
         * previous iteration's total optical depth.
         * We do not reevaluate escape probabilities if the optical depth
         * scale has been overrun due to huge bogus change in solution
         * that would result */
        if( iteration == 1 || lgTauGood( t ) )
        {
                lgGoodTau = true;
        }
        else
        {
                /* do not reevaluate escape probabilities if we have overrun 
                 * the optical depth scale */
                lgGoodTau = false;
        }

        /* atomic number of this element on the C scale*/
        nelem = t->Hi->nelem-1;

        /* damping constant for the line, save it */
        t->Emis->damp = t->Emis->dampXvel / DoppVel.doppler[t->Hi->nelem-1];
        /*if( t->damp <0. )
        {
                fprintf(ioQQQ,"DEBUG hit it \n");
        }*/
        ASSERT( t->Emis->damp > 0. || 
                        (opac.lgCaseB && t->Emis->damp >= 0.) );

        /* static solution - which type of line will determine
         * which redistribution function */
        /* iRedisFun == 1 - alpha resonance line, partial redistribution,
         * ipPRD == 1 */
        if( t->Emis->iRedisFun == ipPRD )
        {
                /* incomplete redistribution with wings */
                if( lgDoEsc )
                {
                        if( lgGoodTau )
                        {
                                t->Emis->Pesc = (realnum)esc_PRD( t->Emis->TauIn , t->Emis->TauTot , t->Emis->damp );

                                /* inward escaping fraction */
                                t->Emis->FracInwd = rt.fracin;
                        }
                }
                nRedis = ipDEST_INCOM;
        }

        /* complete redistribution without wings - t->ipLnRedis is ipCRD == -1 */
        else if( t->Emis->iRedisFun == ipCRD )
        {
                if( lgDoEsc )
                {
                        if( lgGoodTau )
                        {
                                /* >>chng 01 mar -6, escsub will call any of several esc prob routines,
                                 * depending of how core is set.  We always want core-only for this option,
                                 * so call  esca0k2(tau) directly */
                                t->Emis->Pesc = (realnum)esc_CRDcore( t->Emis->TauIn , t->Emis->TauTot );

                                /* inward escaping fraction */
                                t->Emis->FracInwd = rt.fracin;
                        }
                }
                nRedis = ipDEST_K2;
        }

        /* chng 01 feb 27, add CRD with wings, = 2 */
        else if( t->Emis->iRedisFun == ipCRDW )
        {
                /* complete redistribution with damping wings */
                if( lgDoEsc )
                {
                        if( lgGoodTau )
                        {
                                t->Emis->Pesc = (realnum)esc_CRDwing( t->Emis->TauIn , t->Emis->TauTot , t->Emis->damp );

                                /* inward escaping fraction */
                                t->Emis->FracInwd = rt.fracin;
                        }
                }
                nRedis = ipDEST_K2;
        }

        /* chng 03 may 14, special Lya case*/
        else if( t->Emis->iRedisFun == ipLY_A )
        {
                /* incomplete redistribution with wings, for special case of Lya
                 * uses fits to Hummer & Kunasz numerical results 
                 * this routine is different because escape and dest probs
                 * are evaluated together, so no test of lgDoEsc */
                if( lgGoodTau )
                {
                        double dest , esin;

                        /* this will always evaluate escape prob, no matter what lgDoEsc is.
                        * The destruction prob comes back as dest */
                        t->Emis->Pesc = (realnum)RTesc_lya(&esin , &dest , t->Emis->PopOpc , nelem);

                        /* this is current destruction rate */
                        t->Emis->Pdest = (realnum)dest;

                        /*  this is fraction of line which is inward escaping */
                        t->Emis->FracInwd = rt.fracin;
                }
        }
        else
        {
                fprintf( ioQQQ, " RT_line_static called with redistribution function zero\n" );
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }

        /* pumping by incident and diffuse continua */
        /* >>chng 03 may 15, this block had been repeated all four times above, move down
         * here to use common code, only rt.wayin is defined above */
        /* >>chng 03 may 17, do not evaluate when drad is zero, before first
         * zone thickness is established */
        /* >>chng 06 nov 29, reordered if clauses, no point in doing bigger loop if "no induced" is on. */
        /* option to kill induced processes */
        if( !rfield.lgInducProcess )
        {
                t->Emis->pump = 0.;
        }
        else if( (lgDoEsc  || (t->Emis->iRedisFun == ipLY_A) ) )
        {
                double dTau;
                double pumpsave = t->Emis->pump;
                double shield_continuum;

                /* >>chng 04 apr 18, break out line continuum shielding into this function */
                shield_continuum = RT_continuum_shield_fcn( t );

                /* continuum upward pumping rate, A gu/gl abs prob occnum
                 * the "no induced" command causes continuum pumping to be set to 0 */
                /*t->pump = t->Aul * t->gHi / t->gLo * rt.wayin **/
                /* >>chng 05 feb 16, add outward diffuse emission */
                /*t->pump = t->Aul * t->gHi / t->gLo * shield_continuum *
                        rfield.OccNumbIncidCont[t->ipCont-1];*/
                t->Emis->pump = t->Emis->Aul * t->Hi->g / t->Lo->g * shield_continuum *(
                        rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1] );

                /* >>chng 99 dec 02, add correction for line optical depth across zone */
                /* correction for attenuation within line across zone */
                /* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow, 
                 * this caused an overestimate of pumping rates */
                /* >>chng 00 dec 19, nova.in oscillated since correction factor was
                 * very large across zones that were already very optically thick, dr changed,
                 * causing feedback between ionization and thickness
                 * the intention of this correction was to only apply when on the 
                 * verge of optically thick,
                 * change lower bound on product t->dTau * radius.drad_x_fillfac to 1e-3 to avoid
                 * constantly getting unity, also check on relative change optical depth/optical depth*/
                /* this can drive oscillations in very neutral gas, as blr92.in, blr89.in
                 * or nova.in, when higher lyman lines have produce 
                 * t->dTau * radius.drad_x_fillfac about unity, where small changes in
                 * pump can change ionization, dTau, and hence pump */
                /*if( t->dTau * radius.drad_x_fillfac > 1e-8 )*/
                /* >>chng 02 jun 14, nova.in again - upper limit had been 10 - but
                 * there were major jumps in heating rate at this point - increase 10 -> 1e10
                 * not clear why there should be an upper limit at all, other than to save time? */
                /* update dTau since affects pump rate */
                /* >>chng 03 jun 18, include option to not shield across this zone */
                /* >>chng 03 jul 19, add continuum opacity */
                dTau =  (t->Emis->PopOpc * t->Emis->opacity / DoppVel.doppler[nelem] + opac.opacity_abs[t->ipCont-1])*
                         radius.drad_x_fillfac_mean;
                /* >>chng 04 mar 04, do not want to include maser action across this zone */
                dTau = MIN2(1., dTau);
                /*if( 0 && lgShield_this_zone && dTau * radius.drad_x_fillfac > 1e-3 && dTau * radius.drad_x_fillfac < 1e10 )*/
                /* >>chng 04 mar 04, use mean of past few zone thicknesses rather than
                 * current zone thickness since this can induce wild oscillations, see
                 * nova.in test case */
                if( lgShield_this_zone && dTau > 1e-3 && dTau < 1e10 )
                {
                        /* during very first search dTau is -1e38*/
                        /* >>chng 04 mar 04, use mean, see above comment */
                        /*t->pump *= log(1. + dTau * radius.drad_x_fillfac ) / (dTau * radius.drad_x_fillfac);*/
                        t->Emis->pump *= log(1. + dTau ) / dTau;
                }

                /* >>chng 00 Oct 03, add pumping by local diffuse continuum,
                 * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in  parsedont.c
                 * with no DIFFuse LINE PUMPing command */
                if( dTau*rfield.DiffPumpOn > SMALLFLOAT )
                {
                        /* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so
                         * that smallest and largest number are multiplied first to prevent overflow, PvH */
                        double scale;
                        if( iteration > 1 )
                        {
                                if( lgGoodTau )
                                {
                                        /* >>chng 03 jul 19, add check for tau = 1 in continuum */
                                        /* scale for tau = 1 in the continuum */
                                        double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]);
                                        /* inward looking depth, dTau is opacity in line, cm^-1,
                                        * so multiplying by this scale is same as dividing by opacity in 
                                        * source function */
                                        double scale_in =  MIN3( 1./ dTau , radius.depth    , scale_con);
                                        double scale_out = MIN3( 1./ dTau , radius.Depth2Go , scale_con );
                                        t->Emis->pump += t->Emis->Aul * t->Hi->g / t->Lo->g * 0.5 *((scale_in+scale_out)*
                                                rfield.OccNumbDiffCont[t->ipCont-1]);
                                }
                        }
                        else 
                        {
                                /* >>chng 03 jul 19, add check for tau = 1 in continuum */
                                /* scale for tau = 1 in the continuum */
                                double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]);
                                /* first iteration, only inward looking depths really known */
                                scale = 1./ dTau;
                                scale = MIN3( scale , radius.depth , scale_con );
                                t->Emis->pump += t->Emis->Aul * t->Hi->g / t->Lo->g * (scale*
                                        rfield.OccNumbDiffCont[t->ipCont-1]);
                        }
                }
                /* this stop oscillations from developing in very deep blr models 
                 * like blr89 and blr98, but must not be applied too early since
                 * line optical depths change greatly in first few zones */
                if( nzone>50 )
                        t->Emis->pump = (pumpsave + t->Emis->pump) / 2.;
        }


        /* escape following scattering off an electron */
        /* this is turned off with the no scattering escape command */
        if( !rt.lgElecScatEscape )
        {
                t->Emis->Pelec_esc = 0.;
        }
        else
        {
                /* in the transition structure PopOpc is the pop relative to the ion,
                 * but has been converted to a physical population before this routine
                 * was called.  No need to convert here */
                double opac_line = t->Emis->PopOpc * t->Emis->opacity/DoppVel.doppler[nelem];

                if( opac_line > SMALLFLOAT )
                {
                        /* the old escape probability */
                        double eesc_save = t->Emis->Pelec_esc;
                        /* the opacity in electron scattering */
                        double es = dense.eden*6.65e-25;
                        fixit(); // should we use total opacity here instead of opac_line?
                        /* this is equation 5 of 
                         *>>refer       line    desp    Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752*/
                        double opacity_ratio = es/(es+opac_line) * (1.-t->Emis->Pesc);
                        /* keep total probability not more than 1. */
                        t->Emis->Pelec_esc = (realnum)opacity_ratio * MAX2(0.f,1.f-t->Emis->Pesc-t->Emis->Pdest);
                        /*KLUDGE >>chng 03 feb 02 
                         * take mean of old and new if electron scattering escape */
                        if( !conv.lgSearch )
                                t->Emis->Pelec_esc = t->Emis->Pelec_esc*0.75f + (realnum)eesc_save*0.25f;  
                }
                else
                        t->Emis->Pelec_esc = 0.;
        }

        /* >>>chng 99 dec 18, did not have the test for lgGoodTau, and so clobbered
         * dest prob when overrun */
        /* only do this if not Lya special case, since dest prob already done */
        if( lgGoodTau && t->Emis->iRedisFun!=ipLY_A && t->Emis->opacity > 0. )
        {
                double DestSave = t->Emis->Pdest;
                t->Emis->Pdest = (realnum)RT_DestProb(
                        /* the abundance of the species */
                        t->Emis->PopOpc,
                        /* line center opacity in funny units (needs a vel) */
                        t->Emis->opacity,
                        /* array index for line in continuum array,
                         * on f not c scale */
                        t->ipCont,
                        /* line width in velocity units */
                        DoppVel.doppler[nelem],
                        /* escape probability */
                        t->Emis->Pesc,
                        /* redistribution function */
                        nRedis);

                /* >>chng 03 may 10, this had been 3/4 new + 1/4 old, but had to be combined
                 ^ with half and half when this routine was called, to converge very high Z high U
                 * quasar models.  bring all these factors together here */
                /*t->Pdest = 0.75f*t->Pdest + 0.25f * DestSave;*/
#               define OLDFAC   (0.625)
                if( nzone>1 )
                        t->Emis->Pdest = (realnum)((1.-OLDFAC)*t->Emis->Pdest + OLDFAC * DestSave);
        }

        /* zero escape probability and pumping term if energy is below plasma frequency */
        if( t->EnergyErg / EN1RYD < rfield.plsfrq )
        {
                t->Emis->Pesc = 0.;
                t->Emis->Pdest = 0.;
                t->Emis->pump = 0.;
        }

        return;
}

#if 0
/*RT_line_wind do line radiative transfer for wind geometry */
STATIC void RT_line_wind(transition * t , bool lgDoEsc);
/*RT_line_wind do line radiative transfer for wind geometry */
STATIC void RT_line_wind(transition * t , bool lgDoEsc)
{
        realnum velocity;

        DEBUG_ENTRY( "RT_line_wind()" );

        /* 
         * this is the routine that computes much of the radiative transfer
         * physics for lines in wind case.  Lines diveded into two types,
         * high quality (level1) and lower quality (atom_level2) g-bar lines
         */

        /* >>chng 01 aug 18, do not evaluation if no population */
        /* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
         * that rates are evaluated at start of a new iteration */
        if( dense.xIonDense[ t->nelem-1][t->IonStg-1] == 0.  && !lgDoEsc )
        {
                /* zero population, return after setting everything with side effects */
                t->Pesc = 1.f;

                /* inward escaping fraction */
                t->FracInwd = 0.5;

                /* pumping rate */
                t->pump = 0.;

                /* destruction probability */
                t->Pdest = 0.;
                t->Pelec_esc = 0.;
                return;
        }

        /*confirm pops and cs ok */
        ASSERT( t->ColOvTot <= 1. && t->ColOvTot >= 0. );

        /* t->nelem is atomic number of species */
        ASSERT( t->nelem >= 1 && t->ipCont > 0 );

        /* line width, thermal + turbulence */
        velocity = DoppVel.doppler[t->nelem-1];

        t->damp = 0.;

        /* use only one-sided escape probabilities
         * get escape probability, this gets fractions too */
        if( lgDoEsc )
        {
                t->Pesc = esc_CRDwing_1side(t->TauIn,t->damp);

                /* inward fraction */
                t->FracInwd = 0.5;

                /* continuum pumping */
                if( rfield.lgInducProcess )
                {
                        /* >>chng 05 feb 16, include diffuse outward emission */
                        /*t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*
                          rfield.OccNumbIncidCont[t->ipCont-1];*/
                        t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*(
                          rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1]);
                }
                else
                {
                        t->pump = 0.;
                }
        }

        /* get destruction probability */
        t->Pdest = 
                RT_DestProb(
                t->PopOpc,
            /* its line absorption cross section */
                t->opacity,
                /* index for line in continuum array, on f not c scale */
                t->ipCont,
            velocity,
                t->Pesc, 
                ipDEST_K2);
        return;
}
#endif

---