cont_gammas.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/* GammaBn evaluate photoionization rate for single shell with induced recomb */
/* GammaPrt special version of gamma function to print strong contributors */
/* GammaK evaluate photoionization rate for single shell */
/* GammaPrtRate print photo rates for all shells of a ion and element */
/*GammaPrtShells for the element nelem and ion, print total photo rate, subshells,
 * and call GamaPrt for important subshells */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "thermal.h"
#include "secondaries.h"
#include "opacity.h"
#include "rfield.h"
#include "ionbal.h"
#include "atmdat.h"
#include "heavy.h"
#include "gammas.h"
 
/*
 * these are the routines that evaluate the photoionization rates, gammas,
 * throughout cloudy.  a considerable amount of time is spent in these routines,
 * so they should be compiled at the highest possible efficientcy.  
 */
 
/* removed these two unused routines around r3500 (late October, 2009). */
/* GammaBnPL evaluate photoionization rate for single shell with induced recomb */
/* GammaPL evaluate photoionization rate for power law photo cross section */
 
 
/*>>chng 99 apr 16, ConInterOut had not been included in the threshold point for any
 * of these routines.  added it.  also moved thresholds above loop for a few */
 
double GammaBn(
  /* index for low energy */
  long int ipLoEnr, 
  /* index for high energy */
  long int ipHiEnr, 
  /* offset within the opacity stack */
  long int ipOpac, 
  /* threshold (Ryd) for arbitrary photoionization */
  double thresh, 
  /* induced rec rate */
  double *ainduc, 
  /* rec cooling */
  double *rcool,
  t_phoHeat *photoHeat)
{
        long int i, 
          ilo, 
          iup, 
          limit;
        double GamHi, 
          bnfun_v, 
          emin, 
          g, 
          phisig, 
          RateInducRec, 
          RateInducRecCool, 
          prod;
 
        DEBUG_ENTRY( "GammaBn()" );
 
/**********************************************************************
 *
 * special version of GAMFUN for arbitrary threshold, and induc rec
 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool
 *
 **********************************************************************/
 
        /* special version of GAMFUN for arbitrary threshold, and induc rec
         * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */
        if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
        {
                bnfun_v = 0.;
                photoHeat->HeatNet = 0.;
                photoHeat->HeatLowEnr = 0.;
                photoHeat->HeatHiEnr = 0.;
                *ainduc = 0.;
                *rcool = 0.;
                return( bnfun_v );
        }
 
        ASSERT( ipLoEnr >= 0 && ipHiEnr >= 0 );
 
        /* >>chng 96 oct 25, first photo elec energy may have been negative
         * possible for first any point to be below threshold due to
         * finite resolution of continuum mesh */
        emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
        /* >>chng 99 jun 08 heating systematically too small, change to correct
         * threshold, and protect first cell */
        emin = thresh;
 
        photoHeat->HeatNet = 0.;
        g = 0.;
        RateInducRec = 0.;
        RateInducRecCool = 0.;
 
        /* do first point without otscon, which may have diffuse cont,
         * extra minus one after ipOpac is correct since not present in i = */
        i = ipLoEnr;
        /* >>chng 99 apr 16, add ConInterOut */
        /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
        phisig = (rfield.flux[0][i-1] + rfield.otslin[i-1] + 
                rfield.ConInterOut[i-1]*rfield.lgOutOnly)*
                opac.OpacStack[i-ipLoEnr+ipOpac-1];
 
        g += phisig;
        photoHeat->HeatNet += phisig*rfield.anu[i-1];
 
        /* integral part of induced rec rate */
        prod = phisig*rfield.ContBoltz[i-1];
        RateInducRec += prod;
        /* incuded rec cooling */
        RateInducRecCool += prod*(rfield.anu[i-1] - emin);
 
        iup = MIN2(ipHiEnr,rfield.nflux);
        limit = MIN2(iup,secondaries.ipSecIon-1);
 
        /* these is no -1 after the ipLoEnr in the following, due to c off by one counting */
        for( i=ipLoEnr; i < limit; i++ )
        {
                /* SummedCon defined in RT_OTS_Update, 
                 * includes flux, otscon, otslin, ConInterOut, outlin */
                phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
 
                g += phisig;
                photoHeat->HeatNet += phisig*rfield.anu[i];
 
                /* phi-sigma times exp(-hnu/kT) */
                prod = phisig*rfield.ContBoltz[i];
                /* induced recombination rate */
                RateInducRec += prod;
                /* incuded rec cooling */
                RateInducRecCool += prod*(rfield.anu[i] - emin);
        }
 
        /* heating in Rydbergs, no secondary ionization */
        /* >>chng 02 mar 31, added MAX2 due to roundoff error
         * creating slightly negative number, caught downstream as insanity */
        photoHeat->HeatNet = MAX2(0.,photoHeat->HeatNet - g*thresh);
 
        /* we will save this heat - the part that does not secondary ionize */
        photoHeat->HeatLowEnr = photoHeat->HeatNet;
 
        /* now do part where secondary ionization is possible */
        photoHeat->HeatHiEnr = 0.;
        GamHi = 0.;
 
        /* make sure we don't count twice when ipSecIon = ipLoEnr */
        ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
        for( i=ilo-1; i < iup; i++ )
        {
                phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
 
                GamHi += phisig;
                photoHeat->HeatHiEnr += phisig*rfield.anu[i];
 
                /* integral part of induced rec rate */
                prod = phisig*rfield.ContBoltz[i];
                RateInducRec += prod;
                /* incuded rec cooling */
                RateInducRecCool += prod*(rfield.anu[i] - emin);
        }
 
        /* this is total photo rate, low and high energy parts */
        bnfun_v = g + GamHi;
 
        /* heating in Rydbergs, uncorrected for secondary ionization */
        photoHeat->HeatHiEnr = photoHeat->HeatHiEnr - GamHi*thresh;
 
        /* add corrected high energy heat to total */
        photoHeat->HeatNet += photoHeat->HeatHiEnr*secondaries.HeatEfficPrimary;
 
        /* final product is heating in erg */
        photoHeat->HeatNet *= EN1RYD;
        photoHeat->HeatHiEnr *= EN1RYD;
        photoHeat->HeatLowEnr *= EN1RYD;
 
        /* this is an option to turn off induced processes */
        if( rfield.lgInducProcess )
        {
                *rcool = RateInducRecCool*EN1RYD;
                *ainduc = RateInducRec;
        }
        else
        {
                /* the "no induced" command was given */
                *rcool = 0.;
                *ainduc = 0.;
        }
 
        ASSERT( bnfun_v >= 0. );
        ASSERT( photoHeat->HeatNet>= 0. );
        return( bnfun_v );
}
 
