/home66/gary/public_html/cloudy/c08_branch/source/cool_dima.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 */
/*CoolDima compute cooling due to level 2 lines */
/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "dense.h"
#include "phycon.h"
#include "conv.h"
#include "thermal.h"
#include "opacity.h"
#include "lines_service.h"
#include "rfield.h"
#include "mewecoef.h"
#include "atoms.h"
#include "cooling.h"

/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
STATIC double ColStrGBar(transition * t , realnum cs1 );

void CoolDima(void)
{
        long int i, 
          ion,
          nelem;
        double cs;
        double OTSLevel2,
                sumots,
                CoolLevel2;
        static double TeEvalCS = -1.;
        static long int nzoneEval=0;

        DEBUG_ENTRY( "CoolDima()" );

        /* >>chng 03 nov 29, add this option to skip rest of routine */
        /* no level2 command sets nWindLine to -1 */
        if( nWindLine<0 )
                return;

        /* only force evaluation of cooling when significant, or first call
         * in this zone, or not into first zone */
        if( nzone != nzoneEval || conv.lgLevel2_Cool_Imp || !nzone )
        {
                nzoneEval = nzone;

                /* check whether we need to reevaluate the collision strenghts.
                 * Do so if large change in temperature  or any stage of ionization has been lowered. */
                if( (conv.lgSearch || conv.lgIonStageTrimed || 
                        fabs(phycon.te-TeEvalCS)/phycon.te > 0.05) )
                {
                        for( i=0; i < nWindLine; i++ )
                        {
                                if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
                                {
                                        /* only evaluate cs if positive abundance */
                                        ion = TauLine2[i].Hi->IonStg;
                                        nelem = TauLine2[i].Hi->nelem;

                                        if( dense.xIonDense[nelem-1][ion-1] > 0. )
                                        {
                                                /* now generate the collision strength */
                                                cs = ColStrGBar(&TauLine2[i] , cs1_flag_lev2[i] );
                                        }
                                        else
                                        {
                                                cs = 1.;
                                        }
                                        /* now put the cs into the line array */
                                        PutCS(cs,&TauLine2[i] );
                                }
                        }
                        TeEvalCS = phycon.te;
                }

                /* now always update the cooling since this adds ots flux */
                for( i=0; i < nWindLine; i++ )
                {
                        /* we need to call this even if nothing present since then sets zero
                         * ignore all lines with negative atomic number
                         * >>chng 97 aug 30, get rid of test for non-pos ipLnNelem
                         * pointer tells atom_level2 to use existing WineEtc labels */
                        /* only call non-hydrogenic, non-he-like lines */
                         /* >>chng 02 sug 11, do not include he-like in sum */
                        if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
                        {
                                atom_level2(&TauLine2[i] );
                        }
                }
        }
        else
        {
                /* now use old cooling and ots flux */
                for( i=0; i < nWindLine; i++ )
                {
                        if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
                        {
                                /* recall that these should be trivial cooling, ots rates */
                                CoolAdd( "    ", TauLine2[i].WLAng , TauLine2[i].Coll.cool);
                                /*>>chng 03 oct 04, move to RT_OTS */
                                /*RT_OTS_AddLine( TauLine2[i].ots , TauLine2[i].ipCont);*/
                        }
                }
        }

        /* find ratio of level 1 to level 2 ots to see how important level 2 is.
         * if level 2 ots not important then will update dest probs only first time
         * in each zone.  If significant then continuously update */
        OTSLevel2 = 0.;
        CoolLevel2 = 0.;
        if( nzone > 1 )
        {
                for( i=0; i < nWindLine; i++ )
                {
                        long int ip = TauLine2[i].ipCont-1;
                        CoolLevel2 += TauLine2[i].Coll.cool;
                        if( opac.opacity_abs[ip] > SMALLFLOAT )
                        {
                                OTSLevel2 += TauLine2[i].Emis->ots / ( (realnum)opac.opacity_abs[ip]);
                                ASSERT( fabs(OTSLevel2) < 1e37 );
                        }
                }

                /* now sum all ots rates */
                sumots = 0.;
                for( i=0; i < rfield.nflux; i++ )
                {
                        sumots += rfield.otslin[i];
                }
                /* is the ots contribution significant ? */
                if( OTSLevel2/SDIV(sumots) > 1e-4 )
                {
                        /* true if level 2 lines were contributors to the ots rates, set in dimacool */
                        conv.lgLevel2_OTS_Imp = true;
                }
                else
                {
                        /* true if level 2 lines were contributors to the ots rates, set in dimacool */
                        conv.lgLevel2_OTS_Imp = false;
                }

                /* set flags saying if important */
                if( CoolLevel2/SDIV(thermal.ctot) > 1e-4 )
                {
                        /* true if level 2 lines were contributors to the cooling, set in dimacool */
                        conv.lgLevel2_Cool_Imp = true;
                }
                else
                {
                        /* true if level 2 lines were contributors to the cooling, set in dimacool */
                        conv.lgLevel2_Cool_Imp = false;
                }

                {
                        /* debugging code for Lya problems */
                        /*@-redef@*/
                        enum {DEBUG_LOC=false};
                        /*@+redef@*/
                        if( DEBUG_LOC )
                        {
                                fprintf(ioQQQ,"%li\t%.2e\t%.2e\n", 
                                        nzone, OTSLevel2/SDIV(sumots)  , CoolLevel2/SDIV(thermal.ctot) );
                        }
                }
        }
        else
        {
                /* true if level 2 lines were contributors to the ots rates, set in dimacool */
                conv.lgLevel2_OTS_Imp = true;
                /* true if level 2 lines were contributors to the cooling, set in dimacool */
                conv.lgLevel2_Cool_Imp = true;
        }
        return;
}

