/home66/gary/public_html/cloudy/c08_branch/source/cool_nitr.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 */
/*CoolNitr evaluate total cooling due to nitrogen */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "phycon.h"
#include "ligbar.h"
#include "thermal.h"
#include "lines_service.h"
#include "embesq.h"
#include "atoms.h"
#include "cooling.h"
#include "nitro.h"
/* xNI_coll_stren: Author: Terry Yun 2006 */
/* This program interpolate the collision strength of NI 
 * We are only interested in first 5 levels, therefor there are 10 transitions.
 * Used the Excel to find the polynomial fit, eqn: CS = a*t+b*t^1.5+c*t^0.5
 * The range of temp is 0 - 50,000 K */
/* #define PRINT */
static const int N1_SIZE = 10;

static double aNI[N1_SIZE] = {2.755e-5,4.123e-5,7.536e-6,1.486e-5,4.516e-5,-2.935e-6,4.000e-6,3.751e-6,-2.176e-6,1.024e-5};
static double bNI[N1_SIZE] = {-8.150e-8,-1.220e-7,-2.226e-8,-4.390e-8,-1.130e-7,8.000e-9,-1.1447e-8,-1.061e-8,5.610e-9,-3.227e-8};
static double cNI[N1_SIZE] = {2.140e-4,3.272e-4,-4.944e-6,3.473e-6,-8.772e-4,1.654e-3,1.675e-3,1.123e-3,3.867e-3,3.376e-4};

static double NI[6][6][3]; /* coefficients a b c into array */

/*xNI_coll_stren interpolate the collision strength of NI  */
STATIC double xNI_coll_stren(int init, int final)
{
        /* incorporating the N0 collision strengths give in
         *>>refer       N1      cs      Tayal, S.S., 2006, ApJS, 163, 207 */

        int i, j;
        double CS; /* collision strength */
        int index = 0;
        double temp_max = 50e3;
        double temp_min = 0;
        double temp = 0;

        /* assigning array initially to zero */ 
        for(i=0; i<6; i++)
        {
                for(j=0; j<6; j++)
                {
                        NI[i][j][0] = 0; /* coeff a */
                        NI[i][j][1] = 0; /* coeff b */
                        NI[i][j][2] = 0; /* coeff c */
                }
        }

        /* reading coeffs into 3D array */   
        /* The index is in physics scale */
        for(i=1; i<6; i++)
        {
                for(j=i+1; j<6; j++)
                {
                        NI[i][j][0] = aNI[index];
                        NI[i][j][1] = bNI[index];
                        NI[i][j][2] = cNI[index];
                        index++;
                }
        }

#       ifdef PRINT
        /* physical index goes from 1 to 5 */
        for(i=1; i<6; i++)
        {
                for(j=i+1; j<6; j++)
                {
                        printf("NI %i%i a:%e ", i, j, NI[i][j][0]);
                        printf("%i%i b:%e\n", i, j, NI[i][j][1]);
                        /* printf("%i%i c:%lf\n", i, j, NI[i][j][2]); */
                }
        }
#       endif


        /* Invalid entries returns '-1':the initial indices are smaller than the final indices */
        if(init >= final)
        {
                CS = -1;
        }

        /* Invalid returns '-1': the indices are greater than 5 or smaller than 0 */
        else if(init < 1 || init > 5 || final < 1 || final > 5)
        {
                CS = -1;
        }

        else
        {
                temp = MAX2(temp, temp_min); /* return lager temp */
                temp = MIN2(temp, temp_max); /* return smaller temp */
                /* eqn: CS = a*t+b*t^1.5+c*t^0.5 */
                CS = NI[init][final][0]*phycon.te + NI[init][final][1]*phycon.te32 + NI[init][final][2]*phycon.sqrte;
        }

        return CS;
}


