/home66/gary/public_html/cloudy/c08_branch/source/hydro_vs_rates.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 */
/* CS_VS80 compute thermally averaged collision strength for collisional deexcitation 
 * of hydrogenic atoms, from
 * >>refer      H1      collision       Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
 *hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients 
 * for quantum number n */
/*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for ions 
 * with all n, and Z*/
#include "cddefines.h"
#include "phycon.h"
#include "physconst.h"
#include "iso.h"
#include "hydro_vs_rates.h"
#include "lines_service.h"

STATIC double hydro_vs_coll_str( double energy );
STATIC double Therm_ave_coll_str_int_VS80( double EOverKT );

static long global_ipISO, global_nelem, global_ipHi, global_ipLo, global_Collider;
static double global_temp, global_Aul;
/* These are masses relative to the proton mass of the electron, proton, and alpha particle. */
static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};

/*
 Neither of these rates can be modified to account for impact by non-electrons because they 
 are fits to thermally-averaged rates for electron impact...it can not be unravelled to 
 expose a projectile energy that can then be scaled to account for other projectiles.
 Instead, we have included the original cross section formula (eq 14) from 
 Vriens & Smeets (1980) below.
*/

/* VS80 stands for Vriens and Smeets 1980 */
/* This routine calculates thermally-averaged collision strengths. */
double CS_VS80( long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider )
{
        double coll_str;

        global_ipISO = ipISO;
        global_nelem = nelem;
        global_ipHi = ipHi;
        global_ipLo = ipLo;
        global_temp = temp;
        global_Collider = Collider;
        global_Aul = Aul;

        /* evaluate collision strength, two options,
         * do thermal average (very slow) if set with
         * SET COLLISION STRENGTH AVERAGE command,
         * default just do point at kT
         * tests show no impact on test suite, much faster */
        if( iso.lgCollStrenThermAver )
        {
                /* do average over Maxwellian */
                coll_str =  qg32( 0.0, 1.0, Therm_ave_coll_str_int_VS80);
                coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VS80);
        }
        else
        {
                /* the argument to the function is equivalent to evaluating the function at
                 * hnu = kT */
                coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD );
        }

        ASSERT( coll_str >= 0. );
        return coll_str;
}

/* The integrand for calculating the thermal average of collision strengths */
STATIC double Therm_ave_coll_str_int_VS80( double EOverKT )
{
        double integrand;

        integrand = exp( -1.*EOverKT ) * hydro_vs_coll_str( EOverKT * EVRYD*global_temp/TE1RYD );

        return integrand;
}

/*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom, 
 * from Vriens and Smeets */
STATIC double hydro_vs_coll_str( double energy )
{
        double Apn, bp, Bpn, delta;
        double Epi, Epn;
        double gamma, p, n;
        double ryd, s;
        double cross_section, coll_str, gLo, gHi, abs_osc_str, Aul;
        long ipISO, nelem, ipHi, ipLo, Collider;

        DEBUG_ENTRY( "hydro_vs_coll_str()" );

        ipISO = global_ipISO;
        nelem = global_nelem;
        ipHi = global_ipHi;
        ipLo = global_ipLo;
        Collider = global_Collider;
        Aul = global_Aul;

        gLo = StatesElem[ipISO][nelem][ipLo].g;
        gHi = StatesElem[ipISO][nelem][ipHi].g;

        /* This comes from equations 14,15, and 16 of Vriens and Smeets. */
        /* >>refer he-like cs Vriens, L. \& Smeets, A. H. M. Phys. Rev. A, 22, 940 */ 
        /* Computes the Vriens and Smeets
         * rate coeff. for collisional dexcitation between any two levels of H.
         * valid for all nhigh, nlow
         * at end converts this excitation rate to collision strength   */

        n = (double)StatesElem[ipISO][nelem][ipHi].n;
        p = (double)StatesElem[ipISO][nelem][ipLo].n;

        ryd = EVRYD;
        s = fabs(n-p);

        Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
        Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];

        /* This is an absorption oscillator strength. */
        abs_osc_str = GetGF( Aul, Epn*RYD_INF/EVRYD, gHi)/gLo;
        Apn = 2.*ryd/Epn*abs_osc_str;

        bp = 1.4*log(p)/p - .7/p - .51/p/p + 1.16/p/p/p - 0.55/p/p/p/p;

        Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));

        delta = exp(-Bpn/Apn) - 0.4*Epn/ryd;

        gamma = ryd*(8. + 23.*POW2(s/p))/
                ( 8. + 1.1*n*s + 0.8/s/s + 0.4*(pow(n,1.5))/(pow(s,0.5))*fabs(s-1.0) );

        /* Scale the energy to get an equivalent electron energy. */
        energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];

        /* cross section in units of PI*a_o^2 */
        if( energy/2./ryd+delta <= 0 /*|| energy < Epn*/ )
        {
                cross_section = 0.;
        }
        else
        {
                cross_section = 2.*ryd/(energy + gamma)*(Apn*log(energy/2./ryd+delta) + Bpn);
                cross_section = MAX2( cross_section, 0. );
        }

        /* convert to collision strength */
        coll_str = cross_section * gLo * (energy/EVRYD);

        if( nelem==ipCARBON )
                ASSERT( coll_str >= 0. );

        return( coll_str );
}

