cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
helike_recom.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 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 "helike_recom.h"
19 #include "hydro_bauman.h"
20 #include "iso.h"
21 #include "thirdparty.h"
22 #include "atmdat.h"
23 #include "freebound.h"
24 #include "integrate.h"
25 
26 /* The three of these are used in the computation of Recomb_Seaton59 */
27 STATIC double ExponentialInt( double v );
28 STATIC double X1Int( double u );
29 STATIC double X2Int( double u );
30 
31 STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s);
32 STATIC double GetHS98CrossSection( long n, long l, long s, double EgammaRyd );
33 
34 static double Xn_S59;
35 
36 /*He_cross_section returns cross section (cm^-2),
37  * given EgammaRyd, the photon energy in Ryd,
38  * ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
39  * nelem is charge, equal to 1 for Helium
40  * this is a wrapper for cross_section */
41 double He_cross_section( double EgammaRyd , double EthRyd, long n, long l, long S, long nelem )
42 {
43  // get cross section in megabarns
44  double cs = cross_section( EgammaRyd, EthRyd, nelem, n, l, S );
45 
46  // rescale low-lying He values to Hummer & Storey 98, Table 1 Extrapolated
47  if( nelem==ipHELIUM && n <=5 && l<=2 )
48  {
49  static const double rescaled[31] = {
50  7.394,
51  5.485, 9.219, 15.985, 15.985, 15.985, 13.504,
52  8.018, 14.417, 28.501, 18.486, 18.132, 27.009,
53  10.721, 20.235, 41.568, 36.717, 35.766, -1.000, -1.000, 41.787, // };
54  13.527, 26.539, 55.692, 55.010, 53.514, -1.000, -1.000, -1.000, -1.000, 58.120 };
55  long ipLev = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[n][l][S];
56  ASSERT( rescaled[ipLev] > 0. );
57  cs *= rescaled[ipLev]/cross_section( EthRyd, EthRyd, nelem, n, l, S );
58  }
59 
60  // convert to cm^-2
61  return cs * (1.e-18);
62 }
63 
64 /*cross_section calculates the photoionization cross_section for a given level and photon energy
65  * this routine returns megabarns */
66 STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long S)
67 {
68  /* These fit parameters (E0, sigma, y_a, P, y_w, yzero, and yone) all come from the following work: */
69  /* >>refer He pcs Verner, D. A., Verner, E. M., \& Ferland , G. J. 1996,
70  * >>refercon Atomic Data and Nuclear Data Tables, Vol. 64, p.1 */
71  static const double E0[29] = {
72  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,
73  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,
74  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};
75  static const double sigma[29] = {
76  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,
77  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,
78  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};
79  static const double y_a[29] = {
80  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,
81  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,
82  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};
83  static const double P[29] = {
84  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,
85  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,
86  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};
87  static const double y_w[29] =
88  {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};
89  static const double yzero[29] =
90  {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};
91  static const double yone[29] =
92  {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};
93 
94  double pcs,Egamma,y,F,x;
95  double rel_photon_energy;
96 
97  Egamma = EgammaRyd * EVRYD;
98 
99  /* >>chng 02 apr 24, more protection against calling with too small an energy */
100  /* evaluating H-like photo cs at He energies, may be below threshold -
101  * prevent this from happening */
102  rel_photon_energy = EgammaRyd / EthRyd;
103  rel_photon_energy = MAX2( rel_photon_energy , 1. + FLT_EPSILON*2. );
104 
105  long s=0;//portland group 11.2 trips without =0, does not recognized that TotalInsanity does not return
106  if( S == 1 )
107  s=0;
108  else if( S == 3 )
109  s=1;
110  else
111  TotalInsanity();
112 
113  if( nelem==ipHELIUM && n<=25 && l<=4 )
114  {
115  // use Hummer & Storey 1998 cross-sections
116  pcs = GetHS98CrossSection( n, l, s, EgammaRyd );
117  }
118  else if( nelem==ipHELIUM && n>25 && l<=2 )
119  {
120  static const double scale[3][2] = {
121  {1.4673,1.3671},
122  {1.5458,1.5011},
123  {1.4912,1.5144}};
124 
125  long ipLev = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[25][l][S];
126  double EthRyd_25 = iso_sp[ipHE_LIKE][nelem].fb[ipLev].xIsoLevNIonRyd;
127  pcs = GetHS98CrossSection( 25, l, s, EthRyd_25 * rel_photon_energy );
128  pcs *= pow((double)n/25., scale[l][s]);
129  }
130  else if( n==1 )
131  {
132  /* >>refer Helike PCS Verner, D.A., Ferland, G.J., Korista, K.T., & Yakovlev, D.G.
133  * >>refercon 1996a, ApJ 465,487 */
134  /* All helike (but not Helium itself) ground state cross-sections calculated here. */
135  x = Egamma/E0[nelem-1] - yzero[nelem-1];
136  y = sqrt(x*x + yone[nelem-1]*yone[nelem-1]);
137  F = ((x-1)*(x-1)+y_w[nelem-1]*y_w[nelem-1])
138  * pow(y,0.5*P[nelem-1]-5.5) * pow((1+sqrt(y/y_a[nelem-1])),-P[nelem-1]);
139  pcs = sigma[nelem-1]*F;
140  }
141  else if( nelem>=ipLITHIUM && nelem<=ipCALCIUM && n<11 && OP_Helike_NumPts[nelem][n][l][s]>0 )
142  {
143  // use TOPbase cross-sections
144  long numDataPoints = OP_Helike_NumPts[nelem][n][l][s];
145  ASSERT( numDataPoints > 0 );
146  ASSERT( OP_Helike_Xsectn[nelem][n][l][s] != NULL );
147 
148  if( EgammaRyd < OP_Helike_Energy[nelem][n][l][s][numDataPoints-1] )
149  {
150  pcs = linint( OP_Helike_Energy[nelem][n][l][s], OP_Helike_Xsectn[nelem][n][l][s], numDataPoints, EgammaRyd );
151  }
152  else
153  {
154  // use a E^-3 tail
155  pcs = OP_Helike_Xsectn[nelem][n][l][s][numDataPoints-1] * POW3( OP_Helike_Energy[nelem][n][l][s][numDataPoints-1]/EgammaRyd );
156  }
157  }
158  else
159  {
160  /* To everything else we apply a hydrogenic routine. */
161  pcs = (1.e18)*H_photo_cs(rel_photon_energy , n, l, nelem);
162  }
163 
164  ASSERT( pcs > 0. && pcs < 1.E10 );
165 
166  return pcs;
167 }
168 
169 STATIC double GetHS98CrossSection( long n, long l, long s, double EgammaRyd )
170 {
171  double pcs;
172  ASSERT( n<=25 );
173  ASSERT( l<=4 );
174  ASSERT( s==0 || s==1 );
175 
176  // use Hummer & Storey 1998 cross-sections
177  if( EgammaRyd < HS_He1_Energy[n][l][s][NUM_HS98_DATA_POINTS-1] )
178  {
179  pcs = linint( HS_He1_Energy[n][l][s], HS_He1_Xsectn[n][l][s], NUM_HS98_DATA_POINTS, EgammaRyd );
180  }
181  else
182  {
183  // use a E^-3 tail
184  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 );
185  }
186 
187  return pcs;
188 }
189 
190 #if 1
191 /* >>refer He-like RR Seaton, M.J. 1959, MNRAS 119, 81S */
192 double Recomb_Seaton59( long nelem, double temp, long n)
193 {
194  double lambda = TE1RYD * nelem * nelem / temp;
195  /* smallest x ever used here should be lowest Z, highest T, highest n...
196  * using helium, logt = 10., and n = 1000, we get xmin = 1.5789E-11. */
197  double x = lambda / n / n;
198  double AlphaN;
199  double SzeroOfX = 0.;
200  double SoneOfX = 0.;
201  double StwoOfX = 0.;
202  double SnOfLambda = 0.;
203  double lowerlimit, upperlimit, step;
204 
205  fixit("the variant below should be faster and more accurate, but needs more testing");
206  Xn_S59 = x;
207 
208  /* Equation 12 */
209  lowerlimit = x;
210  step = 3. * x;
211  upperlimit = lowerlimit + step;
212  SzeroOfX = qg32( lowerlimit, upperlimit, ExponentialInt);
213 
214  do
215  {
216  lowerlimit = upperlimit;
217  step *= 2;
218  upperlimit = lowerlimit + step;
219  SzeroOfX += qg32( lowerlimit, upperlimit, ExponentialInt);
220  } while ( upperlimit < 20. );
221 
222  /* This must be placed inside integral...too big to be
223  * handled separately.
224  SzeroOfX *= exp( x ); */
225 
226  /* Equations 13 and 14 */
227  lowerlimit = 0.;
228  step = 0.5;
229  upperlimit = lowerlimit + step;
230  SoneOfX = qg32(lowerlimit, upperlimit, X1Int);
231  StwoOfX = qg32(lowerlimit, upperlimit, X2Int);
232 
233  do
234  {
235  lowerlimit = upperlimit;
236  step *= 2;
237  upperlimit = lowerlimit + step;
238  SoneOfX += qg32( lowerlimit, upperlimit, X1Int);
239  StwoOfX += qg32( lowerlimit, upperlimit, X2Int);
240  } while ( upperlimit < 200. );
241 
242  SoneOfX *= 0.1728 * powpq( x, 1, 3 );
243  StwoOfX *= -0.0496 * powpq( x, 2, 3 );
244 
245  /* Equation 11 */
246  SnOfLambda = SzeroOfX + powpq(1./lambda, 1, 3)*SoneOfX + powpq(1./lambda, 2, 3)*StwoOfX;
247 
248  if( false )
249  {
250  double SSzeroOfX = e1_scaled(x);
251  double gamma13 = tgamma(1./3.);
252  double gamma23 = PI2/(sqrt(3.)*gamma13);
253  double x13 = cbrt(x);
254  double x23 = pow2(x13);
255  double SSoneOfX = 0.1728*((3.*x+1)*igamc_scaled(1./3.,x)*gamma13 - 3.*x13);
256  double SStwoOfX = -0.0496*(((1.5*x+2.)*x+1.)*igamc_scaled(2./3.,x)*gamma23 - 1.5*x23*(x+1.));
257  double SSnOfLambda = SSzeroOfX + powpq(1./lambda, 1, 3)*SSoneOfX + powpq(1./lambda, 2, 3)*SStwoOfX;
258 
259  dprintf( ioQQQ, "%e %e %e %e %e old %e new %e\n", x, SzeroOfX/SSzeroOfX-1., SoneOfX/SSoneOfX-1.,
260  StwoOfX/SStwoOfX-1., SnOfLambda/SSnOfLambda-1., SnOfLambda, SSnOfLambda );
261  }
262 
263  AlphaN = 5.197E-14 * nelem * powpq(x, 3, 2) * SnOfLambda;
264 
265  return AlphaN;
266 }
267 
268 STATIC double ExponentialInt( double v )
269 {
270  double Integrand;
271 
272  Integrand = exp( -v + Xn_S59) / v;
273 
274  return Integrand;
275 }
276 
277 STATIC double X1Int( double u )
278 {
279  double Integrand;
280 
281  Integrand = powpq(1./(u + 1.), 5, 3) * (u - 1.) * exp( -Xn_S59 * u );
282 
283  return Integrand;
284 }
285 
286 STATIC double X2Int( double u )
287 {
288  double Integrand;
289 
290  Integrand = powpq(1./(u + 1.), 7, 3) * (u*u + 4./3.*u + 1.) * exp( -Xn_S59 * u );
291 
292  return Integrand;
293 }
294 #else
295 /* >>refer He-like RR Seaton, M.J. 1959, MNRAS 119, 81S */
296 double Recomb_Seaton59(long nelem, double temp, long n)
297 {
298  double lambda = TE1RYD * nelem * nelem / temp;
299  /* smallest x ever used here should be lowest Z, highest T, highest n...
300  * using helium, logt = 10., and n = 1000, we get xmin = 1.5789E-11. */
301  double x = lambda / n / n;
302 
303  /* Equation 12 */
304  double SzeroOfX = e1_scaled(x);
305 
306  /* Equations 13 and 14 */
307  double gamma13 = tgamma(1./3.);
308  double gamma23 = PI2/(sqrt(3.)*gamma13);
309  double x13 = cbrt(x);
310  double x23 = pow2(x13);
311  double SoneOfX = 0.172826*((3.*x+1)*igamc_scaled(1./3.,x)*gamma13 - 3.*x13);
312  double StwoOfX = -0.0495957*(((1.5*x+2.)*x+1.)*igamc_scaled(2./3.,x)*gamma23 - 1.5*x23*(x+1.));
313 
314  /* Equation 11 */
315  double z13 = cbrt(1./lambda);
316  double SnOfLambda = (StwoOfX*z13 + SoneOfX)*z13 + SzeroOfX;
317 
318  return 5.197e-14 * nelem * powpq(x,3,2) * SnOfLambda;
319 }
320 #endif
const int ipHE_LIKE
Definition: iso.h:65
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
static double Xn_S59
STATIC double ExponentialInt(double v)
FILE * ioQQQ
Definition: cddefines.cpp:7
vector< freeBound > fb
Definition: iso.h:481
STATIC double X1Int(double u)
STATIC double X2Int(double u)
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1102
#define NUM_HS98_DATA_POINTS
Definition: atmdat.h:269
#define STATIC
Definition: cddefines.h:118
STATIC double GetHS98CrossSection(long n, long l, long s, double EgammaRyd)
double **** HS_He1_Xsectn
Definition: atmdat.cpp:60
double ***** OP_Helike_Xsectn
Definition: atmdat.cpp:62
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
#define ASSERT(exp)
Definition: cddefines.h:613
const int ipCALCIUM
Definition: cddefines.h:368
double qg32(double, double, double(*)(double))
Definition: service.cpp:1175
T pow2(T a)
Definition: cddefines.h:981
double **** HS_He1_Energy
Definition: atmdat.cpp:61
double powpq(double x, int p, int q)
Definition: service.cpp:630
const int ipHELIUM
Definition: cddefines.h:350
double ***** OP_Helike_Energy
Definition: atmdat.cpp:63
double linint(const double x[], const double y[], long n, double xval)
long **** OP_Helike_NumPts
Definition: atmdat.cpp:64
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
double e1_scaled(double x)
#define MAX2(a, b)
Definition: cddefines.h:824
double igamc_scaled(double a, double x)
#define S(I_, J_)
double Recomb_Seaton59(long nelem, double temp, long n)
#define fixit(a)
Definition: cddefines.h:417
double He_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long nelem)
#define POW3
Definition: cddefines.h:986
const int ipLITHIUM
Definition: cddefines.h:351