void CoolNitr(void)
{
        realnum p2;
        double a21, 
          a31, 
          a32, 
          cs, 
          cs2s2p, 
          cs2s3p, 
          cs21, 
          cs31, 
          cs32, 
          p3,
          pump_rate;
        long int i;

        double a12, 
          a13, 
          a14, 
          a15, 
          a23, 
          a24, 
          a25, 
          a34, 
          a35, 
          a45, 
          cs12, 
          cs13, 
          cs14, 
          cs15, 
          cs23, 
          cs24, 
          cs25, 
          cs34, 
          cs35, 
          cs45;

        double pop[5];
        static double gAr4[5]={4.,6.,4.,2.,4.};
        static double exAr4[4]={19224.464,8.713,9605.74,0.386};

        static long int *ipN1Pump=NULL,
                nN1Pump=0 , lgFirst_N1=true;

        static bool lgFirst_N3=true;
        static long int *ipN3Pump=NULL,
                nN3Pump=0;

/* #define PRINT */

#       ifdef PRINT
        FILE *ofp;
        double pop1, pop2,pop3,pop4,pop5;
        ofp = fopen("c:\\projects\\cloudy\\trunk\\mydata\\pop.out", "w");
        if(ofp == NULL)
        {
                printf("Can't open the file\n");
                cdEXIT( EXIT_FAILURE );
        }
#       endif

        DEBUG_ENTRY( "CoolNitr()" );

        /* >>chng 00 dec 11, add made-up collision strengths 
         * These have made-up collision strengths  */
#       if 0
        PutCS(1.,&TauLines[ipT671]);
        atom_level2(&TauLines[ipT671]);

        PutCS(1.,&TauLines[ipT315]);
        atom_level2(&TauLines[ipT315]);

        PutCS(1.,&TauLines[ipT333]);
        atom_level2(&TauLines[ipT333]);

        PutCS(1.,&TauLines[ipT324]);
        atom_level2(&TauLines[ipT324]);

        PutCS(1.,&TauLines[ipT374g]);
        atom_level2(&TauLines[ipT374g]);

        PutCS(1.,&TauLines[ipT374x]);
        atom_level2(&TauLines[ipT374x]);
#       endif

        /* find pump rate for [N I] 5199, do one time initialization if first 
         * call and level 2 lines are on - NB lgFirst is used in two places
         * in this routine - here we use it, and the last place where it is 
         * used it will be set false */
        if( lgFirst_N1 && nWindLine && dense.lgElmtOn[ipNITROGEN] )
        {
                lgFirst_N1 = false;
                nN1Pump = 0;
                for( i=0; i<nWindLine; ++i )
                {
                        /* don't test on nelem==ipNITROGEN since lines on physics, not C, scale */
                        if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==1  )
                        {
                                ++nN1Pump;
                        }
                }
                if( nN1Pump<0 )
                        TotalInsanity();
                else if( nN1Pump > 0 )
                        /* create the space - can't malloc 0 bytes */
                        ipN1Pump = (long *)MALLOC((unsigned)(nN1Pump)*sizeof(long) );
                nN1Pump = 0;
                for( i=0; i<nWindLine; ++i )
                {
                        /* don't test on nelem==ipNITROGEN since lines on physics, not C, scale */
                        if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==1  )
                        {
#                               if      0
                                DumpLine( &TauLine2[i] );
#                               endif
                                ipN1Pump[nN1Pump] = i;
                                ++nN1Pump;
                        }
                }
        }
        else
                /* level 2 lines are not enabled */
                nN1Pump = 0;

        /* now sum pump rates - will be used in prt statement together
         * with quench rate derive here to estimate continuum pumping of
         [N I] 5199 */
        nitro.pump_rate_N1 = 0.;
        for( i=0; i<nN1Pump; ++i )
        {
                nitro.pump_rate_N1 += TauLine2[ipN1Pump[i]].Emis->pump;
#               if      0
                fprintf(ioQQQ,"DEBUG N %li %.3e %.3e\n",
                        i,
                        TauLine2[ipN1Pump[i]].WLAng , TauLine2[ipN1Pump[i]].pump );
#               endif
        }

        /* >>chng 01 sep 02, go back to old values */
        /* nitrogen I 5200, cs and A from 
         * collision strengths from 
         * >>referold   n1      cs      Dopita, M.A., Mason, D.J., Robb, W.D. 1976, ApJ, 207, 102*/
        cs21 = 1.32e-4*phycon.te/(phycon.te10*phycon.te01);
        cs31 = 3.60e-5*phycon.te/phycon.te10*phycon.te02;
        cs32 = 5.23e-4*phycon.te70*phycon.te10/phycon.te02;

        /* N.B. following cs also appears in nitrogen ionization balance */
        /* >>chng revise collision strengths to
         * >>refer      n1      cs      Tayal, S.S., 2000, ADNDT, 76, 191 
        cs21 = 5.25e-13*phycon.tesqrd*phycon.te70*phycon.te03;
        cs31 = 3.79e-5*phycon.te70*phycon.te03;
        cs32 = 3.14e-13*phycon.tesqrd*phycon.te*phycon.te20*phycon.te07/phycon.te003;*/

        /* following does not have photoionization from excited state
         * since already counted in photo balance equation */
        /* trans prob from 
         * >>refer      n1      as      Butler, K., & Zeippen, C.J. 1984, A&A, 141, 274
         */
        a21 = 1.28e-5;
        a31 = 5.31e-3;
        a32 = 6.81e-2;
        /* >>chng 01 sep 02, include atoms.d5200r term here too */

        /* trans prob from 
         * >>refer      n1      as      Butler, K., & Zeippen, C.J. 1984, A&A, 141, 274
         */
        cs12 = xNI_coll_stren(1, 2);
        a12 = 7.170e-6;

        cs13 = xNI_coll_stren(1, 3);
        a13 = 3.994e-5;

        cs14 = xNI_coll_stren(1, 4);
        a14 = 5.397e-3;

        cs15 = xNI_coll_stren(1,5);
        a15 = 1.329e-2;

        /* this is fraction of 2D excitations that emit a photon in the 5200 multiplet
         * this is used for collisional quenching in prt lines routine */
        nitro.quench_5200 = (a12+a13) / ((a12+a13) + (cs12 + cs13) * dense.cdsqte );

        /* collision strengths and transition probabilities for the lowest 5
         * levels of N^0 */
        cs23 = xNI_coll_stren(2, 3);
        a23 = 2.478e-8;

        cs24 = xNI_coll_stren(2, 4);
        a24 = 3.135e-2;

        cs25 = xNI_coll_stren(2, 5);
        a25 = 5.693e-2;

        cs34 = xNI_coll_stren(3, 4);
        a34 = 4.92e-2;

        cs35 = xNI_coll_stren(3, 5);
        a35 = 2.709e-2;

        cs45 = xNI_coll_stren(4, 5);
        a45 = 1e-15;

        atom_pop5(gAr4,exAr4,cs12,cs13,cs14,cs15,cs23,cs24,cs25,cs34,cs35,
                  cs45,a12,a13,a14,a15,a23,a24,a25,a34,a35,a45,pop,dense.xIonDense[ipNITROGEN][0]);

