root/branches/newmole/source/cont_gammas.cpp

/* 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 */
/* GammaBn evaluate photoionization rate for single shell with induced recomb */
/* GammaBnPL 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 */
/* GammaPL evaluate photoionization rate for power law photo cross section */
/* 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.  
 */

/*>>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)
{
        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.;
                thermal.HeatNet = 0.;
                thermal.HeatLowEnr = 0.;
                thermal.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;

        thermal.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;
        thermal.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;
                thermal.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 */
        thermal.HeatNet = MAX2(0.,thermal.HeatNet - g*thresh);

        /* we will save this heat - the part that does not secondary ionize */
        thermal.HeatLowEnr = thermal.HeatNet;

        /* now do part where secondary ionization is possible */
        thermal.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;
                thermal.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 */
        thermal.HeatHiEnr = thermal.HeatHiEnr - GamHi*thresh;

        /* add corrected high energy heat to total */
        thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;

        /* final product is heating in erg */
        thermal.HeatNet *= EN1RYD;
        thermal.HeatHiEnr *= EN1RYD;
        thermal.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( thermal.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))*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 punch 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)
{
        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.;
                thermal.HeatNet = 0.;
                thermal.HeatHiEnr = 0.;
                thermal.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.;
        thermal.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;
        thermal.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;
                thermal.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 */
        /*thermal.HeatNet += -gamk_v*eauger;*/
        ASSERT( thermal.HeatNet >= 0. );
        thermal.HeatNet = MAX2(0.,thermal.HeatNet - gamk_v*eauger);

        /* this is the total low-energy heating, in ryd, that cannot secondary ionize */
        thermal.HeatLowEnr = thermal.HeatNet;

        /* now do part where secondary ionization possible */
        thermal.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;
                thermal.HeatHiEnr += phisig*rfield.anu[i];
        }

        /* add on the high energy part */
        gamk_v += GamHi;
        /* correct for work function */
        thermal.HeatHiEnr -= GamHi*eauger;

        /* net heating include high energy heat, with correction for
         * secondary ionization */
        thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;

        /* final product is heating in erg */
        thermal.HeatNet *= EN1RYD;
        thermal.HeatLowEnr *= EN1RYD;
        thermal.HeatHiEnr *= EN1RYD;

        ASSERT( gamk_v >= 0. );
        ASSERT( thermal.HeatNet>= 0. );
        return( gamk_v );
}

#if 0
//this is never used
/*GammaPL evaluate photoionization rate for power law photoionization cross section */

double GammaPL(
        /* this is prin quantum number for level */
        long int n, 
        /* nelem = 0 for H */
        long int nelem)
{
        long int i, 
          ilo, 
          iup, 
          limit, 
          ipLoEnr, 
          ipHiEnr;
        realnum GamHi, 
          csThresh, 
          gamPL_v, 
          phisig;

        DEBUG_ENTRY( "GammaPL()" );

        /* evaluate photoioinzation rate and photoelectric heating
         * returns photoionization rate as GammaPL, heating is H
         * n and nelem and principal quantum number and charge
         */

        /* validate the input */
        ASSERT( nelem >= 0);
        ASSERT( nelem < LIMELM);
        ASSERT( n >= 0);

        /* get pointer to photo threshold */
        ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n];
        ASSERT( ipLoEnr>0 );

        if( ipLoEnr >= rfield.nflux )
        {
                gamPL_v = 0.;
                thermal.HeatNet = 0.;
                thermal.HeatHiEnr = 0.;
                thermal.HeatLowEnr = 0.;
                return( gamPL_v );
        }

        ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
        iup = MIN2(ipHiEnr,rfield.nflux);

        csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1]/ POW2(nelem+1.f);

        /* low energies where secondary ionization cannot occur
         * will do threshold point, ipLoEnr, later */
        gamPL_v = 0.;
        thermal.HeatNet = 0.;

        /* 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 ) /rfield.anu3[ipLoEnr-1];
        gamPL_v += phisig;

        /* 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]/rfield.anu3[i];
                gamPL_v += phisig;
                thermal.HeatNet += phisig*rfield.anu[i];
        }

        /* convert photo rate into proper units with CS at threshold */
        gamPL_v *= csThresh;

        /* correct heating for CS at threshold and work function */
        thermal.HeatNet *= csThresh;
        thermal.HeatNet += -gamPL_v*rfield.anu[ipLoEnr-1];
        thermal.HeatLowEnr = thermal.HeatNet;

        /* now do part where secondary ionization possible */
        thermal.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]/rfield.anu3[i];
                GamHi += phisig;
                thermal.HeatHiEnr += phisig*rfield.anu[i];
        }

        /* add on the high energy part */
        gamPL_v += GamHi*csThresh;

        /* correct heating for work function and convert to true cross section */
        thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*rfield.anu[ipLoEnr-1])*csThresh;

        /* net heating includes high energy heating corrected for secondary ionizations */
        thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;

        /* final product is heating in erg */
        thermal.HeatNet *= EN1RYD;
        thermal.HeatHiEnr *= EN1RYD;
        thermal.HeatLowEnr *= EN1RYD;

        ASSERT( gamPL_v>= 0.);
        return( gamPL_v );
}
#endif

