cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /* CS_VS80 compute thermally averaged collision strength for collisional deexcitation 00004 * of hydrogenic atoms, from 00005 * >>refer H1 collision Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 00006 *hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients 00007 * for quantum number n */ 00008 /*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for ions 00009 * with all n, and Z*/ 00010 #include "cddefines.h" 00011 #include "phycon.h" 00012 #include "physconst.h" 00013 #include "iso.h" 00014 #include "hydro_vs_rates.h" 00015 #include "lines_service.h" 00016 00017 STATIC double hydro_vs_coll_str( double energy ); 00018 STATIC double Therm_ave_coll_str_int_VS80( double EOverKT ); 00019 00020 static long global_ipISO, global_nelem, global_ipHi, global_ipLo, global_Collider; 00021 static double global_temp, global_Aul; 00022 /* These are masses relative to the proton mass of the electron, proton, and alpha particle. */ 00023 static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0}; 00024 00025 /* 00026 Neither of these rates can be modified to account for impact by non-electrons because they 00027 are fits to thermally-averaged rates for electron impact...it can not be unravelled to 00028 expose a projectile energy that can then be scaled to account for other projectiles. 00029 Instead, we have included the original cross section formula (eq 14) from 00030 Vriens & Smeets (1980) below. 00031 */ 00032 00033 /* VS80 stands for Vriens and Smeets 1980 */ 00034 /* This routine calculates thermally-averaged collision strengths. */ 00035 double CS_VS80( long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider ) 00036 { 00037 double coll_str; 00038 00039 global_ipISO = ipISO; 00040 global_nelem = nelem; 00041 global_ipHi = ipHi; 00042 global_ipLo = ipLo; 00043 global_temp = temp; 00044 global_Collider = Collider; 00045 global_Aul = Aul; 00046 00047 /* evaluate collision strength, two options, 00048 * do thermal average (very slow) if set with 00049 * SET COLLISION STRENGTH AVERAGE command, 00050 * default just do point at kT 00051 * tests show no impact on test suite, much faster */ 00052 if( iso.lgCollStrenThermAver ) 00053 { 00054 /* do average over Maxwellian */ 00055 coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_VS80); 00056 coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VS80); 00057 } 00058 else 00059 { 00060 /* the argument to the function is equivalent to evaluating the function at 00061 * hnu = kT */ 00062 coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD ); 00063 } 00064 00065 ASSERT( coll_str >= 0. ); 00066 return coll_str; 00067 } 00068 00069 /* The integrand for calculating the thermal average of collision strengths */ 00070 STATIC double Therm_ave_coll_str_int_VS80( double EOverKT ) 00071 { 00072 double integrand; 00073 00074 integrand = exp( -1.*EOverKT ) * hydro_vs_coll_str( EOverKT * EVRYD*global_temp/TE1RYD ); 00075 00076 return integrand; 00077 } 00078 00079 /*hydro_vs_deexcit compute collision strength for collisional deexcitation for hydrogen atom, 00080 * from Vriens and Smeets */ 00081 STATIC double hydro_vs_coll_str( double energy ) 00082 { 00083 double Apn, bp, Bpn, delta; 00084 double Epi, Epn; 00085 double gamma, p, n; 00086 double ryd, s; 00087 double cross_section, coll_str, gLo, gHi, abs_osc_str, Aul; 00088 long ipISO, nelem, ipHi, ipLo, Collider; 00089 00090 DEBUG_ENTRY( "hydro_vs_coll_str()" ); 00091 00092 ipISO = global_ipISO; 00093 nelem = global_nelem; 00094 ipHi = global_ipHi; 00095 ipLo = global_ipLo; 00096 Collider = global_Collider; 00097 Aul = global_Aul; 00098 00099 gLo = StatesElem[ipISO][nelem][ipLo].g; 00100 gHi = StatesElem[ipISO][nelem][ipHi].g; 00101 00102 /* This comes from equations 14,15, and 16 of Vriens and Smeets. */ 00103 /* >>refer he-like cs Vriens, L. \& Smeets, A. H. M. Phys. Rev. A, 22, 940 */ 00104 /* Computes the Vriens and Smeets 00105 * rate coeff. for collisional dexcitation between any two levels of H. 00106 * valid for all nhigh, nlow 00107 * at end converts this excitation rate to collision strength */ 00108 00109 n = (double)StatesElem[ipISO][nelem][ipHi].n; 00110 p = (double)StatesElem[ipISO][nelem][ipLo].n; 00111 00112 ryd = EVRYD; 00113 s = fabs(n-p); 00114 00115 Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); 00116 Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]; 00117 00118 /* This is an absorption oscillator strength. */ 00119 abs_osc_str = GetGF( Aul, Epn*RYD_INF/EVRYD, gHi)/gLo; 00120 Apn = 2.*ryd/Epn*abs_osc_str; 00121 00122 bp = 1.4*log(p)/p - .7/p - .51/p/p + 1.16/p/p/p - 0.55/p/p/p/p; 00123 00124 Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4)); 00125 00126 delta = exp(-Bpn/Apn) - 0.4*Epn/ryd; 00127 00128 gamma = ryd*(8. + 23.*POW2(s/p))/ 00129 ( 8. + 1.1*n*s + 0.8/s/s + 0.4*(pow(n,1.5))/(pow(s,0.5))*fabs(s-1.0) ); 00130 00131 /* Scale the energy to get an equivalent electron energy. */ 00132 energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider]; 00133 00134 /* cross section in units of PI*a_o^2 */ 00135 if( energy/2./ryd+delta <= 0 /*|| energy < Epn*/ ) 00136 { 00137 cross_section = 0.; 00138 } 00139 else 00140 { 00141 cross_section = 2.*ryd/(energy + gamma)*(Apn*log(energy/2./ryd+delta) + Bpn); 00142 cross_section = MAX2( cross_section, 0. ); 00143 } 00144 00145 /* convert to collision strength */ 00146 coll_str = cross_section * gLo * (energy/EVRYD); 00147 00148 if( nelem==ipCARBON ) 00149 ASSERT( coll_str >= 0. ); 00150 00151 return( coll_str ); 00152 } 00153 00154 /*hydro_vs_ioniz generate hydrogenic collisional ionization rate coefficients 00155 * for quantum number n */ 00156 double hydro_vs_ioniz( long int ipISO, long int nelem, long int level ) 00157 { 00158 double coef, 00159 denom, 00160 epi, 00161 t_eV; 00162 00163 DEBUG_ENTRY( "hydro_vs_ioniz()" ); 00164 00165 ASSERT( nelem <= 1 ); 00166 00167 /* a function written to calculate the rate coefficients 00168 * for hydrogen collisional ionizations from 00169 * Jason Ferguson, summer 94 00170 * valid for all n 00171 * >>refer H1 collision Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 00172 * */ 00173 00174 /* this is kT is eV */ 00175 t_eV = phycon.te/EVDEGK; 00176 00177 /* this is the ionization energy relative to kT, dimensionless */ 00178 epi = (iso.xIsoLevNIonRyd[ipISO][nelem][level] * EVRYD) / t_eV; 00179 00180 /* this is the denominator of equation 8 of VS80. */ 00181 denom = pow(epi,2.33) + 4.38*pow(epi,1.72) + 1.32*epi; 00182 00183 /* this is equation 8 of VS80 */ 00184 coef = 9.56e-6 * pow(t_eV,-1.5) * StatesElem[ipISO][nelem][level].ConBoltz / denom; 00185 00186 ASSERT( coef >= 0. ); 00187 return( coef ); 00188 } 00189 00190 /*Hion_coll_ioniz_ratecoef calculate hydrogenic ionization rates for all n, and Z*/ 00191 double Hion_coll_ioniz_ratecoef( 00192 /* the isoelectronic sequence */ 00193 long int ipISO , 00194 /* element, >=1 since only used for ions 00195 * nelem = 1 is helium the least possible charge */ 00196 long int nelem, 00197 /* principal quantum number, > 1 00198 * since only used for excited states */ 00199 long int level) 00200 { 00201 long int n; 00202 double H, 00203 HydColIon_v, 00204 Rnp, 00205 charge, 00206 chim, 00207 eone, 00208 etwo, 00209 ethree, 00210 g, 00211 rate, 00212 rate2, 00213 /* ryd, */ 00214 t1, 00215 t2, 00216 t3, 00217 t4, 00218 tev, 00219 xn, 00220 y; 00221 static const double arrH[4] = {1.48,3.64,5.93,8.32}; 00222 static const double arrRnp[8] = {2.20,1.90,1.73,1.65,1.60,1.56,1.54,1.52}; 00223 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}; 00224 00225 static double small = 0.; 00226 00227 DEBUG_ENTRY( "Hion_coll_ioniz_ratecoef()" ); 00228 /*calculate hydrogenic ionization rates for all n, and Z 00229 * >>refer HI cs Allen 1973, Astro. Quan. for low Te. 00230 * >>refer HI cs Sampson and Zhang 1988, ApJ, 335, 516 for High Te. 00231 * */ 00232 00234 charge = nelem - ipISO; 00235 /* this routine only for ions, nelem=0 is H, nelem=1 he, etc */ 00236 ASSERT( charge > 0); 00237 /* this is quantum level, not for ground state (n=1), only 00238 * n=2 or higher */ 00239 n = StatesElem[ipISO][nelem][level].n; 00240 ASSERT( n>1 ); 00241 00242 if( n > 4 ) 00243 { 00244 H = 2.15*n; 00245 } 00246 else 00247 { 00248 H = arrH[n-1]; 00249 } 00250 00251 if( n > 8 ) 00252 { 00253 Rnp = 1.52; 00254 } 00255 else 00256 { 00257 Rnp = arrRnp[n-1]; 00258 } 00259 00260 if( n > 10 ) 00261 { 00262 g = arrg[9]; 00263 } 00264 else 00265 { 00266 g = arrg[n-1]; 00267 } 00268 00269 tev = EVRYD/TE1RYD*phycon.te; 00270 /* ryd = EVRYD; */ 00271 xn = (double)n; 00272 /* >>chng 19 dec 02, use proper energy! */ 00273 /*chim = POW2(charge+1.)*ryd/xn/xn;*/ 00274 chim = EVRYD * iso.xIsoLevNIonRyd[ipISO][nelem][level]; 00275 y = chim/tev; 00276 00277 eone = ee1(y); 00278 etwo = StatesElem[ipISO][nelem][level].ConBoltz - y*eone; 00279 ethree = (StatesElem[ipISO][nelem][level].ConBoltz - y*etwo)/2.; 00280 00281 t1 = 1/xn*eone; 00282 t2 = 1./3./xn*(StatesElem[ipISO][nelem][level].ConBoltz - y*ethree); 00283 t3 = 3.*H/xn/(3. - Rnp)*(y*etwo - 2.*y*eone + StatesElem[ipISO][nelem][level].ConBoltz); 00284 t4 = 3.36*y*(eone - etwo); 00285 rate = 7.69415e-9*(double)phycon.sqrte*9.28278e-3*powi(xn/(charge+1),4)*g*y; 00286 rate *= t1 - t2 + t3 + t4; 00288 rate = MAX2(rate,small); 00289 00290 rate2 = 2.1e-8*phycon.sqrte/chim/chim*dsexp(2.302585*5040.* 00291 chim/(double)phycon.te); 00292 00293 rate2 = MAX2(rate2,small); 00294 00295 /* Take the lowest of the two, they fit nicely together... */ 00296 /* >>chng 10 Sept 02, sometimes one of these is zero and the other is positive. 00297 * in that case take the bigger one. */ 00298 if( rate==0. || rate2==0. ) 00299 HydColIon_v = MAX2(rate,rate2); 00300 else 00301 HydColIon_v = MIN2(rate,rate2); 00302 00303 ASSERT( HydColIon_v >= 0. ); 00304 return( HydColIon_v ); 00305 } 00306 00307 /*hydro_vs_deexcit compute collisional deexcitation rates for hydrogen atom, 00308 * from Vriens and Smeets 1980 */ 00309 double hydro_vs_deexcit( long ipISO, long nelem, long ipHi, long ipLo, double Aul ) 00310 { 00311 double Anp, bn, Bnp, delta_np; 00312 double Eni, Enp; 00313 double Gamma_np, p, n, g_p, g_n; 00314 double ryd, s, kT_eV, rate, col_str, abs_osc_str; 00315 00316 DEBUG_ENTRY( "hydro_vs_coll_str()" ); 00317 00318 kT_eV = EVRYD * phycon.te / TE1RYD; 00319 00320 /* This comes from equations 24 of Vriens and Smeets. */ 00321 /* >>refer he-like cs Vriens, L. \& Smeets, A. H. M. Phys. Rev. A, 22, 940 */ 00322 /* Computes the Vriens and Smeets 00323 * rate coeff. for collisional dexcitation between any two levels of H. 00324 * valid for all nhigh, nlow 00325 * at end converts this excitation rate to collision strength */ 00326 00327 n = (double)StatesElem[ipISO][nelem][ipLo].n; 00328 p = (double)StatesElem[ipISO][nelem][ipHi].n; 00329 00330 ASSERT( n!=p ); 00331 00332 g_n = StatesElem[ipISO][nelem][ipLo].g; 00333 g_p = StatesElem[ipISO][nelem][ipHi].g; 00334 00335 ryd = EVRYD; 00336 s = fabs(n-p); 00337 00338 Enp = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); 00339 Eni = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]; 00340 00341 ASSERT( Enp > 0. ); 00342 00343 /* This is an absorption oscillator strength. */ 00344 abs_osc_str = GetGF( Aul, Enp*RYD_INF/EVRYD, g_p)/g_n; 00345 Anp = 2.*ryd/Enp*abs_osc_str; 00346 00347 bn = 1.4*log(n)/n - .7/n - .51/n/n + 1.16/n/n/n - 0.55/n/n/n/n; 00348 00349 Bnp = 4.*ryd*ryd/p/p/p*(1./Enp/Enp + 4./3.*Eni/POW3(Enp) + bn*Eni*Eni/powi(Enp,4)); 00350 00351 delta_np = exp(-Bnp/Anp) + 0.1*Enp/ryd; 00352 00353 Gamma_np = ryd*log(1. + n*n*n*kT_eV/ryd) * (3. + 11.*s*s/n/n) / 00354 ( 6. + 1.6*p*s + 0.3/s/s + 0.8*(pow(p,1.5))/(pow(s,0.5))*fabs(s-0.6) ); 00355 00356 if( 0.3*kT_eV/ryd+delta_np <= 0 ) 00357 { 00358 rate = 0.; 00359 } 00360 else 00361 { 00362 rate = 1.6E-7 * pow(kT_eV, 0.5) * g_n / g_p / ( kT_eV + Gamma_np ) * 00363 ( Anp * log(0.3*kT_eV/ryd + delta_np) + Bnp ); 00364 } 00365 00366 col_str = rate / COLL_CONST * phycon.sqrte * StatesElem[ipISO][nelem][ipHi].g; 00367 00368 return( col_str ); 00369 }