cloudy trunk
Loading...
Searching...
No Matches
hydro_vs_rates.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/* CS_VS80 compute thermally averaged collision strength for collisional deexcitation
4 * of hydrogenic atoms, from
5 * >>refer H1 collision Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
6 *hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients */
7/*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for ions
8 * with all n, and Z*/
9#include "cddefines.h"
10#include "dense.h"
11#include "phycon.h"
12#include "physconst.h"
13#include "iso.h"
14#include "hydro_vs_rates.h"
15#include "lines_service.h"
16#include "taulines.h"
17
18STATIC double hydro_vs_coll_str( double energy, long ipISO, long nelem, long ipHi, long ipLo, long Collider, double Aul );
19
20namespace {
21 class my_Integrand : public std::unary_function<double, double>
22 {
23 public:
24 long ipISO, nelem, ipHi, ipLo, Collider;
25 double Aul;
26 double temp;
27
28 double operator() (double EOverKT)
29 {
30 double col_str = hydro_vs_coll_str( EOverKT * EVRYD * temp/TE1RYD, ipISO, nelem, ipHi, ipLo, Collider, Aul );
31 return exp( -1.*EOverKT ) * col_str;
32 }
33 };
34}
35
36/* These are masses relative to the proton mass of the electron, proton, and alpha particle. */
37static const double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
38
39/*
40 Neither of these rates can be modified to account for impact by non-electrons because they
41 are fits to thermally-averaged rates for electron impact...it can not be unravelled to
42 expose a projectile energy that can then be scaled to account for other projectiles.
43 Instead, we have included the original cross section formula (eq 14) from
44 Vriens & Smeets (1980) below.
45*/
46
47/* VS80 stands for Vriens and Smeets 1980 */
48/* This routine calculates thermally-averaged collision strengths. */
49double CS_VS80( long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider )
50{
51 double coll_str;
52
53 if( Collider == ipELECTRON )
54 {
55 coll_str = hydro_vs_deexcit( ipISO, nelem, ipHi, ipLo, Aul);
56 }
57 else
58 {
59 /* evaluate collision strength, two options,
60 * do thermal average (very slow) if set with
61 * SET COLLISION STRENGTH AVERAGE command,
62 * default just do point at kT
63 * tests show no impact on test suite, much faster */
64 if( iso_ctrl.lgCollStrenThermAver )
65 {
66 my_Integrand func;
67
68 func.ipISO = ipISO;
69 func.nelem = nelem;
70 func.ipHi = ipHi;
71 func.ipLo = ipLo;
72 func.temp = temp;
73 func.Collider = Collider;
74 func.Aul = Aul;
75
77 /* do average over Maxwellian */
78 coll_str = VS80.sum( 0., 1., func );
79 coll_str += VS80.sum( 1., 10., func );
80 }
81 else
82 {
83 /* the argument to the function is equivalent to evaluating the function at
84 * hnu = kT */
85 coll_str = hydro_vs_coll_str( EVRYD*temp/TE1RYD, ipISO, nelem, ipHi, ipLo, Collider, Aul );
86 }
87 }
88
89 ASSERT( coll_str >= 0. );
90 return coll_str;
91}
92
93/*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom,
94 * from Vriens and Smeets */
95STATIC double hydro_vs_coll_str( double energy, long ipISO, long nelem, long ipHi, long ipLo, long Collider, double Aul )
96{
97 double Apn, bp, Bpn, delta;
98 double Epi, Epn;
99 double gamma, p, n;
100 double ryd, s;
101 double cross_section, coll_str, gLo, gHi, abs_osc_str, reduced_mass;
102
103 DEBUG_ENTRY( "hydro_vs_coll_str()" );
104
105 // number of colliders is 4 in def, should be macro
106 ASSERT( Collider >= 0.&& Collider <4 );
107 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
108 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
109
110 gLo = iso_sp[ipISO][nelem].st[ipLo].g();
111 gHi = iso_sp[ipISO][nelem].st[ipHi].g();
112
113 /* This comes from equations 14,15, and 16 of Vriens and Smeets. */
114 /* >>refer he-like cs Vriens, L. \& Smeets, A. H. M. Phys. Rev. A, 22, 940 */
115 /* Computes the Vriens and Smeets
116 * rate coeff. for collisional dexcitation between any two levels of H.
117 * valid for all nhigh, nlow
118 * at end converts this excitation rate to collision strength */
119
120 n = (double)iso_sp[ipISO][nelem].st[ipHi].n();
121 p = (double)iso_sp[ipISO][nelem].st[ipLo].n();
122
123 ryd = EVRYD;
124 s = fabs(n-p);
125 ASSERT( s > 0. );
126
127 Epn = EVRYD*(iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd - iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
128 Epi = EVRYD*iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd;
129
130 /* This is an absorption oscillator strength. */
131 abs_osc_str = GetGF( Aul, Epn*RYD_INF/EVRYD, gHi)/gLo;
132 Apn = 2.*ryd/Epn*abs_osc_str;
133
134 bp = 1.4*log(p)/p - .7/p - .51/p/p + 1.16/p/p/p - 0.55/p/p/p/p;
135
136 Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));
137
138 delta = exp(-Bpn/Apn) - 0.4*Epn/ryd;
139
140 gamma = ryd*(8. + 23.*POW2(s/p))/
141 ( 8. + 1.1*n*s + 0.8/pow2(s) + 0.4*sqrt(pow3(n))/sqrt(s)*fabs(s-1.0) );
142
143 /* Scale the energy to get an equivalent electron energy. */
144 energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
145
146 /* cross section in units of PI*a_o^2 */
147 if( energy/2./ryd+delta <= 0 /*|| energy < Epn*/ )
148 {
149 cross_section = 0.;
150 }
151 else
152 {
153 cross_section = 2.*ryd/(energy + gamma)*(Apn*log(energy/2./ryd+delta) + Bpn);
155 }
156
157 /* convert to collision strength */
158 coll_str = ConvCrossSect2CollStr( cross_section * PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM, gLo, energy/EVRYD, reduced_mass );
159
160 ASSERT( coll_str >= 0. );
161
162 return( coll_str );
163}
164
165/*hydro_vs_coll_recomb generate hydrogenic collisional recombination rate coefficients */
166double hydro_vs_coll_recomb( double ionization_energy_Ryd, double Te, double stat_level, double stat_ion )
167{
168 double coef,
169 denom,
170 epi,
171 t_eV;
172
173 DEBUG_ENTRY( "hydro_vs_coll_recomb()" );
174
175 /* This is equation 9 of
176 * >>refer H1 collision recomb Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 */
177
178 /* this is kT in eV */
179 t_eV = Te/EVDEGK;
180
181 /* this is the ionization energy relative to kT, dimensionless */
182 epi = ionization_energy_Ryd * EVRYD / t_eV;
183
184 /* this is the denominator of equation 8 of VS80. */
185 denom = pow(epi,2.33) + 4.38*pow(epi,1.72) + 1.32*epi;
186
187 /* this is equation 9 of VS80 */
188 coef = 3.17e-27 / pow3(t_eV) * stat_level / stat_ion / denom;
189
190 ASSERT( coef >= 0. );
191 return( coef );
192}
193
194
195/*hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients */
196double hydro_vs_ioniz( double ionization_energy_Ryd, double Te )
197{
198 double coef,
199 denom,
200 epi,
201 t_eV;
202
203 DEBUG_ENTRY( "hydro_vs_ioniz()" );
204
205 /* a function written to calculate the rate coefficients
206 * for hydrogen collisional ionizations from
207 * Jason Ferguson, summer 94
208 * valid for all n
209 * >>refer H1 collision Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
210 * */
211
212 /* this is kT in eV */
213 t_eV = Te/EVDEGK;
214
215 /* this is the ionization energy relative to kT, dimensionless */
216 epi = ionization_energy_Ryd * EVRYD / t_eV;
217
218 /* this is the denominator of equation 8 of VS80. */
219 denom = pow(epi,2.33) + 4.38*pow(epi,1.72) + 1.32*epi;
220
221 /* this is equation 8 of VS80 */
222 coef = 9.56e-6 / sqrt(pow3(t_eV)) * dsexp( epi ) / denom;
223
224 ASSERT( coef >= 0. );
225 return( coef );
226}
227
228/*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for all n, and Z*/
230 /* the isoelectronic sequence */
231 long int ipISO ,
232 /* element, >=1 since only used for ions
233 * nelem = 1 is helium the least possible charge */
234 long int nelem,
235 /* principal quantum number, > 1
236 * since only used for excited states */
237 long int n,
238 double ionization_energy_Ryd,
239 double Te )
240{
241 double H,
242 HydColIon_v,
243 Rnp,
244 charge,
245 chim,
246 eone,
247 etwo,
248 ethree,
249 g,
250 rate,
251 rate2,
252 boltz,
253 t1,
254 t2,
255 t3,
256 t4,
257 tev,
258 xn,
259 y;
260 static const double arrH[4] = {1.48,3.64,5.93,8.32};
261 static const double arrRnp[8] = {2.20,1.90,1.73,1.65,1.60,1.56,1.54,1.52};
262 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};
263
264 static double small = 0.;
265
266 DEBUG_ENTRY( "Hion_coll_ioniz_ratecoef()" );
267 /*calculate hydrogenic ionization rates for all n, and Z
268 * >>refer HI cs Allen 1973, Astro. Quan. for low Te.
269 * >>refer HI cs Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
270 * */
271
272 charge = nelem - ipISO;
273 /* this routine only for ions, nelem=0 is H, nelem=1 he, etc */
274 ASSERT( charge > 0);
275 ASSERT( n>1 );
276
277 if( n > 4 )
278 {
279 H = 2.15*n;
280 }
281 else
282 {
283 H = arrH[n-1];
284 }
285
286 if( n > 8 )
287 {
288 Rnp = 1.52;
289 }
290 else
291 {
292 Rnp = arrRnp[n-1];
293 }
294
295 if( n > 10 )
296 {
297 g = arrg[9];
298 }
299 else
300 {
301 g = arrg[n-1];
302 }
303
304 tev = EVRYD/TE1RYD*Te;
305 xn = (double)n;
306 chim = EVRYD * ionization_energy_Ryd;
307 y = chim/tev;
308 boltz = dsexp( chim/tev );
309
310 eone = ee1(y);
311 etwo = boltz - y*eone;
312 ethree = (boltz - y*etwo)/2.;
313
314 t1 = 1/xn*eone;
315 t2 = 1./3./xn*(boltz - y*ethree);
316 t3 = 3.*H/xn/(3. - Rnp)*(y*etwo - 2.*y*eone + boltz);
317 t4 = 3.36*y*(eone - etwo);
318 rate = 7.69415e-9*sqrt(Te)*9.28278e-3*powi(xn/(charge+1),4)*g*y;
319 rate *= t1 - t2 + t3 + t4;
320 rate2 = 2.1e-8*sqrt(Te)/chim/chim*dsexp(2.302585*5040.*chim/Te);
321
322 /* don't let the rates go negative */
323 rate = MAX2(rate,small);
324 rate2 = MAX2(rate2,small);
325
326 /* Take the lowest of the two, they fit nicely together... */
327 /* >>chng 10 Sept 02, sometimes one of these is zero and the other is positive.
328 * in that case take the bigger one. */
329 if( rate==0. || rate2==0. )
330 HydColIon_v = MAX2(rate,rate2);
331 else
332 HydColIon_v = MIN2(rate,rate2);
333
334 ASSERT( HydColIon_v >= 0. );
335 return( HydColIon_v );
336}
337
338/*hydro_vs_deexcit compute collisional deexcitation rates for hydrogen atom,
339 * from Vriens and Smeets 1980 */
340double hydro_vs_deexcit( long ipISO, long nelem, long ipHi, long ipLo, double Aul )
341{
342 double Anp, bn, Bnp, delta_np;
343 double Eni, Enp;
344 double Gamma_np, p, n, g_p, g_n;
345 double ryd, s, kT_eV, rate, col_str, abs_osc_str;
346
347 DEBUG_ENTRY( "hydro_vs_deexcit()" );
348
349 kT_eV = EVRYD * phycon.te / TE1RYD;
350
351 /* This comes from equations 24 of Vriens and Smeets. */
352 /* >>refer he-like cs Vriens, L. \& Smeets, A. H. M. Phys. Rev. A, 22, 940 */
353 /* Computes the Vriens and Smeets
354 * rate coeff. for collisional dexcitation between any two levels of H.
355 * valid for all nhigh, nlow
356 * at end converts this excitation rate to collision strength */
357
358 n = (double)iso_sp[ipISO][nelem].st[ipLo].n();
359 p = (double)iso_sp[ipISO][nelem].st[ipHi].n();
360
361 ASSERT( n!=p );
362
363 g_n = iso_sp[ipISO][nelem].st[ipLo].g();
364 g_p = iso_sp[ipISO][nelem].st[ipHi].g();
365
366 ryd = EVRYD;
367 s = fabs(n-p);
368
369 Enp = EVRYD*(iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd - iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
370 Eni = EVRYD*iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd;
371
372 ASSERT( Enp > 0. );
373
374 /* This is an absorption oscillator strength. */
375 abs_osc_str = GetGF( Aul, Enp*RYD_INF/EVRYD, g_p)/g_n;
376 Anp = 2.*ryd/Enp*abs_osc_str;
377
378 bn = 1.4*log(n)/n - .7/n - .51/n/n + 1.16/n/n/n - 0.55/n/n/n/n;
379
380 Bnp = 4.*ryd*ryd/p/p/p*(1./Enp/Enp + 4./3.*Eni/POW3(Enp) + bn*Eni*Eni/powi(Enp,4));
381
382 delta_np = exp(-Bnp/Anp) + 0.1*Enp/ryd;
383
384 Gamma_np = ryd*log(1. + n*n*n*kT_eV/ryd) * (3. + 11.*s*s/n/n) /
385 ( 6. + 1.6*p*s + 0.3/pow2(s) + 0.8*sqrt(pow3(p))/sqrt(s)*fabs(s-0.6) );
386
387 if( 0.3*kT_eV/ryd+delta_np <= 0 )
388 {
389 rate = 0.;
390 }
391 else
392 {
393 rate = 1.6E-7 * sqrt(kT_eV) * g_n / g_p / ( kT_eV + Gamma_np ) *
394 ( Anp * log(0.3*kT_eV/ryd + delta_np) + Bnp );
395 }
396
397 col_str = rate / COLL_CONST * phycon.sqrte * iso_sp[ipISO][nelem].st[ipHi].g();
398
399 return( col_str );
400}
401
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define POW3
Definition cddefines.h:936
T pow2(T a)
Definition cddefines.h:931
#define MIN2
Definition cddefines.h:761
double ee1(double x)
Definition service.cpp:312
#define POW2
Definition cddefines.h:929
#define MAX2
Definition cddefines.h:782
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
T pow3(T a)
Definition cddefines.h:938
double powi(double, long int)
Definition service.cpp:604
double sum(double min, double max, Integrand func)
Definition cddefines.h:1531
@ ipELECTRON
Definition collision.h:9
t_dense dense
Definition dense.cpp:24
static const double ColliderMass[4]
Definition helike_cs.cpp:87
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
double Hion_coll_ioniz_ratecoef(long int ipISO, long int nelem, long int n, double ionization_energy_Ryd, double Te)
double hydro_vs_deexcit(long ipISO, long nelem, long ipHi, long ipLo, double Aul)
double CS_VS80(long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider)
double hydro_vs_ioniz(double ionization_energy_Ryd, double Te)
STATIC double hydro_vs_coll_str(double energy, long ipISO, long nelem, long ipHi, long ipLo, long Collider, double Aul)
double hydro_vs_coll_recomb(double ionization_energy_Ryd, double Te, double stat_level, double stat_ion)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
double ConvCrossSect2CollStr(double CrsSectCM2, double gLo, double E_ProjectileRyd, double reduced_mass_grams)
double GetGF(double trans_prob, double enercm, double gup)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOHR_RADIUS_CM
Definition physconst.h:222
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double RYD_INF
Definition physconst.h:115
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double PROTON_MASS
Definition physconst.h:94
UNUSED const double EVDEGK
Definition physconst.h:186
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double COLL_CONST
Definition physconst.h:229
UNUSED const double TE1RYD
Definition physconst.h:183
static double ** col_str
Definition species2.cpp:29
static double * g
Definition species2.cpp:28