00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "phycon.h"
00012 #include "physconst.h"
00013 #include "hydrooscilstr.h"
00014 #include "iso.h"
00015 #include "hydro_vs_rates.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;
00022
00023
00024 static double ColliderMass[3] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0};
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 double hydro_vs_deexcit(
00037 long int ipISO,
00038 long int nelem,
00039 long int ipHi,
00040 long int ipLo
00041
00042 )
00043 {
00044 double Anp,
00045 bn,
00046 Bnp,
00047 delta,
00048
00049
00050 Eni,
00051 Enp,
00052 gamma,
00053 hdexct_v,
00054 xlow,
00055 p,
00056 rate,
00057 ryd,
00058 s,
00059 tev;
00060
00061 DEBUG_ENTRY( "hydro_vs_deexcit()" );
00062
00063
00064
00065
00066
00067
00068
00069 p = (double)iso.quant_desig[ipISO][nelem][ipHi].n;
00070
00071 xlow = (double)iso.quant_desig[ipISO][nelem][ipLo].n;
00072 tev = EVRYD/TE1RYD*phycon.te;
00073 ryd = EVRYD;
00074 s = fabs(p-xlow);
00075
00076
00077
00078
00079 Enp = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
00080 Eni = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 Anp = 2./Enp*HydroOscilStr(xlow,p);
00091
00092 bn = 1.4*log(xlow)/xlow - .7/xlow - .51/xlow/xlow + 1.16/xlow/xlow/xlow - 0.55/xlow/xlow/xlow/xlow;
00093
00094 Bnp = 4.*ryd*ryd/p/p/p*(1./Enp/Enp + 4./3.*Eni/POW3(Enp) + bn*Eni*Eni/powi(Enp,4));
00095
00096 delta = exp(-Bnp/Anp) + 0.06*s*s/p/xlow/xlow;
00097
00098 gamma = ryd*log(1.+POW3(xlow)*tev/ryd)*(3. + 11.*POW2(s/xlow))/
00099 (6. + 1.6*p*s + .3/s/s + 0.8*(pow(p,1.5))/(pow(s,.5))*fabs(s-0.6));
00100 rate = (1.6e-7)*sqrt(tev)/(tev + gamma)*(Anp*log(0.3*tev/ryd+delta) + Bnp);
00101
00102 rate *= iso.stat[ipISO][nelem][ipLo]/iso.stat[ipISO][nelem][ipHi];
00103
00104
00105 hdexct_v = rate;
00106
00107
00108
00109
00110 hdexct_v = iso.stat[ipISO][nelem][ipHi]*rate*phycon.sqrte/COLL_CONST;
00111
00112 ASSERT( hdexct_v >= 0. );
00113
00114 DEBUG_EXIT( "hydro_vs_deexcit()" );
00115 return( hdexct_v );
00116 }
00117
00118
00119
00120 double hydro_vs_excit(
00121 long int ipISO,
00122 long int nelem,
00123 long int ipHi,
00124 long int ipLo
00125
00126 )
00127 {
00128 double Apn,
00129 bp,
00130 Bpn,
00131 delta,
00132 epn,
00133 Epi,
00134 Epn,
00135 gamma,
00136
00137 xlow,
00138 n,
00139 rate,
00140 ryd,
00141 s,
00142 tev;
00143
00144 DEBUG_ENTRY( "hydro_vs_excit()" );
00145
00146
00147
00148
00149
00150
00151
00152 n = (double)iso.quant_desig[ipISO][nelem][ipHi].n;
00153
00154 xlow = (double)iso.quant_desig[ipISO][nelem][ipLo].n;
00155
00156 tev = EVRYD/TE1RYD*phycon.te;
00157 ryd = EVRYD;
00158 s = fabs(n-xlow);
00159
00160
00161
00162
00163 Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
00164 Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];
00165
00166 epn = Epn / tev;
00167
00168
00169 Apn = 2./Epn*HydroOscilStr(xlow,n);
00170
00171 bp = 1.4*log(xlow)/xlow - .7/xlow - .51/xlow/xlow + 1.16/xlow/xlow/xlow - 0.55/xlow/xlow/xlow/xlow;
00172
00173 Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));
00174
00175 delta = exp(-Bpn/Apn) + 0.06*s*s/n/xlow/xlow;
00176
00177 gamma = ryd*log(1.+POW3(xlow)*tev/ryd)*(3. + 11.*POW2(s/xlow))/
00178 (6. + 1.6*n*s + .3/s/s + 0.8*(pow(n,1.5))/(pow(s,.5))*fabs(s-0.6));
00179
00180 rate = (1.6e-7)*sqrt(tev)/(tev + gamma)*exp(-epn)*(Apn*log(0.3*tev/ryd+delta) +
00181 Bpn);
00182
00183
00184
00185
00186
00187 ASSERT( rate >= 0. );
00188
00189 DEBUG_EXIT( "hydro_vs_excit()" );
00190 return( rate );
00191 }
00192
00193
00194
00195
00196 double CS_VS80( long int ipISO, long int nelem, long int ipHi, long int ipLo, double temp, long int Collider )
00197 {
00198 double coll_str;
00199
00200 global_ipISO = ipISO;
00201 global_nelem = nelem;
00202 global_ipHi = ipHi;
00203 global_ipLo = ipLo;
00204 global_temp = temp;
00205 global_Collider = Collider;
00206
00207 #if 0
00208 if( helike.lgCS_therm_ave == false )
00209 {
00210
00211
00212
00213
00214 if( (dense.eden > 10000.) && (dense.eden < 1E10 ) )
00215 {
00216 coll_str = qg32( 0.0, 6.0, Therm_ave_coll_str_int_VS80);
00217 }
00218 else
00219 {
00220
00221 coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD, false );
00222 }
00223 }
00224 else
00225 #endif
00226
00227
00228
00229
00230
00231
00232 if( iso.lgCollStrenThermAver )
00233 {
00234
00235 coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_VS80);
00236 coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VS80);
00237 }
00238 else
00239 {
00240
00241
00242 coll_str = hydro_vs_coll_str( EVRYD*global_temp/TE1RYD );
00243 }
00244
00245 ASSERT( coll_str >= 0. );
00246 return coll_str;
00247 }
00248
00249
00250 static double Therm_ave_coll_str_int_VS80( double EOverKT )
00251 {
00252 double integrand;
00253
00254 integrand = exp( -1.*EOverKT ) * hydro_vs_coll_str( EOverKT * EVRYD*global_temp/TE1RYD );
00255
00256 return integrand;
00257 }
00258
00259
00260
00261 static double hydro_vs_coll_str( double energy )
00262 {
00263 double Apn, bp, Bpn, delta;
00264 double Epi, Epn;
00265 double gamma, p, n;
00266 double ryd, s ;
00267 double cross_section, coll_str, stat_weight;
00268 long ipISO, nelem, ipHi, ipLo, Collider;
00269
00270 DEBUG_ENTRY( "hydro_vs_coll_str()" );
00271
00272 ipISO = global_ipISO;
00273 nelem = global_nelem;
00274 ipHi = global_ipHi;
00275 ipLo = global_ipLo;
00276
00277 Collider = global_Collider;
00278
00279 stat_weight = iso.stat[ipISO][nelem][ipLo];
00280
00281
00282
00283
00284
00285
00286
00287
00289
00290 n = (double)iso.quant_desig[ipISO][nelem][ipHi].n;
00291
00292 p = (double)iso.quant_desig[ipISO][nelem][ipLo].n;
00293
00294
00295 ryd = EVRYD;
00296 s = fabs(n-p);
00297
00298 Epn = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
00299 Epi = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][ipLo];
00300
00301
00302
00303
00304
00305
00306
00307 Apn = 2.*ryd/Epn*HydroOscilStr(p,n);
00308
00309 bp = 1.4*log(p)/p - .7/p - .51/p/p + 1.16/p/p/p - 0.55/p/p/p/p;
00310
00311 Bpn = 4.*ryd*ryd/n/n/n*(1./Epn/Epn + 4./3.*Epi/POW3(Epn) + bp*Epi*Epi/powi(Epn,4));
00312
00313 delta = exp(-Bpn/Apn) - 0.4*Epn/ryd;
00314
00315 gamma = ryd*(8. + 23.*POW2(s/p))/
00316 (8. + 1.1*n*s + 0.8/s/s + 0.4*(pow(n,1.5))/(pow(s,0.5))*fabs(s-1.0));
00317
00318
00319 energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
00320
00321
00322 if( energy/2/ryd+delta <= 0 )
00323 {
00324 cross_section = 0.;
00325 }
00326 else
00327 {
00328 cross_section = 2.*ryd/(energy + gamma)*(Apn*log(energy/2/ryd+delta) + Bpn);
00329 cross_section = MAX2( cross_section, 0. );
00330 }
00331
00332
00333 coll_str = cross_section * POW2( BOHR_RADIUS_CM*energy/EVRYD*RYD_INF ) * stat_weight;
00334
00335
00336 if( nelem==ipCARBON )
00337 ASSERT( coll_str >= 0. );
00338 DEBUG_EXIT( "hydro_vs_coll_str()" );
00339
00340 return( coll_str );
00341 }
00342
00343
00344
00345 double hydro_vs_ioniz( long int ipISO, long int nelem, long int level )
00346 {
00347 double coef,
00348 denom,
00349 eps,
00350 smeets_v,
00351 tryd;
00352
00353 DEBUG_ENTRY( "hydro_vs_ioniz()" );
00354
00355 ASSERT( nelem <= 1 );
00356
00357
00358
00359
00360
00361
00362
00363 tryd = EVRYD/TE1RYD*phycon.te;
00364 eps = TE1RYD/phycon.te*iso.xIsoLevNIonRyd[ipISO][nelem][level];
00365 denom = pow(eps,2.33) + 4.38*pow(eps,1.72) + 1.32*eps;
00366
00367 coef = 9.56e-6/pow(tryd,1.5)/denom*iso.ConBoltz[ipISO][nelem][level];
00368
00369 smeets_v = coef;
00370
00371 ASSERT( smeets_v >= 0. );
00372
00373 DEBUG_EXIT( "hydro_vs_ioniz()" );
00374 return( smeets_v );
00375 }
00376
00377
00378 double Hion_coll_ioniz_ratecoef(
00379
00380 long int ipISO ,
00381
00382
00383 long int nelem,
00384
00385
00386 long int level)
00387 {
00388 long int n;
00389 double H,
00390 HydColIon_v,
00391 Rnp,
00392 charge,
00393 chim,
00394 eone,
00395 etwo,
00396 ethree,
00397 g,
00398 rate,
00399 rate2,
00400
00401 t1,
00402 t2,
00403 t3,
00404 t4,
00405 tev,
00406 xn,
00407 y;
00408 static const double arrH[4] = {1.48,3.64,5.93,8.32};
00409 static const double arrRnp[8] = {2.20,1.90,1.73,1.65,1.60,1.56,1.54,1.52};
00410 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};
00411
00412 static double small = 0.;
00413
00414 DEBUG_ENTRY( "Hion_coll_ioniz_ratecoef()" );
00415
00416
00417
00418
00419
00421 charge = nelem - ipISO;
00422
00423 ASSERT( charge > 0);
00424
00425
00426 n = iso.quant_desig[ipISO][nelem][level].n;
00427 ASSERT( n>1 );
00428
00429 if( n > 4 )
00430 {
00431 H = 2.15*n;
00432 }
00433 else
00434 {
00435 H = arrH[n-1];
00436 }
00437
00438 if( n > 8 )
00439 {
00440 Rnp = 1.52;
00441 }
00442 else
00443 {
00444 Rnp = arrRnp[n-1];
00445 }
00446
00447 if( n > 10 )
00448 {
00449 g = arrg[9];
00450 }
00451 else
00452 {
00453 g = arrg[n-1];
00454 }
00455
00456 tev = EVRYD/TE1RYD*phycon.te;
00457
00458 xn = (double)n;
00459
00460
00461 chim = EVRYD * iso.xIsoLevNIonRyd[ipISO][nelem][level];
00462 y = chim/tev;
00463
00464 eone = ee1(y);
00465 etwo = iso.ConBoltz[ipISO][nelem][level] - y*eone;
00466 ethree = (iso.ConBoltz[ipISO][nelem][level] - y*etwo)/2.;
00467
00468 t1 = 1/xn*eone;
00469 t2 = 1./3./xn*(iso.ConBoltz[ipISO][nelem][level] - y*ethree);
00470 t3 = 3.*H/xn/(3. - Rnp)*(y*etwo - 2.*y*eone + iso.ConBoltz[ipISO][nelem][level]);
00471 t4 = 3.36*y*(eone - etwo);
00472 rate = 7.69415e-9*(double)phycon.sqrte*9.28278e-3*powi(xn/(charge+1),4)*g*y;
00473 rate *= t1 - t2 + t3 + t4;
00475 rate = MAX2(rate,small);
00476
00477 rate2 = 2.1e-8*phycon.sqrte/chim/chim*dsexp(2.302585*5040.*
00478 chim/(double)phycon.te);
00479
00480 rate2 = MAX2(rate2,small);
00481
00482
00483
00484
00485 if( rate==0. || rate2==0. )
00486 HydColIon_v = MAX2(rate,rate2);
00487 else
00488 HydColIon_v = MIN2(rate,rate2);
00489
00490 ASSERT( HydColIon_v >= 0. );
00491
00492 DEBUG_EXIT( "Hion_coll_ioniz_ratecoef()" );
00493 return( HydColIon_v );
00494 }
00495
00496
00497
00498 double Hion_colldeexc_cs(long int ipHi,
00499 long int ipLo,
00500
00501 long int nelem,
00502 long int ipISO )
00503 {
00504 long int nHi,
00505 nLo;
00506 double A,
00507 An,
00508 D,
00509 H,
00510 HydColDwn_v,
00511 R,
00512 c,
00513 charge,
00514 e1,
00515 exp_int_2,
00516 gHi,
00517 gLo,
00518 rate,
00519 ryd,
00520 sexpy,
00521 tev,
00522 xn,
00523 xnp2,
00524 xp,
00525 y;
00526 static const double arrAn[10] = {1.30,0.59,0.38,0.286,0.229,0.192,0.164,0.141,0.121,0.105};
00527 static const double arrH[4] = {1.48,3.64,5.93,8.32};
00528 static const double arrR[10] = {1.83,1.60,1.53,1.495,1.475,1.46,1.45,1.45,1.46,1.47};
00529
00530 DEBUG_ENTRY( "Hion_colldeexc_cs()" );
00531
00532
00533
00534
00535
00536
00537
00538 charge = nelem + 1 - ipISO;
00539
00540
00541 if( (charge < 0 || ipLo >= ipHi) || ipLo*ipHi < 0 )
00542 {
00543 fprintf( ioQQQ, "Hion_colldeexc_cs called with impossible parameters.\n" );
00544 fprintf( ioQQQ, "%ld %ld %f\n", ipHi, ipLo, charge );
00545 puts( "[Stop in Hion_colldeexc_cs]" );
00546 cdEXIT(EXIT_FAILURE);
00547 }
00548
00549
00550
00551
00552 nHi = ipHi;
00553 gHi = 2.*nHi*nHi;
00554 if( ipLo == 0 )
00555 {
00556 nLo = 1;
00557 gLo = 2.;
00558 }
00559 else if( ipLo == 1 )
00560 {
00561 nLo = 2;
00562 gLo = 2.;
00563 }
00564 else if( ipLo == 2 )
00565 {
00566 nLo = 2;
00567 gLo = 6.;
00568 }
00569 else
00570 {
00571 nLo = ipLo;
00572 gLo = 2.*nLo*nLo;
00573 }
00574
00575
00576 if( charge == 1 )
00577 {
00578
00579
00580 HydColDwn_v = hydro_vs_deexcit(ipISO, nelem, ipHi, ipLo);
00581
00582 DEBUG_EXIT( "Hion_colldeexc_cs()" );
00583 return( HydColDwn_v );
00584 }
00585
00586 if( nLo > 4 )
00587 {
00588 H = 2.15*nLo;
00589 }
00590 else
00591 {
00592 H = arrH[nLo-1];
00593 }
00594 if( nLo > 10 )
00595 {
00596 An = 1./(float)(nLo);
00597 R = 1.50;
00598 }
00599 else
00600 {
00601 An = arrAn[nLo-1];
00602 R = arrR[nLo-1];
00603 }
00604
00605 xp = (double)(nHi);
00606 xn = (double)(nLo);
00607 xnp2 = POW2(xn/xp);
00608 tev = EVRYD/TE1RYD*phycon.te;
00609 ryd = EVRYD;
00610
00611 y = POW2(charge)*ryd*(1./POW2(xn) - 1./POW2(xp))/tev;
00612
00613 sexpy = dsexp(y);
00614
00615
00616 if( sexpy <= 0. )
00617 {
00618
00619 e1 = ee1_safe(y);
00620
00621 A = 4.*powi(xn,4)/(1 - xnp2)*HydroOscilStr(xn,xp);
00622 D = A*H*(pow(1 - xnp2,R) - An*xnp2);
00623 c = 1.12*xn*A*(1 - xnp2);
00624 if( xp - xn == 1. )
00625 {
00626 c *= dsexp(0.006*powi(xn-1.,6)/charge);
00627 }
00628
00629 rate = COLL_CONST/(POW2(xn*charge)*phycon.sqrte);
00630 rate *= (A + y*c - y*D)*e1 + D;
00631
00632
00633
00634
00635
00636 HydColDwn_v = rate*gLo/gHi;
00637 }
00638 else
00639 {
00640 e1 = ee1(y);
00641 exp_int_2 = sexpy - y*e1;
00642
00643 A = 4.*powi(xn,4)/(1 - xnp2)*HydroOscilStr(xn,xp);
00644 D = A*H*(pow(1 - xnp2,R) - An*xnp2);
00645 c = 1.12*xn*A*(1 - xnp2);
00646 if( xp - xn == 1. )
00647 {
00648 c *= dsexp(0.006*powi(xn-1.,6)/charge);
00649 }
00650
00651 rate = COLL_CONST/(POW2(xn*charge)*phycon.sqrte);
00652 rate *= (A + y*c)*e1 + D*exp_int_2;
00653
00654
00655
00656
00657
00658 HydColDwn_v = rate*gLo/gHi/sexpy;
00659 }
00660
00661 HydColDwn_v = MAX2(HydColDwn_v,SMALLFLOAT);
00662
00663
00664 HydColDwn_v = HydColDwn_v*gHi*phycon.sqrte/COLL_CONST;
00665
00666 DEBUG_EXIT( "Hion_colldeexc_cs()" );
00667 return( HydColDwn_v );
00668 }