/home66/gary/public_html/cloudy/c08_branch/source/hydro_recom.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 */
/*hlike_radrecomb_from_cross_section - integrate Milne relation */
/*RecomInt - Integral in Milne relation.  Called by qg32.       */
/*H_cross_section returns cross section (cm^-2), 
 * given EgammaRyd, the photon energy in Ryd,
 * ipLevel, the index of the level, 0 is ground,
 * nelem is charge, equal to 0 for Hydrogen */
#include "cddefines.h" 
#include "physconst.h" 
#include "hydro_bauman.h"
#include "iso.h"
#include "helike.h"
#include "dense.h"
#include "opacity.h"
#include "atmdat.h"

static double RecomInt(double EE);

static double kTRyd,EthRyd; 
static long int ipLev,globalZ;

/*H_cross_section returns cross section (cm^-2), 
 * given EgammaRyd, the photon energy in Ryd,
 * ipLevel, the index of the level, 0 is ground,
 * nelem is charge, equal to 0 for Hydrogen */
double H_cross_section( double EgammaRyd , long ipLev , long nelem )
{
        double cs;
        double rel_photon_energy;
        long ipISO = ipH_LIKE;

        /* make sure this routine not called within collapsed high levels */
        ASSERT( ipLev < iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem] );
        ASSERT( ipLev > ipH1s);

        EthRyd = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipLev];

        /* >>chng 02 apr 24, more protection against calling with too small an energy */
        /* evaluating H-like photo cs at He energies, may be below threshold -
         * prevent this from happening */
        rel_photon_energy = EgammaRyd / EthRyd;
        rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );

        cs = H_photo_cs(rel_photon_energy , N_(ipLev), L_(ipLev), nelem + 1 );

        ASSERT( cs > 0. && cs < 1.E-8 );

        return cs;
}

/*hlike_radrecomb_from_cross_section calculates radiative recombination coefficients. */
double hlike_radrecomb_from_cross_section(double temp, long nelem, long ipLo)
{
        double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive;
        double change[5] = {0.,0.,0.,0.,0.};
        long ipISO = ipH_LIKE;

        if( ipLo == 0 )
                return t_ADfA::Inst().H_rad_rec(nelem+1,ipLo,temp);

        EthRyd = iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];

        /* Factors outside integral in Milne relation.  */
        b = MILNE_CONST * StatesElem[ipISO][nelem][ipLo].g * pow(temp,-1.5) / 2.;
        /* kT in Rydbergs.      */
        kTRyd = temp / TE1RYD;
        globalZ = nelem;
        ipLev = ipLo;

        /* Begin integration.   */
        /* First define characteristic step */
        E1 = EthRyd;
        step = MIN2( 0.125*kTRyd, 0.5*E1 );
        E2 = E1 + step;
        /* Perform initial integration, from threshold to threshold + step.     */
        RecomIntegral = qg32( E1, E2, RecomInt);
        /* Repeat the integration, adding each new result to the total, 
         * except that the step size is doubled every time, since values away from 
         * threshold tend to fall off more slowly.      */
        do
        {
                OldRecomIntegral = RecomIntegral;
                E1 = E2;
                step *= 1.25;
                E2 = E1 + step;
                RecomIntegral += qg32( E1, E2, RecomInt);
                change[4] = change[3];
                change[3] = change[2];
                change[2] = change[1];
                change[1] = change[0];
                change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral;
                TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4];
        /* Continue integration until the upper limit exceeds 100kTRyd, an arbitrary
         * point at which the integrand component exp(electron energy/kT) is very small,
         * making neglible cross-sections at photon energies beyond that point,
         * OR when the last five steps resulted in less than a 1 percent change.        */
        } while ( ((E2-EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) );

        /* Calculate recombination coefficient. */
        alpha = b * RecomIntegral;

        alpha = MAX2( alpha, SMALLDOUBLE );

        return alpha;
}

/*RecomInt, used in comput milne relation for he rec - the energy is photon Rydbergs.   */
static double RecomInt(double ERyd)
{
        double x1, temp;

        /* Milne relation integrand     */
        x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - EthRyd ) / kTRyd);
        temp = H_cross_section( ERyd , ipLev, globalZ );
        x1 *= temp;

        return x1;
}

---