/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
STATIC double ColStrGBar(transition * t , realnum cs1 )
{
        long int i, 
          j;
        double ColStrGBar_v, 
          a, 
          b, 
          c, 
          d, 
          e1, 
          gb, 
          x, 
          y;
        double xx, 
          yy;

        DEBUG_ENTRY( "ColStrGBar()" );

        /* Calculation of the collision strengths of multiplets.
         * Neutrals are recalculated from 
         * >>refer cs   gbar    Fisher et al. (1996)
         * >>refer cs   gbar    Gaetz & Salpeter (1983, ApJS 52, 155) and 
         * >>refer cs   gbar    Mewe (1972, A&A 20, 215) 
         * fits for ions. */

        /* routine to implement g-bar data taken from
         *>>refer       cs      gbar    Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985,
         *>>refercon    A&AS, 62, 197 */

        /* zero hydrogenic lines since they are done by iso-sequence */
        if( t->Hi->nelem == t->Hi->IonStg )
        {
                ColStrGBar_v = 0.0;
                return( ColStrGBar_v );
        }

        /*was the block data linked in? */
        ASSERT( MeweCoef.g[1][0] != 0.);

        /* which type of transition is this? cs1 is the flag */

        /* >>chng 01 may 30 - cs1 < 0 means a forced collision strength */
        if( cs1 < 0. )
        {
                ColStrGBar_v = -cs1;
                return( ColStrGBar_v );
        }

        /* >>chng 99 feb 27, change to assert */
        ASSERT( cs1 >= 0.05 );

        /* excitation energy over kT */
        y = t->EnergyK/phycon.te;
        if( cs1 < 1.5 )
        {
                xx = -log10(y);

                if( cs1 < 0.5 )
                {
                        yy = (1.398813573838321 + xx*(0.02943050869177121 + xx*
                          (-0.4439783893114510 + xx*(0.2316073358577902 + xx*(0.001870493481643103 + 
                          xx*(-0.008227246351067403))))))/(1.0 + xx*(-0.6064792600526370 + 
                          xx*(0.1958559534507252 + xx*(-0.02110452007196644 + 
                          xx*(0.01348743933722316 + xx*(-0.0001944731334371711))))));
                }

                else
                {
                        yy = (1.359675968512206 + xx*(0.04636500015069853 + xx*
                          (-0.4491620298246676 + xx*(0.2498199231048967 + xx*(0.005053803073345794 + 
                          xx*(-0.01015647880244268))))))/(1.0 + xx*(-0.5904799485819767 + 
                          xx*(0.1877833737815317 + xx*(-0.01536634911179847 + 
                          xx*(0.01530712091180953 + xx*(-0.0001909176790831023))))));
                }

                ColStrGBar_v = pow((double)10,yy)*t->Emis->gf/(t->EnergyWN * WAVNRYD* 13.6);
        }
        else
        {
                i = (long int)cs1;

                if( i < 26 )
                {
                        e1 = log(1.0+1.0/y) - 0.4/POW2(y + 1.0);
                        a = MeweCoef.g[i-1][0];
                        b = MeweCoef.g[i-1][1];
                        c = MeweCoef.g[i-1][2];
                        d = MeweCoef.g[i-1][3];
                        x = (double)t->Hi->nelem - 3.0;

                        if( i == 14 )
                        {
                                a *= 1.0 - 0.5/x;
                                b = 1.0 - 0.8/x;
                                c *= 1.0 - 1.0/x;
                        }

                        else if( i == 16 )
                        {
                                a *= 1.0 - 0.9/x;
                                b *= 1.0 - 1.7/x;
                                c *= 1.0 - 2.1/x;
                        }

                        else if( i == 18 )
                        {
                                a *= 1.0 + 2.0/x;
                                b *= 1.0 - 0.7/x;
                        }

                        gb = a + (b*y - c*y*y + d)*e1 + c*y;

                        /*  ipLnRyd is exciation energy in Rydbergs */
                        ColStrGBar_v = 14.510395*t->Emis->gf*gb/(t->EnergyWN * WAVNRYD);
                        /* following i>=26 */
                }

                else
                {
                        /* 210 is the dimem of g, so [209] is largest val */
                        if( i < 210 )
                        {
                                j = (long)(MeweCoef.g[i-1][3]);
                                if( j == 1 )
                                {
                                        ColStrGBar_v = t->Lo->g*MeweCoef.g[i-1][0]*
                                          pow(phycon.te/pow(10.,(double)MeweCoef.g[i-1][2]),(double)MeweCoef.g[i-1][1]);
                                }
                                else
                                {
                                        ColStrGBar_v = t->Lo->g*MeweCoef.g[i-1][0]*
                                          sexp(MeweCoef.g[i-1][1]*(pow(10.,(double)MeweCoef.g[i-1][2])/
                                          phycon.te));
                                }
                        }

                        else
                        {
                                /* This is for AlII 1670 line only!
                                 *    ColStrGBar=0.0125*te**0.603 */
                                /* 98 dec 27, this is still in use */
                                ColStrGBar_v = 0.0125*phycon.sqrte*phycon.te10*
                                  phycon.te003;
                        }
                }
        }

        /* following to make sure that negative values not returned */
        ColStrGBar_v = MAX2(ColStrGBar_v,1e-10);
        return( ColStrGBar_v );
}

---