/*hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients 
 * for quantum number n */
double hydro_vs_ioniz( long int ipISO, long int nelem, long int level )
{
        double coef, 
          denom, 
          epi, 
          t_eV;

        DEBUG_ENTRY( "hydro_vs_ioniz()" );

        ASSERT( nelem <= 1 );

        /* a function written to calculate the rate coefficients 
         * for hydrogen collisional ionizations from
         * Jason Ferguson, summer 94
         * valid for all n
         * >>refer      H1      collision       Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
         * */

        /* this is kT is eV */
        t_eV = phycon.te/EVDEGK;

        /* this is the ionization energy relative to kT, dimensionless */
        epi = (iso.xIsoLevNIonRyd[ipISO][nelem][level] * EVRYD) / t_eV;

        /* this is the denominator of equation 8 of VS80. */
        denom = pow(epi,2.33) + 4.38*pow(epi,1.72) + 1.32*epi;

        /* this is equation 8 of VS80 */
        coef = 9.56e-6 * pow(t_eV,-1.5) * StatesElem[ipISO][nelem][level].ConBoltz / denom;

        ASSERT( coef >= 0. );
        return( coef );
}

/*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for all n, and Z*/
double Hion_coll_ioniz_ratecoef(
                /* the isoelectronic sequence */
                long int ipISO ,
                /* element, >=1 since only used for ions 
                 * nelem = 1 is helium the least possible charge */
                long int nelem,
                /* principal quantum number, > 1
                 * since only used for excited states */
                long int level)
{
        long int n;
        double H, 
          HydColIon_v, 
          Rnp, 
          charge,
          chim, 
          eone, 
          etwo, 
          ethree, 
          g, 
          rate, 
          rate2, 
          /* ryd, */
          t1, 
          t2, 
          t3, 
          t4, 
          tev, 
          xn, 
          y;
        static const double arrH[4] = {1.48,3.64,5.93,8.32};
        static const double arrRnp[8] = {2.20,1.90,1.73,1.65,1.60,1.56,1.54,1.52};
        static const double arrg[10] = {0.8675,0.932,0.952,0.960,0.965,0.969,0.972,0.975,0.978,0.981};

        static double small = 0.;

        DEBUG_ENTRY( "Hion_coll_ioniz_ratecoef()" );
        /*calculate hydrogenic ionization rates for all n, and Z
         * >>refer      HI      cs      Allen 1973, Astro. Quan. for low Te.
         * >>refer      HI      cs      Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
         * */

        charge = nelem - ipISO;
        /* this routine only for ions, nelem=0 is H, nelem=1 he, etc */
        ASSERT( charge > 0);
        /* this is quantum level, not for ground state (n=1), only
         * n=2 or higher */
        n = StatesElem[ipISO][nelem][level].n;
        ASSERT( n>1 );

        if( n > 4 )
        {
                H = 2.15*n;
        }
        else
        {
                H = arrH[n-1];
        }

        if( n > 8 )
        {
                Rnp = 1.52;
        }
        else
        {
                Rnp = arrRnp[n-1];
        }

        if( n > 10 )
        {
                g = arrg[9];
        }
        else
        {
                g = arrg[n-1];
        }

        tev = EVRYD/TE1RYD*phycon.te;
        /* ryd = EVRYD; */
        xn = (double)n;
        /* >>chng 19 dec 02, use proper energy! */
        /*chim = POW2(charge+1.)*ryd/xn/xn;*/
        chim = EVRYD * iso.xIsoLevNIonRyd[ipISO][nelem][level];
        y = chim/tev;

        eone = ee1(y);
        etwo = StatesElem[ipISO][nelem][level].ConBoltz - y*eone;
        ethree = (StatesElem[ipISO][nelem][level].ConBoltz - y*etwo)/2.;

        t1 = 1/xn*eone;
        t2 = 1./3./xn*(StatesElem[ipISO][nelem][level].ConBoltz - y*ethree);
        t3 = 3.*H/xn/(3. - Rnp)*(y*etwo - 2.*y*eone + StatesElem[ipISO][nelem][level].ConBoltz);
        t4 = 3.36*y*(eone - etwo);
        rate = 7.69415e-9*(double)phycon.sqrte*9.28278e-3*powi(xn/(charge+1),4)*g*y;
        rate *= t1 - t2 + t3 + t4;
        rate = MAX2(rate,small);

        rate2 = 2.1e-8*phycon.sqrte/chim/chim*dsexp(2.302585*5040.*
          chim/(double)phycon.te);

        rate2 = MAX2(rate2,small);

        /* Take the lowest of the two, they fit nicely together... */
        /* >>chng 10 Sept 02, sometimes one of these is zero and the other is positive.
         * in that case take the bigger one.    */
        if( rate==0. || rate2==0. )
                HydColIon_v = MAX2(rate,rate2);
        else
                HydColIon_v = MIN2(rate,rate2);

        ASSERT( HydColIon_v >= 0. );
        return( HydColIon_v );
}

