/home66/gary/public_html/cloudy/c08_branch/source/atom_oi.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 */
/*atom_oi drive the solution of OI level populations, Ly-beta pumping */
/*oi_level_pops get OI level population with Ly-beta pumping */
#include "cddefines.h"
#include "taulines.h"
#include "doppvel.h"
#include "iso.h"
#include "trace.h"
#include "dense.h"
#include "rt.h"
#include "rfield.h"
#include "phycon.h"
#include "lines_service.h"
#include "thirdparty.h"
#include "atoms.h"

/*oi_level_pops get OI level population with Ly-beta pumping */
STATIC void oi_level_pops(double abundoi, 
  double *coloi);

/*atom_oi drive the solution of OI level populations, Ly-beta pumping */
void atom_oi_calc(double *coloi)
{
        long int i;
        double esin, 
          eslb, 
          esoi, 
          flb, 
          foi, 
          opaclb, 
          opacoi, 
          tin, 
          tout, 
          xlb, 
          xoi;
        static double esab;
        static double aoi = TauLines[ipTO1025].Emis->Aul;
        static double alb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul;

        DEBUG_ENTRY( "atom_oi_calc()" );

        /* A's from Pradhan; OI pump line; Ly beta, 8446 */

        /* Notation;
         * PMPH31 net rate hydrogen level 3 depopulated to 1
         * PMPO15 and PMPO51 are similar, but for oxygen
         * ESCH31 is emission in 1025 transition */

        /* called by LINES before calc really start, protect here
         * also used for cases where OI not present */
        if( dense.xIonDense[ipOXYGEN][0] <= 0. )
        {
                for( i=0; i < 6; i++ )
                {
                        atoms.popoi[i] = 0.;
                }
                /* return zero */
                *coloi = 0.;
                atoms.pmpo15 = 0.;
                atoms.pmpo51 = 0.;
                atoms.pmph31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul*
                        (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc +
                  Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc);

                atoms.esch31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul*
                        (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc +
                        Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc);

                /* all trace output turned on with "trace ly beta command' */
                if( trace.lgTr8446 && trace.lgTrace )
                {
                        fprintf( ioQQQ, 
                                "       P8446 called for first time, finds A*escape prob 3-1 =%10.2e\n", 
                          atoms.esch31 );
                }
                return;
        }

        if( iteration > 1 )
        {
                /* two sided escape prob */
                tin = TauLines[ipTO1025].Emis->TauIn;
                esin = esc_CRDwing_1side(tin,1e-4);
                tout = (TauLines[ipTO1025].Emis->TauTot)*0.9 - 
                  TauLines[ipTO1025].Emis->TauIn;

                if( trace.lgTr8446 && trace.lgTrace )
                {
                        fprintf( ioQQQ, "       P8446 tin, tout=%10.2e%10.2e\n", tin, tout );
                }

                if( tout > 0. )
                {
                        tout = TauLines[ipTO1025].Emis->TauTot - TauLines[ipTO1025].Emis->TauIn;

                        /* do not update esab if we overran optical depth scale */
                        esab = 0.5*(esin + esc_CRDwing_1side(tout,1e-4));
                }
        }
        else
        {
                /* one-sided escape probability */
                esab = esc_CRDwing_1side(TauLines[ipTO1025].Emis->TauIn,1e-4);
        }

        esoi =TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc;
        eslb =Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + 
                Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc;

        /* all trace output turned on with "trace ly beta command' */
        if( trace.lgTr8446 && trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        "       P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n", 
                  DoppVel.doppler[0], DoppVel.doppler[7], eslb, esoi, esab );
        }

        /* find relative opacities for two lines */
        opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/DoppVel.doppler[7];
        opaclb = 1.22e-8*dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/DoppVel.doppler[0];

        /* these are x sub a (OI) and x sub b (ly beta) defined in Elitz+Netz */
        xoi = opacoi/(opacoi + opaclb);
        xlb = opaclb/(opacoi + opaclb);

        /* find relative line-widths, assume same rest freq */
        foi = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipOXYGEN];
        flb = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipHYDROGEN]*
          MAX2(0.,1.- TauLines[ipTO1025].Emis->Pelec_esc - TauLines[ipTO1025].Emis->Pesc);

        if( trace.lgTr8446 && trace.lgTrace )
        {
                fprintf( ioQQQ, 
                        "       P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n", 
                  opaclb, opacoi, xlb, xoi, flb, foi );
        }

        /* pumping of OI by L-beta - this goes into OI matrix as 1-5 rate
         * lgInducProcess set false with no induced command, usually true */
        if( rfield.lgInducProcess )
        {
                atoms.pmpo15 = (realnum)(flb*alb*(StatesElem[ipH_LIKE][ipHYDROGEN][ipH3p].Pop*dense.xIonDense[ipHYDROGEN][1])*xoi*(1. - esab)/
                  dense.xIonDense[ipOXYGEN][0]);
                /* net decay rate from upper level */
                atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - 
                  esab)*foi));
        }
        else
        {
                atoms.pmpo15 = 0.;
                atoms.pmpo51 = 0.;
        }

        /* find level populations for OI */
        oi_level_pops(dense.xIonDense[ipOXYGEN][0],coloi);

        /* escape term from n=3 of H; followed by pump term to n=3 */
        //atoms.pmph31 = alb*(1. - (1. - flb)*(1. - eslb) - xlb*(1. - esab)*flb);

        /*pmph13 = xlb*(1. - esab)*foi*aoi*atoms.popoi[4];*/
        atoms.pmph31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul*
                (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc);

        /* following is actual emission rate, used to predict Ly-beta in lines.for */
        atoms.esch31 = alb*(eslb*(1. - flb) + esab*flb);

        /* all trace output turned on with "trace ly beta command' */
        if( trace.lgTr8446 && trace.lgTrace )
        {
                fprintf( ioQQQ, "       P8446 PMPH31=%10.2e EscP3-1%10.2e ESCH31=%10.2e\n", 
                  atoms.pmph31, atoms.pmph31/alb, atoms.esch31/alb );
        }

        /* continuum pumping due to J=1, 0 sub states.
         * neglect J=2 since L-Beta very optically thick */

        /* lower level populations */
        TauLines[ipT1304].Lo->Pop = atoms.popoi[0];
        TauLines[ipTO1025].Lo->Pop = atoms.popoi[0];
        TauLines[ipT1039].Lo->Pop = atoms.popoi[0];
        TauLines[ipT8446].Lo->Pop = atoms.popoi[1];
        TauLines[ipT4368].Lo->Pop = atoms.popoi[1];
        TauLines[ipTOI13].Lo->Pop = atoms.popoi[2];
        TauLines[ipTOI11].Lo->Pop = atoms.popoi[2];
        TauLines[ipTOI29].Lo->Pop = atoms.popoi[3];
        TauLines[ipTOI46].Lo->Pop = atoms.popoi[4];

        /* upper level populations */
        TauLines[ipT1304].Hi->Pop = atoms.popoi[1];
        TauLines[ipTO1025].Hi->Pop = atoms.popoi[4];
        TauLines[ipT1039].Hi->Pop = atoms.popoi[3];
        TauLines[ipT8446].Hi->Pop = atoms.popoi[2];
        TauLines[ipT4368].Hi->Pop = atoms.popoi[5];
        TauLines[ipTOI13].Hi->Pop = atoms.popoi[3];
        TauLines[ipTOI11].Hi->Pop = atoms.popoi[4];
        TauLines[ipTOI29].Hi->Pop = atoms.popoi[5];
        TauLines[ipTOI46].Hi->Pop = atoms.popoi[5];

        /* opacity factor
         * t1304(ipLnEmis.PopOpc) = (popoi(1)-popoi(2)*3.0) */
        TauLines[ipT1304].Emis->PopOpc = atoms.popoi[0];
        TauLines[ipTO1025].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[4]*0.6);
        TauLines[ipT1039].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[3]*3.0);

        /* t8446(ipLnEmis.PopOpc) = (popoi(2)-popoi(3)*.33) */
        TauLines[ipT8446].Emis->PopOpc = 
                (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33));
        TauLines[ipT4368].Emis->PopOpc = (atoms.popoi[1] - atoms.popoi[5]*.33);
        TauLines[ipTOI13].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[3]*3.0);
        TauLines[ipTOI11].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[4]*0.6);
        TauLines[ipTOI29].Emis->PopOpc = (atoms.popoi[3] - atoms.popoi[5]*.33);
        TauLines[ipTOI46].Emis->PopOpc = (atoms.popoi[4] - atoms.popoi[5]*1.67);
        return;
}