#       ifdef PRINT
        pop1 = pop[0]/dense.xIonDense[ipNITROGEN][0];
        pop2 = pop[1]/dense.xIonDense[ipNITROGEN][0];
        pop3 = pop[2]/dense.xIonDense[ipNITROGEN][0];
        pop4 = pop[3]/dense.xIonDense[ipNITROGEN][0];
        pop5 = pop[4]/dense.xIonDense[ipNITROGEN][0];
        fprintf(ofp, "%lf %lf %lf %lf %lf\n", pop1, pop2, pop3, pop4, pop5);

#       endif

        nitro.xN5200 = pop[1]*a12*3.818e-12;
        nitro.xN5198 = pop[2]*a13*3.821e-12;
        nitro.xN3467 = pop[3]*a14*5.729e-12;
        nitro.xN3466 = pop[4]*a15*5.729e-12;
        nitro.xN10398 = pop[3]*a24*1.910e-12;
        nitro.xN10397 = pop[4]*a25*1.910e-12;
        nitro.xN10408 = pop[3]*a34*1.908e-12;
        nitro.xN10407 = pop[4]*a35*1.908e-12;

        p3 = atom_pop3(4.,10.,6.,cs21,cs31,cs32,a21,a31,a32,2.769e4,1.38e4,
          &p2,dense.xIonDense[ipNITROGEN][0],atoms.d5200r,0.,0.);
        /* save population */
        atoms.p2nit = p2/SDIV(dense.xIonDense[ipNITROGEN][0]);

        nitro.c5200 = p2*a21*3.83e-12;
        CoolAdd("N  1",5200,nitro.c5200);
        thermal.dCooldT += nitro.c5200*(2.769e4*thermal.tsq1 + thermal.halfte);

        nitro.c10400 = p3*a32*1.91e-12;
        CoolAdd("N  1",10400,nitro.c10400);

        nitro.c3466 = p3*a31*5.74e-12;
        CoolAdd("N  1",3466,nitro.c3466);
        thermal.dCooldT += (nitro.c10400 + nitro.c3466)*(4.15e4*
          thermal.tsq1 + thermal.halfte);

        /* N I 1200, cs from trans probability */
        PutCS(4.1,&TauLines[ipT1200]);
        atom_level2(&TauLines[ipT1200]);

        /* N II 1084, cs from trans prob */
        PutCS(5.5,&TauLines[ipT1085]);
        atom_level2(&TauLines[ipT1085]);

        /* [N II] coll data from 
         * >>referold   n2      cs      Stafford, R.P., Bell, K.L, Hibbert, A. & Wijesundera, W.P.,
         * >>rereroldcon 1994, MNRAS 268, 816,
         * at 10,000K (v weak T dep)
         * >>chng 00 dec 11, to Lennon & Burke
         * >>refer      n2      cs      Lennon, D.J., & Burke, V.M., 1994, A&AS, 103, 273-277
         * transit prob from 
         * >>refer      n2      as      Nussbaumer, H., & Rusca, C. 1979, A&A, 72, 129
         * 
         * >>chng 00 dec 11, to Lennon & Burke cs */
        a21 = 3.65e-3;
        a31 = 0.0316;
        a32 = 1.17;
        /* >>chng 02 may 02, put in option to switch between Lennon Burke and Stafford et al. ,
         * default state of this var is true */
