cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hydro_recom.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*hlike_radrecomb_from_cross_section - integrate Milne relation */
4 /*RecomInt - Integral in Milne relation. Called by qg32. */
5 /*H_cross_section returns cross section (cm^-2),
6  * given EgammaRyd, the photon energy in Ryd,
7  * ipLevel, the index of the level, 0 is ground,
8  * nelem is charge, equal to 0 for Hydrogen */
9 #include "cddefines.h"
10 #include "physconst.h"
11 #include "hydro_bauman.h"
12 #include "iso.h"
13 #include "helike.h"
14 #include "dense.h"
15 #include "opacity.h"
16 #include "atmdat.h"
17 
18 static double RecomInt(double EE);
19 
20 static double kTRyd,EthRyd;
21 static long int ipLev,globalZ;
22 
23 /*H_cross_section returns cross section (cm^-2),
24  * given EgammaRyd, the photon energy in Ryd,
25  * ipLevel, the index of the level, 0 is ground,
26  * nelem is charge, equal to 0 for Hydrogen */
27 double H_cross_section( double EgammaRyd , long ipLev , long nelem )
28 {
29  double cs;
30  double rel_photon_energy;
31  long ipISO = ipH_LIKE;
32 
33  /* make sure this routine not called within collapsed high levels */
34  ASSERT( ipLev < iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem] );
35  ASSERT( ipLev > ipH1s);
36 
38 
39  /* >>chng 02 apr 24, more protection against calling with too small an energy */
40  /* evaluating H-like photo cs at He energies, may be below threshold -
41  * prevent this from happening */
42  rel_photon_energy = EgammaRyd / EthRyd;
43  rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );
44 
45  cs = H_photo_cs(rel_photon_energy , N_(ipLev), L_(ipLev), nelem + 1 );
46 
47  ASSERT( cs > 0. && cs < 1.E-8 );
48 
49  return cs;
50 }
51 
52 /*hlike_radrecomb_from_cross_section calculates radiative recombination coefficients. */
53 double hlike_radrecomb_from_cross_section(double temp, long nelem, long ipLo)
54 {
55  double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive;
56  double change[5] = {0.,0.,0.,0.,0.};
57  long ipISO = ipH_LIKE;
58 
59  if( ipLo == 0 )
60  return t_ADfA::Inst().H_rad_rec(nelem+1,ipLo,temp);
61 
62  EthRyd = iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];
63 
64  /* Factors outside integral in Milne relation. */
65  b = MILNE_CONST * StatesElem[ipISO][nelem][ipLo].g * pow(temp,-1.5) / 2.;
66  /* kT in Rydbergs. */
67  kTRyd = temp / TE1RYD;
68  globalZ = nelem;
69  ipLev = ipLo;
70 
71  /* Begin integration. */
72  /* First define characteristic step */
73  E1 = EthRyd;
74  step = MIN2( 0.125*kTRyd, 0.5*E1 );
75  E2 = E1 + step;
76  /* Perform initial integration, from threshold to threshold + step. */
77  RecomIntegral = qg32( E1, E2, RecomInt);
78  /* Repeat the integration, adding each new result to the total,
79  * except that the step size is doubled every time, since values away from
80  * threshold tend to fall off more slowly. */
81  do
82  {
83  OldRecomIntegral = RecomIntegral;
84  E1 = E2;
85  step *= 1.25;
86  E2 = E1 + step;
87  RecomIntegral += qg32( E1, E2, RecomInt);
88  change[4] = change[3];
89  change[3] = change[2];
90  change[2] = change[1];
91  change[1] = change[0];
92  change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral;
93  TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4];
94  /* Continue integration until the upper limit exceeds 100kTRyd, an arbitrary
95  * point at which the integrand component exp(electron energy/kT) is very small,
96  * making neglible cross-sections at photon energies beyond that point,
97  * OR when the last five steps resulted in less than a 1 percent change. */
98  } while ( ((E2-EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) );
99 
100  /* Calculate recombination coefficient. */
101  alpha = b * RecomIntegral;
102 
103  alpha = MAX2( alpha, SMALLDOUBLE );
104 
105  return alpha;
106 }
107 
108 /*RecomInt, used in comput milne relation for he rec - the energy is photon Rydbergs. */
109 static double RecomInt(double ERyd)
110 {
111  double x1, temp;
112 
113  /* Milne relation integrand */
114  x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - EthRyd ) / kTRyd);
115  temp = H_cross_section( ERyd , ipLev, globalZ );
116  x1 *= temp;
117 
118  return x1;
119 }

Generated for cloudy by doxygen 1.8.1.2