/home66/gary/public_html/cloudy/c13_branch/source/atom_level2.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*atom_level2 do level population and cooling for two level atom,
 * side effects:
 * set elements of transition struc
 * cooling via  CoolAdd( chLab, (long)t->WLAng , t->cool());
 * cooling derivative */
#include "cddefines.h"
#include "phycon.h"
#include "transition.h"
#include "dense.h"
#include "rfield.h"
#include "thermal.h"
#include "cooling.h"
#include "atoms.h"

void atom_level2(const TransitionProxy &t)
{
        char chLab[5];
        long int ion, 
          ip, 
          nelem;
        double AbunxIon, 
          a21, 
          boltz, 
          col12, 
          col21, 
          coolng, 
          g1, 
          g2, 
          omega, 
          pfs1,
          pfs2,
          r, 
          rate12, 
          ri21;

        DEBUG_ENTRY( "atom_level2()" );

        /* result is density (cm-3) of excited state times A21
         * result normalized to N1+N2=ABUND
         * routine increments dCooldT, call CoolAdd
         * CDSQTE is EDEN / SQRTE * 8.629E-6
         */

        ion = (*t.Hi()).IonStg();
        nelem = (*t.Hi()).nelem();

        /* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3) */
        AbunxIon = dense.xIonDense[nelem-1][ion-1];

        /* continuum pointer */
        ip = t.ipCont();

        /* approximate Boltzmann factor to see if results zero */
        boltz = rfield.ContBoltz[ip-1];

        /* collision strength for this transition, omega is zero for hydrogenic
         * species which are done in special hydro routines */
        omega = t.Coll().col_str();

        /* ROUGH check whether upper level has significant population,*/
        r = (boltz*dense.cdsqte + t.Emis().pump())/(dense.cdsqte + t.Emis().Aul());

        /* following first test needed for 32 bit cpu on search phase
         * >>chng 96 jul 02, was 1e-30 crashed on dec, change to 1e-25
         * if( AbunxIon.lt.1e-25 .or. boltz.gt.30. ) then
         * >>chng 96 jul 11, to below, since can be strong pumping when
         * Boltzmann factor essentially zero */
        /* omega in following is zero for hydrogenic species, since done
         * in hydro routines, so this should cause us to quit on this test */
        /* >>chng 99 nov 29, from 1e-25 to 1e-30 to keep same result for
         * very low density models, where AbunxIon is very small but still significant*/
        /*if( omega*AbunxIon < 1e-25 || r < 1e-25 )*/
        if( omega*AbunxIon < 1e-30 || r < 1e-25 )
        {
                /* put in pop since possible just too cool */
                (*t.Lo()).Pop() = AbunxIon;
                t.Emis().PopOpc() = AbunxIon;
                (*t.Hi()).Pop() = 0.;
                t.Emis().xIntensity() = 0.;
                t.Coll().cool() = 0.;
                t.Emis().ots() = 0.;
                t.Emis().phots() = 0.;
                t.Emis().ColOvTot() = 0.;
                t.Coll().heat() = 0.;
                /* level populations */
                atoms.PopLevels[0] = AbunxIon;
                atoms.PopLevels[1] = 0.;
                atoms.DepLTELevels[0] = 1.;
                atoms.DepLTELevels[1] = 0.;
                return;
        }

        /* net rate down A21*(escape + destruction) */
        a21 = t.Emis().Aul()*(t.Emis().Pesc()+ t.Emis().Pdest() + t.Emis().Pelec_esc());

        /* put null terminated line label into chLab using line array*/
        chIonLbl(chLab,t);

        /* statistical weights of upper and lower levels */
        g1 = (*t.Lo()).g();
        g2 = (*t.Hi()).g();

        /* now get real Boltzmann factor */
        boltz = t.EnergyK()/phycon.te;

        ASSERT( boltz > 0. );
        boltz = sexp(boltz);

        ASSERT( g1 > 0. && g2 > 0. );

        /* this lacks the upper statistical weight */
        col21 = dense.cdsqte*omega;
        /* upward coll rate s-1 */
        col12 = col21/g1*boltz;
        /* downward coll rate s-1 */
        col21 /= g2;

        /* rate 1 to 2 is both collisions and pumping */
        /* the total excitation rate from lower to upper, collisional and pumping */
        rate12 = col12 + t.Emis().pump();

        /* induced emissions down */
        ri21 = t.Emis().pump()*g1/g2;

        /* this is the ratio of lower to upper level populations */
        r = (a21 + col21 + ri21)/rate12;

        /* upper level pop */
        pfs2 = AbunxIon/(r + 1.);
        atoms.PopLevels[1] = pfs2;
        (*t.Hi()).Pop() = pfs2;

        /* pop of ground */
        pfs1 = pfs2*r;
        atoms.PopLevels[0] = pfs1;

        /* compute ratio Aul/(Aul+Cul) */
        /* level population with correction for stimulated emission */
        (*t.Lo()).Pop() = atoms.PopLevels[0];

        t.Emis().PopOpc() = (atoms.PopLevels[0] - atoms.PopLevels[1]*g1/g2 );

        /* departure coef of excited state rel to ground */
        atoms.DepLTELevels[0] = 1.;
        if( boltz > 1e-20 && atoms.PopLevels[1] > 1e-20 )
        {
                /* this line could have zero boltz factor but radiatively excited
                 * dec alpha does not obey () in fast mode?? */
                atoms.DepLTELevels[1] = (atoms.PopLevels[1]/atoms.PopLevels[0])/
                  (boltz*g2/g1);
        }
        else
        {
                atoms.DepLTELevels[1] = 0.;
        }

        /* number of escaping line photons, used elsewhere for outward beam */
        t.Emis().phots() = t.Emis().Aul()*(t.Emis().Pesc() + t.Emis().Pelec_esc())*pfs2;

        /* intensity of line */
        t.Emis().xIntensity() = t.Emis().phots()*t.EnergyErg();
        //double plower = AbunxIon - pfs2;

        /* ratio of collisional to total (collisional + pumped) excitation */
        t.Emis().ColOvTot() = col12/rate12;

        /* two cases - collisionally excited (usual case) or 
         * radiatively excited - in which case line can be a heat source
         * following are correct heat exchange, they will mix to get correct deriv 
         * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to
         * keep stable solution by effectively dividing up heating and cooling,
         * so that negative cooling does not occur */

        //double Enr12 = plower*col12*t.EnergyErg;
        //double Enr21 = pfs2*col21*t.EnergyErg;

        /* energy exchange due to this level
         * net cooling due to excit minus part of de-excit -
         * note that ColOvTot cancels out in the sum heat - cool */
        //coolng = Enr12 - Enr21*t.Emis().ColOvTot();

        /* this form of coolng is guaranteed to be non-negative */
        coolng = t.EnergyErg()*AbunxIon*col12*(a21 + ri21)/(a21 + col21 + ri21 + rate12);
        ASSERT( coolng >= 0. );

        t.Coll().cool() = coolng;

        /* net heating is remainder */
        t.Coll().heat() = t.EnergyErg()*AbunxIon*col21*(t.Emis().pump())/(a21 + col21 + ri21 + rate12);

        /* expression pre jul 3 95, changed for case where line heating dominates
         * coolng = (plower*col12 - pfs2*col21)*t.t(ipLnEnrErg)
         * t.t(ipLnCool) = cooling */

        /* add to cooling - heating part was taken out above,
         * and is not added in here - it will be added to thermal.heating[0][22]
         * in CoolSum */
        CoolAdd( chLab, t.WLAng() , t.Coll().cool());

        /* derivative of cooling function */
        thermal.dCooldT += coolng * (t.EnergyK() * thermal.tsq1 - thermal.halfte );
        return;
}

---