/*GammaPrtShells for the element nelem and ion, print total photo rate, subshells,
 * and call GamaPrt for important subshells */
void GammaPrtShells( long nelem , long ion )
{
        double sum;
        long int ns;
 
        DEBUG_ENTRY( "GammaPrtShells()" );
 
        fprintf(ioQQQ," GammaPrtShells nz\t%.2f \t%.2li %.2li ", 
                fnzone ,
                nelem,
                ion
                );
        sum = 0;
        for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
        {
                sum += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
        }
        fprintf(ioQQQ,"\ttot\t%.2e", sum);
        /* loop over all shells for this ion */
        for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
        {
                fprintf(ioQQQ,"\t%.2e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
#               if 0
                /* always reevaluate the outer shell, and all shells if lgRedoStatic is set */
                if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
                {
                        /* option to redo the rates only on occasion */
                        iplow = opac.ipElement[nelem][ion][ns][0];
                        iphi = opac.ipElement[nelem][ion][ns][1];
                        ipop = opac.ipElement[nelem][ion][ns][2];
 
                        /* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0
                                * with "no photoionization" command */
                        ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 
                                GammaK(iplow,iphi,
                                ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0),
                                &phoHeat)*ionbal.lgPhotoIoniz_On;
                }
#               endif
        }
        fprintf(ioQQQ,"\n");
        return;
}
 
/**********************************************************************
 *
 * GammaPrt special version of gamma function to print strong contributors
 * this only prints info - does not update heating rates properly since
 * this is only a diagnostic
 *
 **********************************************************************/
 