#       define USE_LENNON_BURKE 1
#       if( USE_LENNON_BURKE )
                cs21 = 2.64;
                cs31 = 0.293;
                cs32 = 0.834;
#       else
                /* Bob Rubin could not find the Stafford et al. numbers - it flag ever
                 * set false will not compile */
                cs21 = 3.02;
                cs31 = 0.372;
                cs32 = 0.505;
#       endif

        /* POP3(   G1,G2,G3,   O12,  O13,  O23,  A21, A31, A32,*/
        p3 = atom_pop3(  9.,5.,1.,  cs21, cs31, cs32, a21, a31, a32 ,
        /*   E12,E23,P2,ABUND,GAM2) */
                21955.,24982.,&p2,dense.xIonDense[ipNITROGEN][1],0.,0.,0.);

        nitro.c5755 = p3*a32*3.46e-12;

        nitro.c6584 = p2*a21*3.03e-12;
        thermal.dCooldT += nitro.c6584*(2.2e4*thermal.tsq1 - thermal.halfte);
        CoolAdd("N  2",5755,nitro.c5755);
        CoolAdd("N  2",6584,nitro.c6584);

        /* nitro.xN2_A3_tot is fraction of excitations that produce a photon 
         * and represents the correction for collisiona deexcitation */
        nitro.xN2_A3_tot = (a31+a32) /(a31+a32 + (cs31+cs32)/1.*dense.cdsqte );
        ASSERT( nitro.xN2_A3_tot <= 1. );

        /* N II fine structure lines, */
        /* >>chng 00 dec 11, to Lennon & Burke CS */
#       if( USE_LENNON_BURKE )
        cs21 = 0.408;
        cs32 = 1.12;
        cs31 = 0.272;
#       else
        /* Bob Rubin could not find the Stafford et al. numbers - it flag ever
         * set false will not compile */
        cs21 = 0.429;
        cs32 = 1.13;
        cs31 = 0.265;
#       endif

        PutCS(cs21,&TauLines[ipT205]);
        PutCS(cs32,&TauLines[ipT122]);
        PutCS(cs31,&TauDummy);
        atom_level3(&TauLines[ipT205],&TauLines[ipT122],&TauDummy);

        /* N II 2140, data 
         * >>refer      n2      cs      Stafford, R.P., Bell, K.L, Hibbert, A. & Wijesundera, W.P. 1994,
         * >>rerercon MNRAS, 268, 816
         * >>refer      n2      cs      Lennon, D.J., & Burke, V.M., 1994, A&AS, 103, 273-277
         * A from 
         * >>refer      n2      as      Brage, T., Hibbert, A., Leckrone, D.S. 1997, ApJ, 478, 423 */
#       if( USE_LENNON_BURKE )
                cs21 = 1.19;
#       else
                /* Bob Rubin could not find the Stafford et al. numbers - it flag ever
                 * set false will not compile */
                cs21 = 1.15;
#       endif

        PutCS(cs21,&TauLines[ipT2140]);
        atom_level2(&TauLines[ipT2140]);

        /* N III 989.8, cs from 
         * >>refer      N3      cs      Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
        PutCS(7.12,&TauLines[ipT990]);
        atom_level2(&TauLines[ipT990]);

        /* 57 micron N III, A=
         * >>refer      n3      as      Froese Fischer, C. 1983, J.Phys. B, 16, 157
         * collision strength from 
         * >>refer      n3      cs      Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
        cs = MIN2(1.90,0.2291*phycon.te10*phycon.te10);
        PutCS(cs,&TauLines[ipT57]);


        /* one time initialization if first call, and level 2 lines are on */
        if( lgFirst_N3 && nWindLine && dense.lgElmtOn[ipNITROGEN] )
        {
                lgFirst_N3 = false;
                nN3Pump = 0;
                for( i=0; i<nWindLine; ++i )
                {
                        /* don't test on nelem==ipIRON since lines on physics, not C, scale */
                        if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==3  )
                        {
                                ++nN3Pump;
                        }
                }
                if( nN3Pump<0 )
                        TotalInsanity();
                else if( nN3Pump > 0 )
                        /* create the space - can't malloc 0 bytes */
                        ipN3Pump = (long *)MALLOC((unsigned)(nN3Pump)*sizeof(long) );
                nN3Pump = 0;
                for( i=0; i<nWindLine; ++i )
                {
                        /* don't test on nelem==ipIRON since lines on physics, not C, scale */
                        if( TauLine2[i].Hi->nelem ==7 && TauLine2[i].Hi->IonStg==3  )
                        {
#                               if      0
                                DumpLine( &TauLine2[i] );
#                               endif
                                ipN3Pump[nN3Pump] = i;
                                ++nN3Pump;
                        }
                }
        }
        else
                /* level 2 lines are not enabled */
                nN3Pump = 0;

        /* now sum pump rates */
        pump_rate = 0.;
        for( i=0; i<nN3Pump; ++i )
        {
                pump_rate += TauLine2[ipN3Pump[i]].Emis->pump;
#               if      0
                fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
                        i,
                        TauLine2[ipN3Pump[i]].WLAng , TauLine2[ipN3Pump[i]].pump );
