root/branches/newmole/source/atom_level2.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 */
/*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 "lines_service.h"
#include "dense.h"
#include "rfield.h"
#include "thermal.h"
#include "cooling.h"
#include "atoms.h"

void atom_level2(transition * t)
{
        char chLab[5];
        long int ion, 
          ip, 
          nelem;
        double AbunxIon, 
          a21, 
          boltz, 
          col12, 
          col21, 
          coolng, 
          g1, 
          g2, 
          omega, 
          pfs1,
          pfs2,
          plower, 
          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;
        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 = Enr21*(1. - t->Emis->ColOvTot);
        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;
}