/*GammaBnPL evaluate hydrogenic photoionization rate for single level with induced recomb */

double GammaBnPL(long int n, 
  long int nelem, /* 0 for H, etc */
  double *ainduc, 
  double *rcool)
{
        long int i, 
          ilo, 
          iup, 
          limit, 
          ipLoEnr, 
          ipHiEnr;
        double GamHi, 
          bnfunPL_v, 
          csThresh, 
          emin, 
          g, 
          phisig, 
          RateInducRec, 
          RateInducRecCool, 
          prod, 
          thresh;

        DEBUG_ENTRY( "GammaBnPL()" );

        /* validate the incoming charge */
        ASSERT( nelem >= 0);
        ASSERT( nelem < LIMELM );

        /* special version of GAMFUN for arbitrary threshold, and induc rec
         * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */

        /* these are the upper and lower energy bounds, which we know from iso-sequence */
        ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n];
        ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];

        ilo = ipLoEnr;
        iup = MIN2(ipHiEnr,rfield.nflux);
        if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
        {
                bnfunPL_v = 0.;
                thermal.HeatNet = 0.;
                thermal.HeatHiEnr = 0.;
                thermal.HeatLowEnr = 0.;
                *ainduc = 0.;
                *rcool = 0.;
                return( bnfunPL_v );
        }

        /* >>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 */
        thresh = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][n];
        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;

        thermal.HeatNet = 0.;
        g = 0.;
        RateInducRec = 0.;
        RateInducRecCool = 0.;

        /*csThresh = t_ADfA::Inst().sth(n)*POW3(thresh) / POW2(nelem+1.f);*/
        csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1] / POW2(nelem+1.f);

        /* do first point without otscon, which may have diffuse cont */
        /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
        phisig = (rfield.flux[0][ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+ 
                rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly)/
                rfield.anu3[ipLoEnr-1];

        g += phisig;
        thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1];

        /* integral part of induced rec rate */
        prod = phisig*rfield.ContBoltz[ipLoEnr-1];
        RateInducRec += prod;
        /* induced rec cooling */
        RateInducRecCool += prod*(rfield.anu[ipLoEnr-1] - emin);

        limit = MIN2(iup,secondaries.ipSecIon-1);
        for( i=ipLoEnr; i < limit; i++ )
        {
                phisig = rfield.SummedCon[i]/rfield.anu3[i];

                g += phisig;
                thermal.HeatNet += 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);
        }

        /* heating in Rydbergs, no secondary ionization */
        thermal.HeatNet = MAX2(0.,(thermal.HeatNet - g*thresh)*csThresh);
        thermal.HeatLowEnr = thermal.HeatNet;

        /* now do part where secondary ionization is possible */
        thermal.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]/rfield.anu3[i];
                GamHi += phisig;
                thermal.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);
        }
        RateInducRec *= csThresh;
        RateInducRecCool *= csThresh;

        /* this is total photo rate, low and high energy parts */
        bnfunPL_v = (g + GamHi)*csThresh;

        /* heating in Rydbergs, uncorrected for secondary ionization */
        thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*thresh)*csThresh;

        /* net heating includes high energy heat corrected for secondary ionization */
        thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary;

        /* final product is heating in erg */
        thermal.HeatNet *= EN1RYD;
        thermal.HeatHiEnr *= EN1RYD;
        thermal.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( bnfunPL_v>= 0.);
        return( bnfunPL_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;
}