/*oi_level_pops get OI level population with Ly-beta pumping */
STATIC void oi_level_pops(double abundoi, 
  double *coloi)
{
        bool lgNegPop;

        long int i, j;

        int32 ipiv[6], ner;

        double a21, 
          a32, 
          a41, 
          a43, 
          a51, 
          a53, 
          a62, 
          a64, 
          a65, 
          c12, 
          c13, 
          c14, 
          c15, 
          c16, 
          c21, 
          c23, 
          c24, 
          c25, 
          c26, 
          c31, 
          c32, 
          c34, 
          c35, 
          c36, 
          c41, 
          c42, 
          c43, 
          c45, 
          c46, 
          c51, 
          c52, 
          c53, 
          c54, 
          c56, 
          c61, 
          c62, 
          c63, 
          c64, 
          c65, 
          cs, 
          deptoi[6], 
          e12, 
          e23, 
          e34, 
          e45, 
          e56, 
          simple;

        double amat[6][6], 
          bvec[6], 
          zz[7][7];

        static realnum g[6]={9.,3.,9.,3.,15.,9};

        DEBUG_ENTRY( "oilevl()" );

        /* following used for linpac matrix inversion */

        /* compute emission from six level OI atom*/

        /* Boltzmann factors for ContBoltz since collisions not dominant in UV tran
         * ipoiex is array lof dEnergy for each level, set in DoPoint */
        e12 = rfield.ContBoltz[atoms.ipoiex[0]-1];
        e23 = rfield.ContBoltz[atoms.ipoiex[1]-1];
        e34 = rfield.ContBoltz[atoms.ipoiex[2]-1];
        e45 = rfield.ContBoltz[atoms.ipoiex[3]-1];
        e56 = rfield.ContBoltz[atoms.ipoiex[4]-1];

        /* total rad rates here have dest by background continuum */
        a21 = TauLines[ipT1304].Emis->Aul*(TauLines[ipT1304].Emis->Pdest+ TauLines[ipT1304].Emis->Pesc + TauLines[ipT1304].Emis->Pelec_esc);
        a41 = TauLines[ipT1039].Emis->Aul*(TauLines[ipT1039].Emis->Pdest+ TauLines[ipT1039].Emis->Pesc + TauLines[ipT1039].Emis->Pelec_esc);
        a51 = TauLines[ipTO1025].Emis->Aul*(TauLines[ipTO1025].Emis->Pdest+ TauLines[ipTO1025].Emis->Pesc + TauLines[ipTO1025].Emis->Pelec_esc);
        a51 = atoms.pmpo51;
        a32 = TauLines[ipT8446].Emis->Aul*(TauLines[ipT8446].Emis->Pdest+ TauLines[ipT8446].Emis->Pesc + TauLines[ipT8446].Emis->Pelec_esc);
        a62 = TauLines[ipT4368].Emis->Aul*(TauLines[ipT4368].Emis->Pdest+ TauLines[ipT4368].Emis->Pesc + TauLines[ipT4368].Emis->Pelec_esc);
        a43 = TauLines[ipTOI13].Emis->Aul*(TauLines[ipTOI13].Emis->Pdest+ TauLines[ipTOI13].Emis->Pesc + TauLines[ipTOI13].Emis->Pelec_esc);
        a53 = TauLines[ipTOI11].Emis->Aul*(TauLines[ipTOI11].Emis->Pdest+ TauLines[ipTOI11].Emis->Pesc + TauLines[ipTOI11].Emis->Pelec_esc);
        a64 = TauLines[ipTOI29].Emis->Aul*(TauLines[ipTOI29].Emis->Pdest+ TauLines[ipTOI29].Emis->Pesc + TauLines[ipTOI29].Emis->Pelec_esc);
        a65 = TauLines[ipTOI46].Emis->Aul*(TauLines[ipTOI46].Emis->Pdest+ TauLines[ipTOI46].Emis->Pesc + TauLines[ipTOI46].Emis->Pelec_esc);

        /* even at density of 10^17 excited states not in lte due
         * to fast transitions down - just need to kill a21 to get to unity at 10^17*/

        /* the 2-1 transition is 1302, cs Wang and McConkey '92 Jphys B 25, 5461 */
        cs = 2.151e-5*phycon.te/phycon.te03;
        PutCS(cs,&TauLines[ipT1304]);

        /* the 5-1 transition is 1027, cs Wang and McConkey '92 Jphys B 25, 5461 */
        cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01;
        PutCS(cs,&TauLines[ipTO1025]);
        c21 = dense.cdsqte*TauLines[ipT1304].Coll.cs/g[1];
        c51 = dense.cdsqte*TauLines[ipTO1025].Coll.cs/g[4];

        /* all following are g-bar approx, g-bar = 0.2 */
        c31 = dense.cdsqte*1.0/g[2];
        PutCS(0.27,&TauLines[ipT1039]);
        c41 = dense.cdsqte*TauLines[ipT1039].Coll.cs/g[3];
        c61 = dense.cdsqte*1./g[5];

        c12 = c21*g[1]/g[0]*e12;
        c13 = c31*g[2]/g[0]*e12*e23;
        c14 = c41*g[3]/g[0]*e12*e23*e34;
        c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
        c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;

        c32 = dense.cdsqte*85./g[2];
        c42 = dense.cdsqte*85./g[3];
        c52 = dense.cdsqte*85./g[4];
        c62 = dense.cdsqte*85./g[5];

        c23 = c32*g[2]/g[1]*e23;
        c24 = c42*g[3]/g[1]*e23*e34;
        c25 = c52*g[4]/g[1]*e23*e34*e45;
        c26 = c62*g[5]/g[1]*e23*e34*e45*e56;

        c43 = dense.cdsqte*70./g[3];
        c53 = dense.cdsqte*312./g[4];
        c63 = dense.cdsqte*1./g[5];

        c34 = c43*g[3]/g[2]*e34;
        c35 = c53*g[4]/g[2]*e34*e45;
        c36 = c63*g[5]/g[2]*e34*e45*e56;

        c54 = dense.cdsqte*50./g[4];
        c64 = dense.cdsqte*415./g[5];

        c45 = c54*g[4]/g[3]*e45;
        c46 = c64*g[5]/g[3]*e45*e56;

        c65 = dense.cdsqte*400./g[5];
        c56 = c65*g[5]/g[4]*e56;

        /* this is check for whether matrix inversion likely to fail */
        simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
        if( simple < 1e-19 )
        {
                atoms.popoi[0] = (realnum)abundoi;
                for( i=1; i < 6; i++ )
                {
                        atoms.popoi[i] = 0.;
                }
                *coloi = 0.;
                return;
        }

        /*--------------------------------------------------------- */

        for( i=0; i < 6; i++ )
        {
                zz[i][0] = 1.0;
                zz[6][i] = 0.;
        }

        /* first equation is sum of populations */
        zz[6][0] = abundoi;

        /* level two, 3s 3So */
        zz[0][1] = -c12;
        zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
        zz[2][1] = -c32 - a32;
        zz[3][1] = -c42;
        zz[4][1] = -c52;
        zz[5][1] = -c62 - a62;

        /* level three */
        zz[0][2] = -c13;
        zz[1][2] = -c23;
        zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
        zz[3][2] = -c43 - a43;
        zz[4][2] = -c53 - a53;
        zz[5][2] = -c63;

        /* level four */
        zz[0][3] = -c14;
        zz[1][3] = -c24;
        zz[2][3] = -c34;
        zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
        zz[4][3] = -c54;
        zz[5][3] = -c64 - a64;

        /* level five */
        zz[0][4] = -c15 - atoms.pmpo15;
        zz[1][4] = -c25;
        zz[2][4] = -c35;
        zz[3][4] = -c45;
        zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
        zz[5][4] = -c65 - a65;

        /* level six */
        zz[0][5] = -c16;
        zz[1][5] = -c26;
        zz[2][5] = -c36;
        zz[3][5] = -c46;
        zz[4][5] = -c56;
        zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;

        /* this one may be more robust */
        for( j=0; j < 6; j++ )
        {
                for( i=0; i < 6; i++ )
                {
                        amat[i][j] = zz[i][j];
                }
                bvec[j] = zz[6][j];
        }

        ner = 0;

        getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner);
        getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner);

        /*DGETRF(6,6,(double*)amat,6,ipiv,&ner);*/
        /*DGETRS('N',6,1,(double*)amat,6,ipiv,bvec,6,&ner);*/
        if( ner != 0 )
        {
                fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* now put results back into z so rest of code treates only
                * one case - as if matin1 had been used */
        for( i=0; i < 6; i++ )
        {
                zz[6][i] = bvec[i];
        }

        lgNegPop = false;
        for( i=0; i < 6; i++ )
        {
                atoms.popoi[i] = (realnum)zz[6][i];
                if( atoms.popoi[i] < 0. )
                        lgNegPop = true;
        }

        /* following used to confirm that all dep coef are unity at
         * density of 1e17, t=10,000, and all A's set to zero */
        if( trace.lgTrace && trace.lgTr8446 )
        {
                deptoi[0] = 1.;
                deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]*
                  e12);
                deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]*
                  e12*e23);
                deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]*
                  e12*e23*e34);
                deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]*
                  e12*e23*e34*e45);
                deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]*
                  e12*e23*e34*e45*e56);

                fprintf( ioQQQ, " oilevl finds levl pop" );
                for(i=0; i < 6; i++)
                        fprintf( ioQQQ, "%11.3e", atoms.popoi[i] );
                fprintf( ioQQQ, "\n" );

                fprintf( ioQQQ, " oilevl finds dep coef" );
                for(i=0; i < 6; i++)
                        fprintf( ioQQQ, "%11.3e", deptoi[i] );
                fprintf( ioQQQ, "\n" );
        }

        /* this happens due to numerical instability in matrix inversion routine */
        if( lgNegPop )
        {
                atoms.nNegOI += 1;

                fprintf( ioQQQ, " OILEVL finds negative population" );
                for(i=0;i < 6; i++)
                        fprintf( ioQQQ, "%10.2e", atoms.popoi[i] );
                fprintf( ioQQQ, "\n" );

                fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );

                atoms.popoi[5] = 0.;
                atoms.popoi[4] = (realnum)(abundoi*(c15 + atoms.pmpo15)/(a51 + a53 + 
                  c51 + c53));
                atoms.popoi[3] = 0.;
                atoms.popoi[2] = (realnum)(atoms.popoi[4]*(a53 + c53)/(a32 + c32));
                atoms.popoi[1] = (realnum)((atoms.popoi[2]*(a32 + c32) + abundoi*
                  c12)/(a21 + c21));

                atoms.popoi[0] = (realnum)abundoi;
                /*  write(QQ,'('' OILEVL resets this to simple pop'',1P,6E10.2)')
                 *  1   popoi */
        }

        /* this is total cooling due to model atom, can be neg (heating) */
        *coloi = 
           (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 + 
           (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 + 
           (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 + 
           (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 + 
           (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 + 
           (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 + 
           (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 + 
           (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 + 
           (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13;

        return;
}

---