/home66/gary/public_html/cloudy/c08_branch/source/ion_nitro.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 */
/*IonNitro ionization balance for nitrogen */
#include "cddefines.h"
#include "dense.h"
#include "thermal.h"
#include "atoms.h"
#include "opacity.h"
#include "trace.h"
#include "conv.h"
#include "gammas.h"
#include "ionbal.h"

void IonNitro(void)
{
        const int NDIM = ipNITROGEN+1;

        static const double dicoef[2][NDIM] = {
                {2.98e-3,7.41e-3,1.13e-2,2.62e-3,7.5e-2,4.61e-2,0.},
                {0.,.0764,.164,.243,.35,.309,0.}
        };
        static const double dite[2][NDIM] = {
                {2.2e5,2.01e5,1.72e5,1.02e5,4.75e6,5.44e6,0.},
                {0.,7.37e4,2.25e5,1.25e5,8.35e5,1.14e6,0.}
        };
        static const double ditcrt[NDIM] = {1.8e4,1.8e4,2.4e4,1.5e4,6.8e5,1.0e6,1e20};
        static const double aa[NDIM] = {0.0,0.0320,-0.8806,0.4134,0.,0.,0.};
        static const double bb[NDIM] = {0.6310,-0.6624,11.2406,-4.6319,0.,0.,0.};
        static const double cc[NDIM] = {0.1990,4.3191,30.7066,25.9172,0.,0.,0.};
        static const double dd[NDIM] = {-0.0197,0.0003,-1.1721,-2.2290,0.,0.,0.};
        static const double ff[NDIM] = {0.4398,0.5946,0.6127,0.2360,0.1,0.,0.};

        bool lgDebug;
        double save_rec;

        DEBUG_ENTRY( "IonNitro()" );

        /* nitrogen, atomic number 7 */
        if( !dense.lgElmtOn[ipNITROGEN] )
        {
                return;
        }

        ion_zero(ipNITROGEN);

        ion_photo(ipNITROGEN,false);

        /* find collisional ionization rates */
        ion_collis(ipNITROGEN);

        /* get recombination coefficients */
        /*lint -e64 type mismatch */
        ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipNITROGEN);
        /*lint +e64 type mismatch */

        /* >>chng 04 jul 09, synch up ion and co solvers.
         * the co solver will have a different N/N+ balance that
         * is strongly affected by the chemistry if we are in neutral gas
         * in this case use co solver's N/N+ balance */
        save_rec = ionbal.RateRecomTot[ipNITROGEN][0];
        /*>>chng 04 apr 24, do not test for ionhigh being 1,
         * no reason for test on upper stage of ionization,
         * if codrive called but not evaluated then hevmol is all zero */
        /* >>chng 04 sep 10, rm check on search phase, no reason for it */
#       if 0
        if( dense.IonLow[ipNITROGEN]==0 && 
                co.hevmol[ipNP] > SMALLFLOAT && 
                ionbal.RateIonizTot[ipNITROGEN][0]*co.hevmol[ipATN]> 0. )
        {
                ionbal.RateRecomTot[ipNITROGEN][0] = 
                        ionbal.RateIonizTot[ipNITROGEN][0]*
                        co.hevmol[ipATN]/co.hevmol[ipNP];
        }
#       endif

        /* photoionization from 2D of NI, is atomic nitrogen present? */
        if( dense.xIonDense[ipNITROGEN][0] > 0. )
        {
                atoms.d5200r = (realnum)GammaK(opac.in1[0],opac.in1[1],opac.in1[2],1.);
                /*aeff = atoms.d5200r + 1.28e-5;*/
                /* >>chng 01 sep 02, use values from coolnitr */
                /*cs5200 = 1.32e-4*phycon.te/(phycon.te10*phycon.te01);*/
                /* uses 1 for total pop of atom - so is normalized to unity */
                /*atoms.p2nit = (realnum)(atom_pop2(cs5200,4.,10.,aeff,2.769e4,1.)/aeff);*/
                /* if very first call then set to zero */
                if( !conv.nTotalIoniz ) 
                        atoms.p2nit = 0.;
                /* valence shell photoionization, [0][6] => atomic nitrogen */
                ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0]*
                  (1. - atoms.p2nit) + atoms.p2nit*atoms.d5200r;
                ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1]*
                  (1. - atoms.p2nit) + thermal.HeatNet*atoms.p2nit;
        }
        else
        {
                atoms.p2nit = 0.;
                atoms.d5200r = 0.;
        }

        /* solve for ionization balance */
#       if 0
        broken();/* rm following true */
        if(nzone == 1 )
                lgDebug = true;
        else
                lgDebug = false;
#       endif
        lgDebug = false;
        ion_solver(ipNITROGEN,lgDebug);

        /* reset the var we just hosed */
        if( save_rec > 0.)
                ionbal.RateRecomTot[ipNITROGEN][0] = save_rec;

        if( trace.lgTrace && trace.lgHeavyBug )
        {
                fprintf( ioQQQ, "     IonNitro retun; frac=" );
                for( int i=0; i < 8; i++ )
                {
                        fprintf( ioQQQ, "%10.3e", dense.xIonDense[ipNITROGEN][i]/
                          dense.gas_phase[ipNITROGEN] );
                }
                fprintf( ioQQQ, "\n" );
        }
        return;
}

---