/home66/gary/public_html/cloudy/c08_branch/source/mole_co_atom.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 */
/*CO_PopsEmisCool evaluate rotation levels populations, emission, and cooling */
/*cdCO_colden return column density in H2, negative -1 if cannot find state,
 * header is cddrive */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "thermal.h"
#include "hmi.h"
#include "radius.h"
#include "atoms.h"
#include "phycon.h"
#include "rt.h"
#include "lines_service.h"
#include "cddrive.h"
#include "mole.h"
/*lint -e778 constant expression evaluatess to 0 in operation '-' */
/*=================================================================*/

/* will be used to save CO column densities */
static double *col12 , *col13;

/* evaluate CO rotation cooling */
void CO_PopsEmisCool(
                transition ** Rotate , 
                long int nRotate ,
                realnum abundan, 
                const char * chLabel ,
                realnum * Cooling ,
                realnum * dCoolingDT  )
{

        /* will need to MALLOC space for these but only on first call */
        static double 
                **AulEscp ,
                **col_str ,
                **AulDest, 
                /* AulPump[low][high] is rate (s^-1) from lower to upper level */
                **AulPump,
                **CollRate,
                *pops,
                *create,
                *destroy,
                *depart,
                /* statistical weight */
                *stat ,
                /* excitation energies in kelvin */
                *excit;

        /*static long int **ipdest;*/

        static bool lgFirst=true;
        long int i,
        j,
        ilo , 
        ihi;
        static long int nUsed;
        bool lgDeBug;
        int nNegPop;
        bool lgZeroPop;
        double rot_cooling , dCoolDT;
        static long int ndimMalloced = 0;

        DEBUG_ENTRY( "CO_PopsEmisCool()" );

        if( lgFirst )
        {
                /* will never do this again */
                lgFirst = false;
                /* remember how much space we malloced in case ever called with more needed */
                ndimMalloced = nRotate;
                /* allocate the 1D arrays*/
                excit = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
                stat = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
                pops = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
                create = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
                destroy = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
                depart = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
                /* create space for the 2D arrays */
                AulPump = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
                CollRate = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
                AulDest = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
                AulEscp = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
                col_str = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
                for( i=0; i<(nRotate+1); ++i )
                {
                        AulPump[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
                        CollRate[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
                        AulDest[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
                        AulEscp[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
                        col_str[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
                }
        }
        /* this is test for call with too many rotation levels to handle - logic needs
         * for largest rotor to be called first */
        if( nRotate > ndimMalloced )
        {
                fprintf(ioQQQ," CO_PopsEmisCool has been called with the number of rotor levels greater than space allocated.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* all elements are used, and must be set to zero if zero */
        for( i=0; i < (nRotate+1); i++ )
        {
                create[i] = 0.;
                destroy[i] = 0.;
                for( j=0; j < (nRotate+1); j++ )
                {
                        col_str[j][i] = 0.;
                        AulEscp[j][i] = 0.;
                        AulDest[j][i] = 0.;
                        AulPump[j][i] = 0.;
                }
        }

        /* the statistical weights of the levels */
        for( j=0; j < nRotate; j++ )
        {
                /* statistical weights for each level */
                stat[j] = (*Rotate)[j].Lo->g;
        }
        /* this is the highest level, which is one more than the highest line */
        stat[nRotate] = (*Rotate)[nRotate-1].Hi->g;

        /* set up the excitation potentials of each level relative to ground -
         * the struc saves the energy of the line only */
        excit[0] = 0.;
        for( j=1; j < nRotate; j++ )
        {
                /* excitation energy of each level relative to ground, in K */
                excit[j] = excit[j-1] + (*Rotate)[j-1].EnergyK;
        }

        /* this is the highest level, which is one more than the highest line */
        excit[nRotate] = excit[nRotate-1] + (*Rotate)[nRotate-1].EnergyK;

        nUsed = nRotate;

        /* this determines the largest molecule that can be inverted at this
         * temperature.  Need Boltzmann factors to be positive for all levels
         * nUsed is the index of the highest level, so the number of levels (passed
         * to the solver) is nUsed+1 
         * excit[nUsed]/phycon.te cannot be much larger than 20, or matrix inversion will fail */
        /* >>chng 03 sep 18, allow to go to 1, a two level atom */
        /*while ( (excit[nUsed] > phycon.te*20.) && (nUsed > 1) )*/
        /* >>chng 03 oct 03, chng 20 to 25 */
        /*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/
        /* >>chng 03 oct 03, keep at least 5 levels */
        /*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/
        while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 5) )
        {
                --nUsed;
        }

        for( j=0; j < nRotate; j++ )
        {
                /*data[j][j+1] = (*Rotate)[j].Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc);*/
                AulEscp[j+1][j] = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc);
                AulDest[j+1][j] = (*Rotate)[j].Emis->Aul*(*Rotate)[j].Emis->Pdest;
                /* next 2 not not used by levelN since flag passed saying to use CollRate instead */
                (*Rotate)[j].Coll.cs = 1.;
                /*data[j+1][j] = (*Rotate)[j].Coll.cs*hmi.Hmolec[ipMH2g]/dense.eden;*/
                col_str[j+1][j] = (*Rotate)[j].Coll.cs*hmi.H2_total/dense.eden;
                /* photon pumping rate */
                AulPump[j][j+1] = (*Rotate)[j].Emis->pump;
                /* the continuum indices are on the f, not c, scale, and will be passed to 
                 * atom_levelN, which works on f, not c, scale */
                /*ipdest[j][j+1] = (*Rotate)[j].ipCont;*/
        }

        /* next loop for collisional transitions that have delta J >1 */
        for( ilo=0; ilo < nRotate-1; ilo++ )
        {
                /* need to have upper limit to to nRotate because 
                 * number ov levels is 1 gt number of lines*/
                for( ihi=ilo+2; ihi <= nRotate; ihi++ )
                {
                        /* this is not used by levelN since flag passed saying to use CollRate instead */
                        /*data[ihi][ilo] = 1. *hmi.Hmolec[ipMH2g]/dense.eden;*/
                        /*data[ihi][ilo] = 1. *hmi.H2_total/dense.eden;*/
                        col_str[ihi][ilo] = 1. *hmi.H2_total/dense.eden;
                        /* these are escape and dest rates, which are zero for a rigid rotator */
                        /*data[ilo][ihi] = 0;*/
                        AulEscp[ihi][ilo] = 0;
                        AulDest[ihi][ilo] = 0;
                }
        }

        /* now evaluate the H2 collision rates */
        /* recall one more level than lines */
        for( ilo=0; ilo < nRotate; ilo++ )
        {
                /* need to have upper limit to to nRotate because 
                 * number of levels is 1 gt number of lines*/
                for( ihi=ilo+1; ihi <= MIN2(nRotate,ilo+5); ihi++ )
                {
                        /* >>refer      CO-H2   collision       de Jong, T., Chu, S-I., & Dalgarno, A. 1975, ApJ, 199, 69 */
                        double a[5]={1.66, 2.80, 1.19, 1.00, 1.30};
                        double b[5]={1.67, 1.47, 1.85, 1.55, 2.24};
                        /* >>refer      CO-He   collision       McKee, C.F., Storey, J.W.V., Watson, D.M., & Green, S., 1982, ApJ, 259, 647 */
                        /* this reference gives He-CO collisions,
                         * their table 1 says He collisions are 1.37 slower than H2 collisions*/
                        /* >>chng 03 sep 19, from ground H2 to total H2 */
                        double collid = hmi.H2_total + dense.xIonDense[ipHELIUM][0]/1.37;
                        /* de Jong et al. explicitly consider temperatures as low as 20K,
                         * don't extrapolate significantly lower than this */
                        double TeUsed = MAX2(10., phycon.te );

                        /* first do deexcitation rate, equation 17 of deJong et al. */
                        CollRate[ihi][ilo] = a[ihi-ilo-1]*1.e-10*(*Rotate)[ilo].Lo->g/(*Rotate)[ilo].Hi->g*
                                (1.+(excit[ihi]-excit[ilo])/TeUsed) *
                                sexp( b[ihi-ilo-1]*sqrt((excit[ihi]-excit[ilo])/TeUsed) )*collid;
                        /* this is mainly so that rest of code gets valid cs in case crit dense printed */
                        if( ihi == ilo+1 )
                        {
                                /* save downward collision rate */
                                (*Rotate)[ilo].Coll.ColUL = (realnum)CollRate[ihi][ilo];
                                /* convert into fake electron collision strength */
                                LineConvRate2CS( &(*Rotate)[ilo] , (*Rotate)[ilo].Coll.ColUL);
                        }

                        /* now get excitation rate */
                        CollRate[ilo][ihi] = CollRate[ihi][ilo]*
                                sexp( (excit[ihi]-excit[ilo])/phycon.te )*
                                (*Rotate)[ilo].Hi->g / (*Rotate)[ilo].Lo->g;

                        /* debug print statement */
                        /*fprintf(ioQQQ," %li %li %.2e %.2e \n",ilo, ihi, 
                                CollRate[ihi][ilo]/hmi.H2_total  , CollRate[ilo][ihi]/hmi.H2_total );*/
                }
                /* finish off with zeros */
                for( ihi=ilo+6; ihi <= nRotate; ihi++ )
                {
                        CollRate[ihi][ilo] = 0.;
                        CollRate[ilo][ihi] = 0.;
                }
        }

        lgDeBug = false;

        atom_levelN(
                /* number of levels is number of lines plus one */
                /* set nUsed so that CO is evaluated even at very low temperatures */
                nUsed+1, /*nRotate+1,*/
                abundan,
                stat,
                excit,
                'K',
                pops,
                depart,
                &AulEscp,
                &col_str,
                &AulDest,
                &AulPump,
                &CollRate,
                create,
                destroy,
                /* say that we have evaluated the collision rates already */
                true,
                /*&ipdest,*/
                &rot_cooling,
                &dCoolDT,
                chLabel,
                /* nNegPop positive if negative pops occured, negative if too cold */
                &nNegPop,
                &lgZeroPop,
                lgDeBug );/* option to print stuff - set to true for debug printout */

        /* put cooling into place where we can use it later */
        *Cooling = (realnum)rot_cooling;

        /* >>chng 03 sep 18, CO rot cooling temp deriv is too small, and temp
         * changes too big - incr deriv by 5x 
        *dCoolingDT = (realnum)(dCoolDT);*/
        *dCoolingDT = (realnum)(rot_cooling/phycon.te);

#       if 0
        if( lgMainCO )
                fprintf(ioQQQ,"COcool\t%i\t%.2f\t%e\t%e\t%e\t%e\t%e\n",
                nUsed,
                fnzone,
                phycon.te,
                *Cooling/abundan,
                *dCoolingDT ,
                *Cooling,
                thermal.htot );
#       endif

        /* zero out higher populations for case where full CO levels not done */
        for( i=nUsed+1; i<=nRotate; ++i )
        {
                pops[i] = 0.;
                depart[i] = 0.;
        }

        /* can only define first LIMLEVELN elements, the vector's length */
        for( i=0; i< MIN2(LIMLEVELN,nRotate); ++i )
        {
                atoms.PopLevels[i] = pops[i];
                atoms.DepLTELevels[i] = depart[i];
        }

        if( nNegPop > 0 )
        {
                fprintf(ioQQQ,"CO_PopsEmisCool called atom_levelN which returned negative populations.\n");
        }

        /* now establish information that is passed out to rest of code's infrastructure */
        for( j=0; j<nRotate; ++j )
        {
                double EnrLU, EnrUL;
                /* lower upper populations, stim emission correct population */
                (*Rotate)[j].Lo->Pop = pops[j];
                (*Rotate)[j].Hi->Pop = pops[j+1];
                (*Rotate)[j].Emis->PopOpc = (pops[j] - pops[j+1]*(*Rotate)[j].Lo->g/(*Rotate)[j].Hi->g);

                /* number of photons in the line */
                (*Rotate)[j].Emis->phots = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc)*pops[j+1];

#               if 0
                /* >>chng 03 oct 04, moved to RT_OTS */
                /* local Emis.ots rates, added to line ots array */
                (*Rotate)[j].Emis->ots = (*Rotate)[j].Aul*(*Rotate)[j].Emis->Pdest*(realnum)(*Rotate)[j].Hi->Pop;
                RT_OTS_AddLine( (*Rotate)[j].Emis->ots , (*Rotate)[j].ipCont);
#               endif

                /* the intensity in the line */
                (*Rotate)[j].Emis->xIntensity = (*Rotate)[j].Emis->phots*(*Rotate)[j].EnergyErg;

                /* ratio of collisional to total excitation */
                (*Rotate)[j].Emis->ColOvTot = (realnum)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].Emis->pump) );

                /* two cases - collisionally excited (usual case) or 
                 * radiatively excited - in which case line can be a heat source
                 * following are correct heat exchange, if will mix to get correct deriv */
                EnrLU = (*Rotate)[j].Lo->Pop*CollRate[j][j+1]*(*Rotate)[j].EnergyErg;
                EnrUL = (*Rotate)[j].Hi->Pop*CollRate[j+1][j]*(*Rotate)[j].EnergyErg;
                /* energy exchange due to this level
                 * net cooling due to excit minus part of de-excit */
                (*Rotate)[j].Coll.cool = EnrLU - EnrUL*(*Rotate)[j].Emis->ColOvTot;
                /* net heating is remainder */
                (*Rotate)[j].Coll.heat = EnrUL*(1.f - (*Rotate)[j].Emis->ColOvTot);
                /* do not add to cooling, since done with evaluated cooling from atom_levelN */
                /*CoolAdd( chLabel, (long)t10->WLAng , t10->cool);*/
        }

        /* generate flag if co cooling important and highest level is fainter
         * than second highest level */
        if( rot_cooling / thermal.ctot > 0.1 && 
                 (*Rotate)[nUsed-1].Emis->xIntensity > (*Rotate)[nUsed-2].Emis->xIntensity &&
                 /*>>chng 03 oct 03, add check whether molecule has been trimed down 
                  * due to small Boltzmann factor */ 
                 /* >>chng 03 oct 20, following had; after ), so CoolCaped always true */
                 nUsed == nRotate )
        {
                co.lgCOCoolCaped = true;
        }

        return;
}

/*CO_Colden maintain H2 column densities within X */
void CO_Colden( const char *chLabel )
{
        long int iRot;

        DEBUG_ENTRY( "CO_Colden()" );

        if( strcmp(chLabel,"ZERO") == 0 )
        {
                static bool lgFIRST=true;
                if( lgFIRST )
                {
                        lgFIRST = false;
                        col12 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) );
                        col13 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) );
                }
                /* the column density (cm-2) of ortho and para H2 */
                /* zero out formation rates and column densites */
                for( iRot=0; iRot<=nCORotate; ++iRot )
                {
                        /* zero it out */
                        col12[iRot] = 0.;
                        col13[iRot] = 0.;
                }
        }
        else if( strcmp(chLabel,"ADD ") == 0 )
        {
                /*  add together column densities */
                for( iRot=0; iRot<nCORotate; ++iRot )
                {
                        /* zero it out */
                        col12[iRot] += C12O16Rotate[iRot].Lo->Pop*radius.drad_x_fillfac;
                        col13[iRot] += C13O16Rotate[iRot].Lo->Pop*radius.drad_x_fillfac;
                }
        }

        /* we will not print column densities so skip that - if not print then we have a problem */
        else if( strcmp(chLabel,"PRIN") != 0 )
        {
                fprintf( ioQQQ, " CO_Colden does not understand the label %s\n", 
                                 chLabel );
                cdEXIT(EXIT_FAILURE);
        }
}

