00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "dense.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 #include "taulines.h"
00017
00018 STATIC double hydro_vs_coll_str( double energy, long ipISO, long nelem, long ipHi, long ipLo, long Collider, double Aul );
00019
00020 namespace {
00021 class my_Integrand : public std::unary_function<double, double>
00022 {
00023 public:
00024 long ipISO, nelem, ipHi, ipLo, Collider;
00025 double Aul;
00026 double temp;
00027
00028 double operator() (double EOverKT)
00029 {
00030 double col_str = hydro_vs_coll_str( EOverKT * EVRYD * temp/TE1RYD, ipISO, nelem, ipHi, ipLo, Collider, Aul );
00031 return exp( -1.*EOverKT ) * col_str;
00032 }
00033 };
00034 }
00035
00036
00037 static const double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 double CS_VS80( long int ipISO, long int nelem, long int ipHi, long int ipLo, double Aul, double temp, long int Collider )
00050 {
00051 double coll_str;
00052
00053 if( Collider == ipELECTRON )
00054 {
00055 coll_str = hydro_vs_deexcit( ipISO, nelem, ipHi, ipLo, Aul);
00056 }
00057 else
00058 {
00059
00060
00061
00062
00063
00064 if( iso.lgCollStrenThermAver )
00065 {
00066 my_Integrand func;
00067
00068 func.ipISO = ipISO;
00069 func.nelem = nelem;
00070 func.ipHi = ipHi;
00071 func.ipLo = ipLo;
00072 func.temp = temp;
00073 func.Collider = Collider;
00074 func.Aul = Aul;
00075
00076 Integrator<my_Integrand,Gaussian32> VS80;
00077
00078 coll_str = VS80.sum( 0., 1., func );
00079 coll_str += VS80.sum( 1., 10., func );
00080 }
00081 else
00082 {
00083
00084
00085 coll_str = hydro_vs_coll_str( EVRYD*temp/TE1RYD, ipISO, nelem, ipHi, ipLo, Collider, Aul );
00086 }
00087 }
00088
00089 ASSERT( coll_str >= 0. );
00090 return coll_str;
00091 }
00092
00093
00094
00095 STATIC double hydro_vs_coll_str( double energy, long ipISO, long nelem, long ipHi, long ipLo, long Collider, double Aul )
00096 {
00097 double Apn, bp, Bpn, delta;
00098 double Epi, Epn;
00099 double gamma, p, n;
00100 double ryd, s;
00101 double cross_section, coll_str, gLo, gHi, abs_osc_str, reduced_mass;
00102
00103 DEBUG_ENTRY( "hydro_vs_coll_str()" );
00104
00105
00106 ASSERT( Collider >= 0.&& Collider <4 );
00107 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
00108 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
00109
00110 gLo = StatesElemNEW[nelem][nelem-ipISO][ipLo].g;
00111 gHi = StatesElemNEW[nelem][nelem-ipISO][ipHi].g;
00112
00113
00114
00115
00116
00117
00118
00119
00120 n = (double)StatesElemNEW[nelem][nelem-ipISO][ipHi].n;
00121 p = (double)StatesElemNEW[nelem][nelem-ipISO][ipLo].n;
00122
00123 ryd = EVRYD;
00124 s = fabs(n-p);
00125 ASSERT( s > 0. );
00126
00127 Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
00128 Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];
00129
00130
00131 abs_osc_str = GetGF( Aul, Epn*RYD_INF/EVRYD, gHi)/gLo;
00132 Apn = 2.*ryd/Epn*abs_osc_str;
00133
00134 bp = 1.4*log(p)/p - .7/p - .51/p/p + 1.16/p/p/p - 0.55/p/p/p/p;
00135
00136 Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));
00137
00138 delta = exp(-Bpn/Apn) - 0.4*Epn/ryd;
00139
00140 gamma = ryd*(8. + 23.*POW2(s/p))/
00141 ( 8. + 1.1*n*s + 0.8/pow2(s) + 0.4*sqrt(pow3(n))/sqrt(s)*fabs(s-1.0) );
00142
00143
00144 energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
00145
00146
00147 if( energy/2./ryd+delta <= 0 )
00148 {
00149 cross_section = 0.;
00150 }
00151 else
00152 {
00153 cross_section = 2.*ryd/(energy + gamma)*(Apn*log(energy/2./ryd+delta) + Bpn);
00154 cross_section = MAX2( cross_section, 0. );
00155 }
00156
00157
00158 coll_str = ConvCrossSect2CollStr( cross_section * PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM, gLo, energy/EVRYD, reduced_mass );
00159
00160 ASSERT( coll_str >= 0. );
00161
00162 return( coll_str );
00163 }
00164
00165
00166 double hydro_vs_coll_recomb( double ionization_energy_Ryd, double Te, double stat_level, double stat_ion )
00167 {
00168 double coef,
00169 denom,
00170 epi,
00171 t_eV;
00172
00173 DEBUG_ENTRY( "hydro_vs_coll_recomb()" );
00174
00175
00176
00177
00178
00179 t_eV = Te/EVDEGK;
00180
00181
00182 epi = ionization_energy_Ryd * EVRYD / t_eV;
00183
00184
00185 denom = pow(epi,2.33) + 4.38*pow(epi,1.72) + 1.32*epi;
00186
00187
00188 coef = 3.17e-27 / pow3(t_eV) * stat_level / stat_ion / denom;
00189
00190 ASSERT( coef >= 0. );
00191 return( coef );
00192 }
00193
00194
00195
00196 double hydro_vs_ioniz( double ionization_energy_Ryd, double Te )
00197 {
00198 double coef,
00199 denom,
00200 epi,
00201 t_eV;
00202
00203 DEBUG_ENTRY( "hydro_vs_ioniz()" );
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 t_eV = Te/EVDEGK;
00214
00215
00216 epi = ionization_energy_Ryd * EVRYD / t_eV;
00217
00218
00219 denom = pow(epi,2.33) + 4.38*pow(epi,1.72) + 1.32*epi;
00220
00221
00222 coef = 9.56e-6 / sqrt(pow3(t_eV)) * dsexp( epi ) / denom;
00223
00224 ASSERT( coef >= 0. );
00225 return( coef );
00226 }
00227
00228
00229 double Hion_coll_ioniz_ratecoef(
00230
00231 long int ipISO ,
00232
00233
00234 long int nelem,
00235
00236
00237 long int n,
00238 double ionization_energy_Ryd,
00239 double Te )
00240 {
00241 double H,
00242 HydColIon_v,
00243 Rnp,
00244 charge,
00245 chim,
00246 eone,
00247 etwo,
00248 ethree,
00249 g,
00250 rate,
00251 rate2,
00252 boltz,
00253 t1,
00254 t2,
00255 t3,
00256 t4,
00257 tev,
00258 xn,
00259 y;
00260 static const double arrH[4] = {1.48,3.64,5.93,8.32};
00261 static const double arrRnp[8] = {2.20,1.90,1.73,1.65,1.60,1.56,1.54,1.52};
00262 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};
00263
00264 static double small = 0.;
00265
00266 DEBUG_ENTRY( "Hion_coll_ioniz_ratecoef()" );
00267
00268
00269
00270
00271
00272 charge = nelem - ipISO;
00273
00274 ASSERT( charge > 0);
00275 ASSERT( n>1 );
00276
00277 if( n > 4 )
00278 {
00279 H = 2.15*n;
00280 }
00281 else
00282 {
00283 H = arrH[n-1];
00284 }
00285
00286 if( n > 8 )
00287 {
00288 Rnp = 1.52;
00289 }
00290 else
00291 {
00292 Rnp = arrRnp[n-1];
00293 }
00294
00295 if( n > 10 )
00296 {
00297 g = arrg[9];
00298 }
00299 else
00300 {
00301 g = arrg[n-1];
00302 }
00303
00304 tev = EVRYD/TE1RYD*Te;
00305 xn = (double)n;
00306 chim = EVRYD * ionization_energy_Ryd;
00307 y = chim/tev;
00308 boltz = dsexp( chim/tev );
00309
00310 eone = ee1(y);
00311 etwo = boltz - y*eone;
00312 ethree = (boltz - y*etwo)/2.;
00313
00314 t1 = 1/xn*eone;
00315 t2 = 1./3./xn*(boltz - y*ethree);
00316 t3 = 3.*H/xn/(3. - Rnp)*(y*etwo - 2.*y*eone + boltz);
00317 t4 = 3.36*y*(eone - etwo);
00318 rate = 7.69415e-9*sqrt(Te)*9.28278e-3*powi(xn/(charge+1),4)*g*y;
00319 rate *= t1 - t2 + t3 + t4;
00320 rate2 = 2.1e-8*sqrt(Te)/chim/chim*dsexp(2.302585*5040.*chim/Te);
00321
00322
00323 rate = MAX2(rate,small);
00324 rate2 = MAX2(rate2,small);
00325
00326
00327
00328
00329 if( rate==0. || rate2==0. )
00330 HydColIon_v = MAX2(rate,rate2);
00331 else
00332 HydColIon_v = MIN2(rate,rate2);
00333
00334 ASSERT( HydColIon_v >= 0. );
00335 return( HydColIon_v );
00336 }
00337
00338
00339
00340 double hydro_vs_deexcit( long ipISO, long nelem, long ipHi, long ipLo, double Aul )
00341 {
00342 double Anp, bn, Bnp, delta_np;
00343 double Eni, Enp;
00344 double Gamma_np, p, n, g_p, g_n;
00345 double ryd, s, kT_eV, rate, col_str, abs_osc_str;
00346
00347 DEBUG_ENTRY( "hydro_vs_deexcit()" );
00348
00349 kT_eV = EVRYD * phycon.te / TE1RYD;
00350
00351
00352
00353
00354
00355
00356
00357
00358 n = (double)StatesElemNEW[nelem][nelem-ipISO][ipLo].n;
00359 p = (double)StatesElemNEW[nelem][nelem-ipISO][ipHi].n;
00360
00361 ASSERT( n!=p );
00362
00363 g_n = StatesElemNEW[nelem][nelem-ipISO][ipLo].g;
00364 g_p = StatesElemNEW[nelem][nelem-ipISO][ipHi].g;
00365
00366 ryd = EVRYD;
00367 s = fabs(n-p);
00368
00369 Enp = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
00370 Eni = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipHi];
00371
00372 ASSERT( Enp > 0. );
00373
00374
00375 abs_osc_str = GetGF( Aul, Enp*RYD_INF/EVRYD, g_p)/g_n;
00376 Anp = 2.*ryd/Enp*abs_osc_str;
00377
00378 bn = 1.4*log(n)/n - .7/n - .51/n/n + 1.16/n/n/n - 0.55/n/n/n/n;
00379
00380 Bnp = 4.*ryd*ryd/p/p/p*(1./Enp/Enp + 4./3.*Eni/POW3(Enp) + bn*Eni*Eni/powi(Enp,4));
00381
00382 delta_np = exp(-Bnp/Anp) + 0.1*Enp/ryd;
00383
00384 Gamma_np = ryd*log(1. + n*n*n*kT_eV/ryd) * (3. + 11.*s*s/n/n) /
00385 ( 6. + 1.6*p*s + 0.3/pow2(s) + 0.8*sqrt(pow3(p))/sqrt(s)*fabs(s-0.6) );
00386
00387 if( 0.3*kT_eV/ryd+delta_np <= 0 )
00388 {
00389 rate = 0.;
00390 }
00391 else
00392 {
00393 rate = 1.6E-7 * sqrt(kT_eV) * g_n / g_p / ( kT_eV + Gamma_np ) *
00394 ( Anp * log(0.3*kT_eV/ryd + delta_np) + Bnp );
00395 }
00396
00397 col_str = rate / COLL_CONST * phycon.sqrte * StatesElemNEW[nelem][nelem-ipISO][ipHi].g;
00398
00399 return( col_str );
00400 }
00401