#               endif
        }

        /* N III] N 3, N  3, 1765 multiplet */
        /*atom_level2(&TauLines[ipT57]);*/
        /*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
        AtomSeqBoron(&TauLines[ipT57], 
          &TauLines[ipN3_1749], 
          &TauLines[ipN3_1747], 
          &TauLines[ipN3_1754], 
          &TauLines[ipN3_1752], 
          &TauLines[ipN3_1751], 
          0.201 , 1.088 , 0.668 , 2.044 , pump_rate,"N  3");
        /*fprintf(ioQQQ," n4 %.3e\n",  (
                TauLines[ipN3_1749].xIntensity + 
                TauLines[ipN3_1747].xIntensity + 
                TauLines[ipN3_1754].xIntensity+ 
                TauLines[ipN3_1752].xIntensity+ 
                TauLines[ipN3_1751].xIntensity) / dense.xIonDense[ipNITROGEN][2] );*/

        /* N IV 1486, N 4, N  4,collisions within 3P just guess
         * cs to ground from 
         * >>refer      n4      cs      Ramsbottom, C.A., Berrington, K.A., Hibbert, A., Bell, K.L. 1994,
         * >>refercon Physica Scripta, 50, 246 */
        if( phycon.te > 1.584e4 )
        {
                cs = 21.346/(phycon.te10*phycon.te10*phycon.te10*phycon.te02);
        }
        else
        {
                cs = 75.221/(phycon.sqrte/phycon.te03/phycon.te02);
        }
        /* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
        cs = MAX2(0.01,cs);
        PutCS(cs,&TauLines[ipT1486]);
        AtomSeqBeryllium(.9,.9,3.0,&TauLines[ipT1486],.0115);
        embesq.em1486 = (realnum)(atoms.PopLevels[3]*0.0115*1.34e-11);
        /*fprintf(ioQQQ," n4 %.3e\n", (TauLines[ipT1486].xIntensity + embesq.em1486 ) / dense.xIonDense[ipNITROGEN][3] );*/

        /* N IV 765, cs from 
         * >>refer      n4      cs      Ramsbottom, C.A., Berrington, K.A., Hibbert, A., Bell, K.L. 1994,
         * >>refercon Physica Scripta, 50, 246 */
        /* >>refer      n4      as      Flemming, J., Brage, T., Bell, K.L., Vaeck, N., Hibbert, A., 
         * >>refercon   Godefroid, M., & Froese Fischer, C., 1995, ApJ, 455, 758*/
        cs = MIN2(4.0,1.864*phycon.te03*phycon.te03);
        PutCS(cs,&TauLines[ipT765]);
        atom_level2(&TauLines[ipT765]);

        /* N V 1240
         * >>refer      n5      cs      Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
        ligbar(7,&TauLines[ipT1239],&TauLines[ipT209],&cs2s2p,&cs2s3p);
        PutCS(cs2s2p,&TauLines[ipT1239]);
        PutCS(cs2s2p*0.5,&TauLines[ipT1243]);
        PutCS(1.0,&TauDummy);
        atom_level3(&TauLines[ipT1243],&TauDummy,&TauLines[ipT1239]);
        /*fprintf(ioQQQ," n5 %.3e\n", (TauLines[ipT1243].xIntensity + TauLines[ipT1239].xIntensity ) / dense.xIonDense[ipNITROGEN][4] );*/

        /* N V 209 */
        PutCS(cs2s3p,&TauLines[ipT209]);
        atom_level2(&TauLines[ipT209]);
        return;
}

---