/*hydro_vs_deexcit compute collisional deexcitation rates for hydrogen atom, 
 * from Vriens and Smeets 1980 */
double hydro_vs_deexcit( long ipISO, long nelem, long ipHi, long ipLo, double Aul )
{
        double Anp, bn, Bnp, delta_np;
        double Eni, Enp;
        double Gamma_np, p, n, g_p, g_n;
        double ryd, s, kT_eV, rate, col_str, abs_osc_str;

        DEBUG_ENTRY( "hydro_vs_coll_str()" );

        kT_eV = EVRYD * phycon.te / TE1RYD;

        /* This comes from equations 24 of Vriens and Smeets. */
        /* >>refer he-like cs Vriens, L. \& Smeets, A. H. M. Phys. Rev. A, 22, 940 */ 
        /* Computes the Vriens and Smeets
         * rate coeff. for collisional dexcitation between any two levels of H.
         * valid for all nhigh, nlow
         * at end converts this excitation rate to collision strength   */

        n = (double)StatesElem[ipISO][nelem][ipLo].n;
        p = (double)StatesElem[ipISO][nelem][ipHi].n;

        ASSERT( n!=p );

        g_n = StatesElem[ipISO][nelem][ipLo].g;
        g_p = StatesElem[ipISO][nelem][ipHi].g;

        ryd = EVRYD;
        s = fabs(n-p);

        Enp = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
        Eni = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipHi];

        ASSERT( Enp > 0. );

        /* This is an absorption oscillator strength. */
        abs_osc_str = GetGF( Aul, Enp*RYD_INF/EVRYD, g_p)/g_n;
        Anp = 2.*ryd/Enp*abs_osc_str;

        bn = 1.4*log(n)/n - .7/n - .51/n/n + 1.16/n/n/n - 0.55/n/n/n/n;

        Bnp = 4.*ryd*ryd/p/p/p*(1./Enp/Enp + 4./3.*Eni/POW3(Enp) + bn*Eni*Eni/powi(Enp,4));

        delta_np = exp(-Bnp/Anp) + 0.1*Enp/ryd;

        Gamma_np = ryd*log(1. + n*n*n*kT_eV/ryd) * (3. + 11.*s*s/n/n) /
                ( 6. + 1.6*p*s + 0.3/s/s + 0.8*(pow(p,1.5))/(pow(s,0.5))*fabs(s-0.6) );

        if( 0.3*kT_eV/ryd+delta_np <= 0 )
        {
                rate = 0.;
        }
        else
        {
                rate = 1.6E-7 * pow(kT_eV, 0.5) * g_n / g_p / ( kT_eV + Gamma_np ) *
                        ( Anp * log(0.3*kT_eV/ryd + delta_np) + Bnp );
        }

        col_str = rate / COLL_CONST * phycon.sqrte * StatesElem[ipISO][nelem][ipHi].g;

        return( col_str );
}

---