cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*hlike_radrecomb_from_cross_section - integrate Milne relation */ 00004 /*RecomInt - Integral in Milne relation. Called by qg32. */ 00005 /*H_cross_section returns cross section (cm^-2), 00006 * given EgammaRyd, the photon energy in Ryd, 00007 * ipLevel, the index of the level, 0 is ground, 00008 * nelem is charge, equal to 0 for Hydrogen */ 00009 #include "cddefines.h" 00010 #include "physconst.h" 00011 #include "hydro_bauman.h" 00012 #include "iso.h" 00013 #include "helike.h" 00014 #include "dense.h" 00015 #include "opacity.h" 00016 #include "atmdat.h" 00017 00018 static double RecomInt(double EE); 00019 00020 static double kTRyd,EthRyd; 00021 static long int ipLev,globalZ; 00022 00023 /*H_cross_section returns cross section (cm^-2), 00024 * given EgammaRyd, the photon energy in Ryd, 00025 * ipLevel, the index of the level, 0 is ground, 00026 * nelem is charge, equal to 0 for Hydrogen */ 00027 double H_cross_section( double EgammaRyd , long ipLev , long nelem ) 00028 { 00029 double cs; 00030 double rel_photon_energy; 00031 long ipISO = ipH_LIKE; 00032 00033 /* make sure this routine not called within collapsed high levels */ 00034 ASSERT( ipLev < iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem] ); 00035 ASSERT( ipLev > ipH1s); 00036 00037 EthRyd = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipLev]; 00038 00039 /* >>chng 02 apr 24, more protection against calling with too small an energy */ 00040 /* evaluating H-like photo cs at He energies, may be below threshold - 00041 * prevent this from happening */ 00042 rel_photon_energy = EgammaRyd / EthRyd; 00043 rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. ); 00044 00045 cs = H_photo_cs(rel_photon_energy , N_(ipLev), L_(ipLev), nelem + 1 ); 00046 00047 ASSERT( cs > 0. && cs < 1.E-8 ); 00048 00049 return cs; 00050 } 00051 00052 /*hlike_radrecomb_from_cross_section calculates radiative recombination coefficients. */ 00053 double hlike_radrecomb_from_cross_section(double temp, long nelem, long ipLo) 00054 { 00055 double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive; 00056 double change[5] = {0.,0.,0.,0.,0.}; 00057 long ipISO = ipH_LIKE; 00058 00059 if( ipLo == 0 ) 00060 return t_ADfA::Inst().H_rad_rec(nelem+1,ipLo,temp); 00061 00062 EthRyd = iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]; 00063 00064 /* Factors outside integral in Milne relation. */ 00065 b = MILNE_CONST * StatesElem[ipISO][nelem][ipLo].g * pow(temp,-1.5) / 2.; 00066 /* kT in Rydbergs. */ 00067 kTRyd = temp / TE1RYD; 00068 globalZ = nelem; 00069 ipLev = ipLo; 00070 00071 /* Begin integration. */ 00072 /* First define characteristic step */ 00073 E1 = EthRyd; 00074 step = MIN2( 0.125*kTRyd, 0.5*E1 ); 00075 E2 = E1 + step; 00076 /* Perform initial integration, from threshold to threshold + step. */ 00077 RecomIntegral = qg32( E1, E2, RecomInt); 00078 /* Repeat the integration, adding each new result to the total, 00079 * except that the step size is doubled every time, since values away from 00080 * threshold tend to fall off more slowly. */ 00081 do 00082 { 00083 OldRecomIntegral = RecomIntegral; 00084 E1 = E2; 00085 step *= 1.25; 00086 E2 = E1 + step; 00087 RecomIntegral += qg32( E1, E2, RecomInt); 00088 change[4] = change[3]; 00089 change[3] = change[2]; 00090 change[2] = change[1]; 00091 change[1] = change[0]; 00092 change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral; 00093 TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4]; 00094 /* Continue integration until the upper limit exceeds 100kTRyd, an arbitrary 00095 * point at which the integrand component exp(electron energy/kT) is very small, 00096 * making neglible cross-sections at photon energies beyond that point, 00097 * OR when the last five steps resulted in less than a 1 percent change. */ 00098 } while ( ((E2-EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) ); 00099 00100 /* Calculate recombination coefficient. */ 00101 alpha = b * RecomIntegral; 00102 00103 alpha = MAX2( alpha, SMALLDOUBLE ); 00104 00105 return alpha; 00106 } 00107 00108 /*RecomInt, used in comput milne relation for he rec - the energy is photon Rydbergs. */ 00109 static double RecomInt(double ERyd) 00110 { 00111 double x1, temp; 00112 00113 /* Milne relation integrand */ 00114 x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - EthRyd ) / kTRyd); 00115 temp = H_cross_section( ERyd , ipLev, globalZ ); 00116 x1 *= temp; 00117 00118 return x1; 00119 }