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

Generated for cloudy by doxygen 1.8.3.1