/*cdCO_colden return column density in H2, negative -1 if cannot find state,
 * header is cddrive */
double cdCO_colden( long isotope , long iRot )
{

        /* make sure incoming parameters are ok */
        if( isotope !=12 && isotope != 13 )
        {
                fprintf(ioQQQ," cdCO_colden can only do 12CO and 13CO\n");
                return -1.;
        }
        if( iRot < 0 || iRot > nCORotate )
        {
                fprintf(ioQQQ," cdCO_colden - rotation quantum number must be 0 or greater, and less than %li\n", 
                                nCORotate);
                return -1.;
        }

        /* the incoming parameters are fully qualified - return the column density */
        if( isotope == 12 )
        {
                return col12[iRot];
        }
        else
        {
                return col13[iRot];
        }
}

/*CO_OTS - add CO lines to ots fields */
void CO_OTS( void )
{

        long int j;

        DEBUG_ENTRY( "CO_OTS()" );

        /* add all CO lines */
        for( j=0; j<nCORotate; ++j )
        {
                C12O16Rotate[j].Emis->ots = C12O16Rotate[j].Emis->Aul*C12O16Rotate[j].Emis->Pdest*
                        C12O16Rotate[j].Hi->Pop;
                RT_OTS_AddLine( C12O16Rotate[j].Emis->ots , C12O16Rotate[j].ipCont);

                C13O16Rotate[j].Emis->ots = C13O16Rotate[j].Emis->Aul*C13O16Rotate[j].Emis->Pdest*
                        C13O16Rotate[j].Hi->Pop;
                RT_OTS_AddLine( C13O16Rotate[j].Emis->ots , C13O16Rotate[j].ipCont);
        }

        return;
}

---