/home66/gary/public_html/cloudy/c08_branch/source/hydrot2low.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 */
/*HydroT2Low called to do hydrogenic level populations when temp too low for matrix,
 * forces total of level populations to add up to iso.xIonSimple[ipISO][nelem],
 * so results of this routine will always agree with that value */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "secondaries.h"
#include "ionbal.h"
#include "dense.h"
#include "trace.h"
#include "phycon.h"
#include "hydrogenic.h"

void HydroT2Low( long int ipISO, long int nelem )
{
        long int iup, 
          level, 
          low;
        double entries, 
          exits, 
          sum, 
          renorm,
          collider;

        DEBUG_ENTRY( "HydroT2Low()" );

        /* check that we were called with valid charge */
        ASSERT( nelem >= 0);
        ASSERT( nelem < LIMELM );
        /* >>chng 05 aug 17, use eden as collider for H and corrected eden for heavier
         * nuclei - in hot PDR H is collisionally ionized, but not by other H0 since
         * collision is homonuclear */
        if( nelem==ipHYDROGEN )
        {
                /* special version for H0 onto H0 */
                collider = dense.EdenHontoHCorr;
        }
        else
        {
                collider = dense.EdenHCorr;
        }

        /* do radiative and downward collisional cascades */
        /* 06 aug 28, from numLevels_max to _local. */
        for( level=iso.numLevels_local[ipISO][nelem]-1; level >= 0; level-- )
        {
                /* captures into this level from the continuum,
                 * radiative rec plus three body rec, units s-1 */
                entries = iso.RateCont2Level[ipISO][nelem][level];

                /* sum over all levels between current level and top of atom */
                /* 06 aug 28, from numLevels_max to _local. */
                for( iup=level + 1; iup < iso.numLevels_local[ipISO][nelem]; iup++ )
                {
                        /* radiative decays from upper level to here, s-1 */
                        double RadDecay;

                        RadDecay = MAX2( iso.SmallA, 
                                Transitions[ipISO][nelem][iup][level].Emis->Aul*
                                (Transitions[ipISO][nelem][iup][level].Emis->Pesc  + 
                                Transitions[ipISO][nelem][iup][level].Emis->Pelec_esc + 
                                Transitions[ipISO][nelem][iup][level].Emis->Pdest) );

                        /* rate upper level (rel to continuum) decays to here, units s-1 */
                        entries += StatesElem[ipISO][nelem][iup].Pop*(RadDecay + 
                          Transitions[ipISO][nelem][iup][level].Coll.ColUL*collider);
                }

                /* now find rate this levels decays to all lower levels and continuum 
                 * continuum ionization rate, s-1, first */
                exits = iso.gamnc[ipISO][nelem][level];
                for( low=ipH1s; low <= (level - 1); low++ )
                {
                        /* radiative decays to all lower levels - s-1 */
                        double RadDecay;

                        RadDecay = MAX2( iso.SmallA, 
                                Transitions[ipISO][nelem][level][low].Emis->Aul*
                                (Transitions[ipISO][nelem][level][low].Emis->Pesc  + 
                                Transitions[ipISO][nelem][level][low].Emis->Pelec_esc + 
                                Transitions[ipISO][nelem][level][low].Emis->Pdest) );

                        /* add rad and collisional decays */
                        exits += RadDecay + 
                                Transitions[ipISO][nelem][level][low].Coll.ColUL*collider;
                }
                if( exits > 1e-25 )
                {
                        StatesElem[ipISO][nelem][level].Pop = entries/exits;

                }
                else
                {
                        StatesElem[ipISO][nelem][level].Pop = 0.;
                }
        }

        /* =================================================================== 
         *
         * at this point all destruction processes have been established 
         *
         * ==================================================================== */

        /* do simple sums for atom to ion, then fill in for rest */
        StatesElem[ipISO][nelem][ipH1s].Pop = 
                ionbal.RateRecomTot[nelem][nelem-ipISO]/
                iso.RateLevel2Cont[ipISO][nelem][ipH1s];

        if( ipISO==ipH_LIKE )
        {

                StatesElem[ipISO][nelem][ipH2s].Pop = 
                        (iso.RadRec_caseB[ipISO][nelem]*0.33*dense.eden + 
                        StatesElem[ipISO][nelem][ipH1s].Pop*secondaries.Hx12[ipISO][nelem][ipH2s] + 
                        StatesElem[ipISO][nelem][ipH2p].Pop*
                        Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider)/
                        (Transitions[ipISO][nelem][ipH2s][ipH1s].Emis->Aul + 
                        Transitions[ipISO][nelem][ipH2s][ipH1s].Coll.ColUL*collider + 
                        Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider*6. );

                StatesElem[ipISO][nelem][ipH2p].Pop = 
                        (iso.RadRec_caseB[ipISO][nelem]*0.67*dense.eden + 
                        StatesElem[ipISO][nelem][ipH1s].Pop*secondaries.Hx12[ipISO][nelem][ipH2p] + 
                        StatesElem[ipISO][nelem][ipH2s].Pop*
                        Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider*6.)/
                        (Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Aul*
                        (Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pesc  + 
                        Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pelec_esc +
                        Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pdest) +
                        Transitions[ipISO][nelem][ipH2p][ipH1s].Coll.ColUL * collider + 
                        Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL * collider);
        }


        /* now renormalize to simple ionization ratio calculated in hydrolevel */
        sum = 0.;
        /* 06 aug 28, from numLevels_max to _local. */
        for( level=0; level < iso.numLevels_local[ipISO][nelem]; level++ )
        {
                sum += StatesElem[ipISO][nelem][level].Pop;
        }

        /* rescale so that inverse of sum is iso.xIonSimple[ipISO][nelem] */
        renorm = 1. / SDIV( iso.xIonSimple[ipISO][nelem]) / SDIV( sum);

        for( level=0; level < iso.numLevels_local[ipISO][nelem]; level++ )
        {
                StatesElem[ipISO][nelem][level].Pop *= renorm;
                ASSERT( StatesElem[ipISO][nelem][level].Pop<BIGFLOAT );
                if( iso.PopLTE[ipISO][nelem][level] > 1e-25 )
                {
                        iso.DepartCoef[ipISO][nelem][level] = 
                          StatesElem[ipISO][nelem][level].Pop/
                          (iso.PopLTE[ipISO][nelem][level]*dense.eden);
                }
                else
                {
                        iso.DepartCoef[ipISO][nelem][level] = 0.;
                }
        }

        if( trace.lgHBug && trace.lgTrace )
        {
                fprintf( ioQQQ, "       LOW TE,=%10.3e HN(1)=%10.3e rec=%10.3e Hgamnc(1s)=%10.3e\n", 
                  phycon.te, StatesElem[ipISO][nelem][ipH1s].Pop, ionbal.RateRecomTot[nelem][nelem-ipISO], 
                  iso.gamnc[ipISO][nelem][ipH1s] );
        }

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        "       HydroT2Low return, LOW TE used, HII/HI=%.3e simple=%.3e IonRate=%.3e RecCo=%.3e\n", 
                  iso.pop_ion_ov_neut[ipISO][nelem], 
                  iso.xIonSimple[ipISO][nelem], 
                  iso.RateLevel2Cont[ipISO][nelem][ipH1s], 
                  ionbal.RateRecomTot[nelem][nelem-ipISO] );
        }
        return;
}

---