void GammaPrt(
          /* first three par are identical to GammaK */
          long int ipLoEnr, 
          long int ipHiEnr, 
          long int ipOpac, 
          /* io unit we will write to */
          FILE * ioFILE, 
          /* total photo rate from previous call to GammaK */
          double total, 
          /* we will print contributors that are more than this rate */
          double threshold)
{
        long int i, 
          j, 
          k;
        double flxcor, 
          phisig;
 
        DEBUG_ENTRY( "GammaPrt()" );
 
        /* special special version of GAMFUN to save step-by-step results */
        /* >>chng 99 apr 02, test for lower greater than higher (when shell
         * does not exist) added.  This caused incorrect photo rates for
         * non-existant shells */
        if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
        { 
                return;
        }
 
        fprintf( ioFILE, " GammaPrt %.2f from ",fnzone);
        fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1]));
        fprintf( ioFILE, " to ");
        fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1]));
        fprintf( ioFILE, "R rates >");
        fprintf( ioFILE,PrintEfmt("%9.2e",threshold));
        fprintf( ioFILE, " of total=");
        fprintf( ioFILE,PrintEfmt("%9.2e",total));
        fprintf( ioFILE, " (frac inc, otslin, otscon, ConInterOut, outlin ConOTS_local_OTS_rate ) chL, C\n");
 
        if( threshold <= 0. || total <= 0. )
        { 
                return;
        }
 
        k = ipLoEnr;
        j = MIN2(ipHiEnr,rfield.nflux);
 
        /* do theshold special, do not pick up otscon */
        i = k-1;
 
        /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
        phisig = (rfield.flux[0][i] + rfield.otslin[i]+ rfield.ConInterOut[i]*rfield.lgOutOnly)*
                opac.OpacStack[i-k+ipOpac];
        if( phisig > threshold || phisig < 0.)
        {
                flxcor = rfield.flux[0][i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly;
 
                /* this really is array index on C scale */
                fprintf( ioFILE, "[%5ld]" , i );
                fprintf( ioFILE, PrintEfmt("%9.2e",rfield.anu[i]));
                fprintf( ioFILE, PrintEfmt("%9.2e",phisig/total));
                fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n", 
                  rfield.flux[0][i]/SDIV(flxcor), 
                  rfield.otslin[i]/SDIV(flxcor), 
                  /* otscon will appear below, but is not counted here, so do not print it (deceiving) */
                  0./SDIV(flxcor), 
                  rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor), 
                  (rfield.outlin[0][i]+rfield.outlin_noplot[i])/SDIV(flxcor), 
                  rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor), 
                  rfield.chLineLabel[i], 
                  rfield.chContLabel[i], 
                  opac.OpacStack[i-k+ipOpac]);
        }
 
        for( i=k; i < j; i++ )
        {
                phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac];
                if( phisig > threshold || phisig < 0.)
                {
                        /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
                        flxcor = rfield.flux[0][i] + rfield.otslin[i] + rfield.otscon[i] + 
                          rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i]*rfield.lgOutOnly;
 
                        fprintf( ioFILE, "%5ld", i );
                        fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
                        fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
                        fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n", 
                          rfield.flux[0][i]/SDIV(flxcor), 
                          rfield.otslin[i]/SDIV(flxcor), 
                          rfield.otscon[i]/SDIV(flxcor), 
                          rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor), 
                          (rfield.outlin[0][i] + rfield.outlin_noplot[i])/SDIV(flxcor), 
                          rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor), 
                          rfield.chLineLabel[i], 
                          rfield.chContLabel[i], 
                          opac.OpacStack[i-k+ipOpac]);
                }
        }
        return;
}
 
/*GammaK evaluate photoionization rate for single shell */
 
/* this routine is a major pace setter for the code
 * carefully check anything put into the following do loop */
 
