cloudy trunk
Loading...
Searching...
No Matches
helike_recom.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*HeRecom - do recomb coef for He, called by HeLike */
4/*cross_section - calculates the photoionization cross_section for a given level and photon energy*/
5/*radrecomb - calculates radiative recombination coefficients. */
6/*He_cross_section returns cross section (cm^-2),
7 * given EgammaRyd, the photon energy in Ryd,
8 * ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
9 * nelem is charge, equal to 1 for Helium
10 * this is a wrapper for cross_section */
11/*Recomb_Seaton59 - find recombination for given n,using Seaton 59 approximation.
12 * The following three are needed by Recomb_Seaton59:
13 * ExponentialInt
14 * X1Int
15 * X2Int */
16
17#include "cddefines.h"
18#include "physconst.h"
19#include "hydro_bauman.h"
20#include "iso.h"
21#include "helike.h"
22#include "helike_recom.h"
23#include "thirdparty.h"
24#include "dense.h"
25#include "opacity.h"
26#include "atmdat.h"
27#include "taulines.h"
28
29/* The three of these are used in the computation of Recomb_Seaton59 */
30STATIC double ExponentialInt( double v );
31STATIC double X1Int( double u );
32STATIC double X2Int( double u );
33
34STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s);
35STATIC double GetHS98CrossSection( long n, long l, long s, double EgammaRyd );
36
37static double Xn_S59;
38
39/*He_cross_section returns cross section (cm^-2),
40 * given EgammaRyd, the photon energy in Ryd,
41 * ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
42 * nelem is charge, equal to 1 for Helium
43 * this is a wrapper for cross_section */
44double He_cross_section( double EgammaRyd , double EthRyd, long n, long l, long S, long nelem )
45{
46 // get cross section in megabarns
47 double cs = cross_section( EgammaRyd, EthRyd, nelem, n, l, S );
48
49 // rescale low-lying He values to Hummer & Storey 98, Table 1 Extrapolated
50 if( nelem==ipHELIUM && n <=5 && l<=2 )
51 {
52 double rescaled[31] = {
53 7.394,
54 5.485, 9.219, 15.985, 15.985, 15.985, 13.504,
55 8.018, 14.417, 28.501, 18.486, 18.132, 27.009,
56 10.721, 20.235, 41.568, 36.717, 35.766, -1.000, -1.000, 41.787, // };
57 13.527, 26.539, 55.692, 55.010, 53.514, -1.000, -1.000, -1.000, -1.000, 58.120 };
58 long ipLev = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[n][l][S];
59 ASSERT( rescaled[ipLev] > 0. );
60 cs *= rescaled[ipLev]/cross_section( EthRyd, EthRyd, nelem, n, l, S );
61 }
62
63 // convert to cm^-2
64 return cs * (1.e-18);
65}
66
67/*cross_section calculates the photoionization cross_section for a given level and photon energy
68 * this routine returns megabarns */
69STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long S)
70{
71 /* These fit parameters (E0, sigma, y_a, P, y_w, yzero, and yone) all come from the following work: */
72 /* >>refer He pcs Verner, D. A., Verner, E. M., \& Ferland , G. J. 1996,
73 * >>refercon Atomic Data and Nuclear Data Tables, Vol. 64, p.1 */
74 double E0[29] = {
75 1.36E+01,2.01E+01,1.76E+01,3.34E+01,4.62E+01,6.94E+01,8.71E+01,1.13E+02,1.59E+02,2.27E+02,
76 2.04E+02,2.74E+02,2.75E+02,3.38E+02,4.39E+02,4.17E+02,4.47E+02,5.18E+02,6.30E+02,6.27E+02,
77 8.66E+02,7.67E+02,9.70E+02,9.66E+02,1.06E+03,1.25E+03,1.35E+03,1.43E+03,1.56E+03};
78 double sigma[29] = {
79 9.49E+02,3.20E+02,5.46E+02,2.85E+02,2.34E+02,1.52E+02,1.33E+02,1.04E+02,6.70E+01,4.00E+01,
80 6.14E+01,4.04E+01,4.75E+01,3.65E+01,2.45E+01,3.14E+01,3.11E+01,2.59E+01,1.94E+01,2.18E+01,
81 1.23E+01,1.76E+01,1.19E+01,1.31E+01,1.20E+01,9.05E+00,8.38E+00,8.06E+00,7.17E+00};
82 double y_a[29] = {
83 1.47E+00,7.39E+00,1.72E+01,2.16E+01,2.18E+01,2.63E+01,2.54E+01,2.66E+01,3.35E+01,5.32E+01,
84 2.78E+01,3.57E+01,2.85E+01,3.25E+01,4.41E+01,3.16E+01,3.04E+01,3.28E+01,3.92E+01,3.45E+01,
85 5.89E+01,3.88E+01,5.35E+01,4.83E+01,5.77E+01,6.79E+01,7.43E+01,7.91E+01,9.10E+01};
86 double P[29] = {
87 3.19E+00,2.92E+00,3.16E+00,2.62E+00,2.58E+00,2.32E+00,2.34E+00,2.26E+00,2.00E+00,1.68E+00,
88 2.16E+00,1.92E+00,2.14E+00,2.00E+00,1.77E+00,2.04E+00,2.09E+00,2.02E+00,1.86E+00,2.00E+00,
89 1.62E+00,1.93E+00,1.70E+00,1.79E+00,1.72E+00,1.61E+00,1.59E+00,1.58E+00,1.54E+00};
90 double y_w[29] =
91 {2.039,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
92 double yzero[29] =
93 {0.4434,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
94 double yone[29] =
95 {2.136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
96
97 double pcs,Egamma,y,F,x;
98 double rel_photon_energy;
99
100 Egamma = EgammaRyd * EVRYD;
101
102 /* >>chng 02 apr 24, more protection against calling with too small an energy */
103 /* evaluating H-like photo cs at He energies, may be below threshold -
104 * prevent this from happening */
105 rel_photon_energy = EgammaRyd / EthRyd;
106 rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );
107
108 long s=0;//portland group 11.2 trips without =0, does not recognized that TotalInsanity does not return
109 if( S == 1 )
110 s=0;
111 else if( S == 3 )
112 s=1;
113 else
115
116 if( nelem==ipHELIUM && n<=25 && l<=4 )
117 {
118 // use Hummer & Storey 1998 cross-sections
119 pcs = GetHS98CrossSection( n, l, s, EgammaRyd );
120 }
121 else if( nelem==ipHELIUM && n>25 && l<=2 )
122 {
123 double scale[3][2] = {
124 {1.4673,1.3671},
125 {1.5458,1.5011},
126 {1.4912,1.5144}};
127
128 long ipLev = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[25][l][S];
129 double EthRyd_25 = iso_sp[ipHE_LIKE][nelem].fb[ipLev].xIsoLevNIonRyd;
130 pcs = GetHS98CrossSection( 25, l, s, EthRyd_25 * rel_photon_energy );
131 pcs *= pow((double)n/25., scale[l][s]);
132 }
133 else if( n==1 )
134 {
135 /* >>refer Helike PCS Verner, D.A., Ferland, G.J., Korista, K.T., & Yakovlev, D.G.
136 * >>refercon 1996a, ApJ 465,487 */
137 /* All helike (but not Helium itself) ground state cross-sections calculated here. */
138 x = Egamma/E0[nelem-1] - yzero[nelem-1];
139 y = sqrt(x*x + yone[nelem-1]*yone[nelem-1]);
140 F = ((x-1)*(x-1)+y_w[nelem-1]*y_w[nelem-1])
141 * pow(y,0.5*P[nelem-1]-5.5) * pow((1+sqrt(y/y_a[nelem-1])),-P[nelem-1]);
142 pcs = sigma[nelem-1]*F;
143 }
144 else if( nelem>=ipLITHIUM && nelem<=ipCALCIUM && n<11 && OP_Helike_NumPts[nelem][n][l][s]>0 )
145 {
146 // use TOPbase cross-sections
147 long numDataPoints = OP_Helike_NumPts[nelem][n][l][s];
148 ASSERT( numDataPoints > 0 );
149 ASSERT( OP_Helike_Xsectn[nelem][n][l][s] != NULL );
150
151 if( EgammaRyd < OP_Helike_Energy[nelem][n][l][s][numDataPoints-1] )
152 {
153 pcs = linint( OP_Helike_Energy[nelem][n][l][s], OP_Helike_Xsectn[nelem][n][l][s], numDataPoints, EgammaRyd );
154 }
155 else
156 {
157 // use a E^-3 tail
158 pcs = OP_Helike_Xsectn[nelem][n][l][s][numDataPoints-1] * POW3( OP_Helike_Energy[nelem][n][l][s][numDataPoints-1]/EgammaRyd );
159 }
160 }
161 else
162 {
163 /* To everything else we apply a hydrogenic routine. */
164 pcs = (1.e18)*H_photo_cs(rel_photon_energy , n, l, nelem);
165 }
166
167 ASSERT( pcs > 0. && pcs < 1.E10 );
168
169 return pcs;
170}
171
172STATIC double GetHS98CrossSection( long n, long l, long s, double EgammaRyd )
173{
174 double pcs;
175 ASSERT( n<=25 );
176 ASSERT( l<=4 );
177 ASSERT( s==0 || s==1 );
178
179 // use Hummer & Storey 1998 cross-sections
180 if( EgammaRyd < HS_He1_Energy[n][l][s][NUM_HS98_DATA_POINTS-1] )
181 {
182 pcs = linint( HS_He1_Energy[n][l][s], HS_He1_Xsectn[n][l][s], NUM_HS98_DATA_POINTS, EgammaRyd );
183 }
184 else
185 {
186 // use a E^-3 tail
187 pcs = HS_He1_Xsectn[n][l][s][NUM_HS98_DATA_POINTS-1] * POW3( HS_He1_Energy[n][l][s][NUM_HS98_DATA_POINTS-1]/EgammaRyd );
188 }
189
190 return pcs;
191}
192
193/* >>refer He-like RR Seaton, M.J. 1959, MNRAS 119, 81S */
194double Recomb_Seaton59( long nelem, double temp, long n)
195{
196 double lambda = TE1RYD * nelem * nelem / temp;
197 /* smallest x ever used here should be lowest Z, highest T, highest n...
198 * using helium, logt = 10., and n = 1000, we get xmin = 1.5789E-11. */
199 double x = lambda / n / n;
200 double AlphaN;
201 double SzeroOfX = 0.;
202 double SoneOfX = 0.;
203 double StwoOfX = 0.;
204 double SnOfLambda = 0.;
205 double lowerlimit, upperlimit, step;
206
207 Xn_S59 = x;
208
209 /* Equation 12 */
210 lowerlimit = x;
211 step = 3. * x;
212 upperlimit = lowerlimit + step;
213 SzeroOfX = qg32( lowerlimit, upperlimit, ExponentialInt);
214
215 do
216 {
217 lowerlimit = upperlimit;
218 step *= 2;
219 upperlimit = lowerlimit + step;
220 SzeroOfX += qg32( lowerlimit, upperlimit, ExponentialInt);
221 } while ( upperlimit < 20. );
222
223 /* This must be placed inside integral...too big to be
224 * handled separately.
225 SzeroOfX *= exp( x ); */
226
227 /* Equations 13 and 14 */
228 lowerlimit = 0.;
229 step = 0.5;
230 upperlimit = lowerlimit + step;
231 SoneOfX = qg32( lowerlimit, upperlimit, X1Int);
232 StwoOfX = qg32( lowerlimit, upperlimit, X2Int);
233
234 do
235 {
236 lowerlimit = upperlimit;
237 step *= 2;
238 upperlimit = lowerlimit + step;
239 SoneOfX += qg32( lowerlimit, upperlimit, X1Int);
240 StwoOfX += qg32( lowerlimit, upperlimit, X2Int);
241 } while ( upperlimit < 200. );
242
243 SoneOfX *= 0.1728 * pow( x, 1./3. );
244 StwoOfX *= -0.0496 * pow( x, 2./3. );
245
246 /* Equation 11 */
247 SnOfLambda = SzeroOfX + pow(1./lambda, 1./3.)*SoneOfX + pow(1./lambda, 2./3.)*StwoOfX;
248
249 AlphaN = 5.197E-14 * nelem * pow(x, 1.5) * SnOfLambda;
250
251 return AlphaN;
252
253}
254
255STATIC double ExponentialInt( double v )
256{
257 double Integrand;
258
259 Integrand = exp( -1. * v + Xn_S59) / v;
260
261 return Integrand;
262}
263
264STATIC double X1Int( double u )
265{
266 double Integrand;
267
268 Integrand = pow(1./(u + 1.), 5./3.) * (u - 1.) * exp( -1. * Xn_S59 * u );
269
270 return Integrand;
271}
272
273STATIC double X2Int( double u )
274{
275 double Integrand;
276
277 Integrand = pow(1./(u + 1.), 7./3.) * (u*u + 4./3.*u + 1.) * exp( -1. * Xn_S59 * u );
278
279 return Integrand;
280}
281
double ***** OP_Helike_Xsectn
Definition atmdat.cpp:10
double **** HS_He1_Xsectn
Definition atmdat.cpp:8
double **** HS_He1_Energy
Definition atmdat.cpp:9
long **** OP_Helike_NumPts
Definition atmdat.cpp:12
double ***** OP_Helike_Energy
Definition atmdat.cpp:11
#define NUM_HS98_DATA_POINTS
Definition atmdat.h:127
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define POW3
Definition cddefines.h:936
double qg32(double, double, double(*)(double))
Definition service.cpp:1053
const int ipLITHIUM
Definition cddefines.h:307
const int ipHELIUM
Definition cddefines.h:306
#define MAX2
Definition cddefines.h:782
const int ipCALCIUM
Definition cddefines.h:324
NORETURN void TotalInsanity(void)
Definition service.cpp:886
STATIC double GetHS98CrossSection(long n, long l, long s, double EgammaRyd)
static double Xn_S59
STATIC double X2Int(double u)
double He_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long nelem)
STATIC double ExponentialInt(double v)
STATIC double X1Int(double u)
double Recomb_Seaton59(long nelem, double temp, long n)
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipHE_LIKE
Definition iso.h:63
#define S(I_, J_)
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double TE1RYD
Definition physconst.h:183
double linint(const double x[], const double y[], long n, double xval)