double GammaK(
        /* ipLoEnr and ipHiEnr are pointers within frequency array to upper and
         * lower bounds of evaluation */
        long int ipLoEnr, 
        long int ipHiEnr, 
        /* ipOpac is offset to map onto OPSV*/
        long int ipOpac, 
        /* yield1 is fraction of ionizations that emit 1 electron only,
         * only used for energy balance */
        double yield1,
        t_phoHeat *photoHeat)
{
        long int i, 
          ilo, 
          ipOffset, 
          iup, 
          limit;
        double GamHi, 
          eauger; 
 
        double gamk_v ,
                phisig; 
 
        DEBUG_ENTRY( "GammaK()" );
 
        /* evaluate photoioinzation rate and photoelectric heating
         * returns photoionization rate as GAMK, heating is H    */
 
        /* >>chng 99 apr 02, test for lower greater than higher (when shell
         * does not exist) added.  This caused incorrect photo rates for
         * non-existant shells */
        if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr)
        {
                gamk_v = 0.;
                photoHeat->HeatNet = 0.;
                photoHeat->HeatHiEnr = 0.;
                photoHeat->HeatLowEnr = 0.;
                return( gamk_v );
        }
 
        iup = MIN2(ipHiEnr,rfield.nflux);
 
        /* anu(i) is threshold, assume each auger electron has this energy
         * less threshold energy, IP lost in initial photoionizaiton
         * yield1 is the fraction of photos that emit 1 electron */
        eauger = rfield.anu[ipLoEnr-1]*yield1;
 
        /* low energies where secondary ionization cannot occur
         * will do threshold point, ipLoEnr, later */
        gamk_v = 0.;
        photoHeat->HeatNet = 0.;
 
        /* set up total offset for pointer addition not in loop */
        ipOffset = ipOpac - ipLoEnr;
 
        /* >>>chng 99 apr 16, this had followed the loop below, moved here to 
         * be like rest of gamma functions */
        /* first do the threshold point
         * do not include otscon, which may contain diffuse continuum */
        /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
        phisig = (rfield.flux[0][ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+ 
                rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly) * opac.OpacStack[ipOpac-1];
        gamk_v += phisig;
        photoHeat->HeatNet += phisig*rfield.anu[ipLoEnr-1];
 
        /* this loop only executed for energies than cannot sec ioniz */
        limit = MIN2(iup,secondaries.ipSecIon-1);
        for( i=ipLoEnr; i < limit; i++ )
        {
                phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
                gamk_v += phisig;
                photoHeat->HeatNet += phisig*rfield.anu[i];
        }
 
        /* correct heating for work function */
        /* >>chng 02 apr 10, from first to sec, due to neg heat in blister.in */
        /* *photoHeat->HeatNet += -gamk_v*eauger;*/
        ASSERT( photoHeat->HeatNet >= 0. );
        photoHeat->HeatNet = MAX2(0.,photoHeat->HeatNet - gamk_v*eauger);
 
        /* this is the total low-energy heating, in ryd, that cannot secondary ionize */
        photoHeat->HeatLowEnr = photoHeat->HeatNet;
 
        /* now do part where secondary ionization possible */
        photoHeat->HeatHiEnr = 0.;
        GamHi = 0.;
        /* make sure we don't count twice when ipSecIon = ipLoEnr */
        ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
        for( i=ilo-1; i < iup; i++ )
        {
                /* produce of flux and cross section */
                phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
                GamHi += phisig;
                photoHeat->HeatHiEnr += phisig*rfield.anu[i];
        }
 
        /* add on the high energy part */
        gamk_v += GamHi;
        /* correct for work function */
        photoHeat->HeatHiEnr -= GamHi*eauger;
 
        /* net heating include high energy heat, with correction for
         * secondary ionization */
        photoHeat->HeatNet += photoHeat->HeatHiEnr*secondaries.HeatEfficPrimary;
 
        /* final product is heating in erg */
        photoHeat->HeatNet *= EN1RYD;
        photoHeat->HeatLowEnr *= EN1RYD;
        photoHeat->HeatHiEnr *= EN1RYD;
 
        ASSERT( gamk_v >= 0. );
        ASSERT( photoHeat->HeatNet>= 0. );
        return( gamk_v );
}
 
/* GammaPrtRate will print resulting rates for ion and element */
void GammaPrtRate(
        /* io unit we will write to */
        FILE * ioFILE, 
        /* stage of ionization on C scale, 0 for atom */
        long int ion ,
        /* 0 for H, etc */
        long int nelem,
        /* true - then print photo sources for each shell */
        bool lgPRT )
{
        long int nshell , ns;
 
        DEBUG_ENTRY( "GammaPrtRate()" );
 
        /* number of shells for this species */
        nshell = Heavy.nsShells[nelem][ion];
 
        /* now print subshell photo rate */
        fprintf(ioFILE , "GammaPrtRate: %li %li",ion , nelem );
        for( ns=nshell-1; ns>=0; --ns )
        {
                fprintf(ioFILE , " %.2e" ,  ionbal.PhotoRate_Shell[nelem][ion][ns][0]);
 
                /* option to print individual contributors to each shell */
                if( lgPRT )
                {
                        fprintf(ioFILE , "\n");
                        GammaPrt(
                                /* first three par are identical to GammaK */
                                opac.ipElement[nelem][ion][ns][0], 
                                opac.ipElement[nelem][ion][ns][1], 
                                opac.ipElement[nelem][ion][ns][2], 
                                /* io unit we will write to */
                                ioFILE, 
                                /* total photo rate from previous call to GammaK */
                                ionbal.PhotoRate_Shell[nelem][ion][ns][0], 
                                /* we will print contributors that are more than this rate */
                                ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.05);
                }
        }
        fprintf(ioFILE , "\n");
        return;
}