00001
00002
00003
00004
00005
00006
00007
00008 #include "cddefines.h"
00009 #include "iso.h"
00010 #include "dense.h"
00011 #include "helike.h"
00012 #include "helike_cs.h"
00013 #include "phycon.h"
00014 #include "physconst.h"
00015 #include "taulines.h"
00016 #include "hydro_vs_rates.h"
00017 #include "trace.h"
00018 #include "ionbal.h"
00019 #include "opacity.h"
00020 #include "heavy.h"
00021 #include "rfield.h"
00022 #include "atmdat.h"
00023 #include "thirdparty.h"
00024
00025
00026
00027 STATIC double S62_Therm_ave_coll_str( double EProjectile_eV );
00028
00029
00030
00031 STATIC double Therm_ave_coll_str_int_VF01( double EProjectileRyd);
00032 STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel );
00033 STATIC double L_mix_integrand_VF01( double alpha );
00034 STATIC double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi);
00035
00036 static long int global_n, global_l, global_l_prime, global_s, global_z, global_Collider;
00037
00038 static double global_bmax, global_red_vel, global_an, global_collider_charge;
00039
00040 static double global_I_energy_eV, global_deltaE, global_temp, global_osc_str, global_stat_weight;
00041
00042 static double kTRyd;
00043
00044
00045 static double ColliderMass[3] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0};
00046
00047
00048
00049 void HeCollid( long int nelem)
00050 {
00051
00052 double factor1 , ConvLTEPOP, *ColIonizPerN, DimaRate, crate;
00053 long int ipLo , ipHi, n;
00054
00055
00056
00057
00058 static double TeUsed[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00059 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0};
00060
00061
00062
00063
00064 static double TeUsedForCS[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00065 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0};
00066
00067 DEBUG_ENTRY( "HeCollid()" );
00068
00069
00070
00071 ColIonizPerN = (double*)MALLOC(sizeof(double)* (unsigned)(iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) );
00072
00073 for( n=0; n <= ( iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] ); ++n)
00074 {
00075 ColIonizPerN[n] = 0.;
00076 }
00077
00078
00079 if( TeUsed[nelem] != phycon.te )
00080 {
00081
00082 TeUsed[nelem] = phycon.te;
00083
00084 if( trace.lgTrace )
00085 {
00086 fprintf( ioQQQ,
00087 " HeCollid called nelem %li - will reeval Boltz fac, LTE dens\n",
00088 nelem );
00089 }
00090
00091
00092 factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
00093 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
00094
00095
00096 ConvLTEPOP = pow(factor1,1.5)/(2.*iso.stat_ion[ipHE_LIKE])/phycon.te32;
00097
00098 iso.lgPopLTE_OK[ipHE_LIKE][nelem] = true;
00099
00100 for( ipLo=ipHe1s1S; ipLo < (iso.numLevels_max[ipHE_LIKE][nelem]); ipLo++ )
00101 {
00102
00103 iso.ConBoltz[ipHE_LIKE][nelem][ipLo] =
00104 dsexp(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]/phycon.te_ryd);
00105
00106
00107
00108
00109
00110
00111
00112
00113 if( iso.ConBoltz[ipHE_LIKE][nelem][ipLo] >= SMALLDOUBLE )
00114 {
00115
00116 iso.PopLTE[ipHE_LIKE][nelem][ipLo] =
00117 iso.stat[ipHE_LIKE][nelem][ipLo] / iso.ConBoltz[ipHE_LIKE][nelem][ipLo]* ConvLTEPOP;
00118 ASSERT( iso.PopLTE[ipHE_LIKE][nelem][ipLo] < BIGDOUBLE );
00119 }
00120 else
00121 {
00122 iso.PopLTE[ipHE_LIKE][nelem][ipLo] = 0.;
00123 }
00124
00125
00126 if( iso.PopLTE[ipHE_LIKE][nelem][ipLo] <= 0. )
00127 {
00128 iso.lgPopLTE_OK[ipHE_LIKE][nelem] = false;
00129 }
00130 }
00131
00132
00133
00134 if( (TeUsedForCS[nelem] == 0.) ||
00135 ( TeUsedForCS[nelem]/phycon.te > 1.15 ) ||
00136 ( TeUsedForCS[nelem]/phycon.te < 0.85 ) )
00137 {
00138
00139
00140
00141
00142
00143
00144
00145 for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00146 {
00147 for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
00148 {
00149
00150 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs =
00151 HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
00152
00153 helike.cs_proton[nelem][ipHi][ipLo] =
00154 HeCSInterp( nelem , ipHi , ipLo, ipPROTON );
00155
00156 helike.cs_heplus[nelem][ipHi][ipLo] =
00157 HeCSInterp( nelem , ipHi , ipLo, ipHE_PLUS );
00158
00159
00160 ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs >= 0. );
00161 ASSERT( helike.cs_proton[nelem][ipHi][ipLo] >= 0. );
00162 ASSERT( helike.cs_heplus[nelem][ipHi][ipLo] >= 0. );
00163 }
00164 }
00165
00166 TeUsedForCS[nelem] = phycon.te;
00167 }
00168
00169 for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00170 {
00171 for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
00172 {
00173 double reduced_mass_proton = dense.AtomicWeight[nelem]*ColliderMass[ipPROTON]/
00174 (dense.AtomicWeight[nelem]+ColliderMass[ipPROTON])*ATOMIC_MASS_UNIT;
00175
00176 double reduced_mass_heplus = dense.AtomicWeight[nelem]*ColliderMass[ipHE_PLUS]/
00177 (dense.AtomicWeight[nelem]+ColliderMass[ipHE_PLUS])*ATOMIC_MASS_UNIT;
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = (float)(
00193 (
00194
00195 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs+
00196
00197 helike.cs_proton[nelem][ipHi][ipLo]*
00198 (float)(dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr)*
00199 pow(ELECTRON_MASS/reduced_mass_proton, 1.5) +
00200
00201 helike.cs_heplus[nelem][ipHi][ipLo]*
00202 (float)(dense.xIonDense[ipHELIUM][1]/dense.EdenHCorr)*
00203 pow(ELECTRON_MASS/reduced_mass_heplus, 1.5)
00204 ) / phycon.sqrte*COLL_CONST/(double)iso.stat[ipHE_LIKE][nelem][ipHi] );
00205
00206
00207
00208
00209 if( N_(ipHi) > iso.n_HighestResolved_max[ipHE_LIKE][nelem] &&
00210 N_(ipLo) <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00211 {
00212 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL *=
00213 (float)( (2./3.)*(log((double)N_(ipHi))+2) );
00214 }
00215
00216
00217
00218
00219 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK > 0. )
00220 {
00221 iso.Boltzmann[ipHE_LIKE][nelem][ipHi][ipLo] =
00222 sexp( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te );
00223 }
00224 else
00225 {
00226
00227 iso.Boltzmann[ipHE_LIKE][nelem][ipHi][ipLo] = 1.;
00228 }
00229 }
00230
00231
00232
00234 ASSERT( ipHi > 0 && ipHi < iso.numLevels_max[ipHE_LIKE][nelem]);
00235 if( nelem==ipHELIUM && helike.lgSetBenjamin )
00236 {
00237
00238 double BenCollIonParamA[31] = {
00239 3.36E-09,9.36E-08,1.33E-07,1.58E-07,1.58E-07,1.58E-07,1.81E-07,5.36E-07,
00240 6.58E-07,7.23E-07,7.81E-07,7.81E-07,7.93E-07,1.62E-06,1.87E-06,2.00E-06,
00241 2.11E-06,2.11E-06,2.11E-06,2.11E-06,2.13E-06,3.61E-06,4.02E-06,4.21E-06,
00242 4.39E-06,4.40E-06,4.40E-06,4.40E-06,4.40E-06,4.40E-06,4.44E-06};
00243
00244 double BenCollIonParamB[31] = {
00245 0.499,0.280,0.253,0.240,0.240,0.240,0.229,0.138,0.120,0.111,0.104,0.104,
00246 0.103,0.035,0.021,0.014,0.008,0.008,0.008,0.008,0.007,-0.047,-0.059,-0.064,
00247 -0.068,-0.068,-0.068,-0.068,-0.068,-0.068,-0.069};
00248
00249 double BenCollIonParamC[31] = {
00250 28.625,5.655,4.733,4.331,4.331,4.331,4.037,2.304,2.072,1.972,1.895,
00251 1.895,1.880,1.298,1.207,1.167,1.135,1.135,1.134,1.134,1.128,0.866,
00252 0.821,0.802,0.785,0.785,0.784,0.784,0.784,0.784,0.781};
00253
00254 ASSERT( ipHi < 31 );
00255
00256 iso.ColIoniz[ipHE_LIKE][nelem][ipHi] = BenCollIonParamA[ipHi] *
00257 pow( (double)phycon.te/10000., BenCollIonParamB[ipHi] ) *
00258 sexp( 10000.*BenCollIonParamC[ipHi]/phycon.te );
00259 }
00260 else
00261 {
00262 if( nelem == ipHELIUM )
00263 {
00264
00265
00266 iso.ColIoniz[ipHE_LIKE][nelem][ipHi] = hydro_vs_ioniz( ipHE_LIKE, nelem, ipHi );
00267 }
00268 else
00269 {
00270
00271
00272
00273
00274
00275 iso.ColIoniz[ipHE_LIKE][nelem][ipHi] = Hion_coll_ioniz_ratecoef(ipHE_LIKE , nelem , ipHi);
00276 }
00277
00278 if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n > iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00279 {
00280
00281
00282
00283 iso.ColIoniz[ipHE_LIKE][nelem][ipHi] *= 2.*N_(ipHi);
00284 }
00285
00286
00287 iso.ColIoniz[ipHE_LIKE][nelem][ipHi] *= iso.lgColl_ionize[ipHE_LIKE];
00288 }
00289
00290
00291 ColIonizPerN[iso.quant_desig[ipHE_LIKE][nelem][ipHi].n] += iso.ColIoniz[ipHE_LIKE][nelem][ipHi];
00292 }
00293
00294
00295
00296 helike.lgErrGenDone = true;
00297
00298
00299
00300
00301
00302
00303
00304
00305 if( !iso.lgLevelsLowered[ipHE_LIKE][nelem] && (nelem != ipHELIUM || !helike.lgSetBenjamin) )
00306 iso.ColIoniz[ipHE_LIKE][nelem][iso.numLevels_max[ipHE_LIKE][nelem]-1] *= 1.E5;
00307
00308
00309 if( iso.lgColl_ionize[ipHE_LIKE] && (nelem!=ipHELIUM || !helike.lgSetBenjamin ) )
00310 {
00311 for( n = 2; n < iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]; n++ )
00312 {
00313 ASSERT( (ColIonizPerN[n] < ColIonizPerN[n+1]) || (ColIonizPerN[n]==0.) );
00314 }
00315 }
00316
00317
00318
00319 if( opac.lgCaseB_HummerStorey )
00320 {
00321
00322
00323
00324 ipLo = ipHe1s1S;
00325 iso.ColIoniz[ipHE_LIKE][nelem][ipLo] = 0.;
00326 for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00327 {
00328 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 0.;
00329 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = 0.;
00330 }
00331 ipLo = ipHe2s1S;
00332 iso.ColIoniz[ipHE_LIKE][nelem][ipLo] = 0.;
00333 for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00334 {
00335 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 0.;
00336 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = 0.;
00337 }
00338 ipLo = ipHe2p1P;
00339 iso.ColIoniz[ipHE_LIKE][nelem][ipLo] = 0.;
00340 for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00341 {
00342 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs = 0.;
00343 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ColUL = 0.;
00344 }
00345 }
00346 }
00347 else
00348 {
00349 if( trace.lgTrace )
00350 {
00351 fprintf( ioQQQ,
00352 " HeCollid called nelem %li - no reeval Boltz fac, LTE dens\n",
00353 nelem );
00354 }
00355 }
00356
00357
00358 DimaRate = t_ADfA::Inst().coll_ion( nelem+1, 2 , phycon.te );
00359
00360 crate = DimaRate*dense.EdenHCorr*iso.lgColl_ionize[ipHE_LIKE];
00361
00362
00363 iso.ColIoniz[ipHE_LIKE][nelem][ipHe1s1S] = DimaRate*iso.lgColl_ionize[ipHE_LIKE];
00364
00365
00366
00367
00368 ionbal.CollIonRate_Ground[nelem][nelem-1][0] = crate;
00369
00370
00371 ionbal.CollIonRate_Ground[nelem][nelem-1][1] = (crate*
00372 rfield.anu[Heavy.ipHeavy[nelem][nelem-1]-1]* EN1RYD);
00373
00374
00375 if( nelem == ipHELIUM && !helike.lgSetBenjamin )
00376 {
00377 double c2sion , c2pion;
00378
00379 c2sion = 8.8e-10*(double)phycon.sqrte*iso.ConBoltz[ipHE_LIKE][nelem][ipHe2s3S];
00380 c2pion = 3.*8.8e-10*(double)phycon.sqrte*iso.ConBoltz[ipHE_LIKE][nelem][ipHe2p3P0];
00381 iso.ColIoniz[ipHE_LIKE][nelem][ipHe2s3S] = c2sion*iso.lgColl_ionize[ipHE_LIKE];
00382
00383
00384
00385 iso.ColIoniz[ipHE_LIKE][nelem][ipHe2p3P0] = c2pion*iso.lgColl_ionize[ipHE_LIKE]/3.;
00386 iso.ColIoniz[ipHE_LIKE][nelem][ipHe2p3P1] = c2pion*iso.lgColl_ionize[ipHE_LIKE]/3.;
00387 iso.ColIoniz[ipHE_LIKE][nelem][ipHe2p3P2] = c2pion*iso.lgColl_ionize[ipHE_LIKE]/3.;
00388 }
00389 {
00390
00391
00392 enum {PRINTIT=false};
00393
00394 if( PRINTIT )
00395 {
00396
00397 nelem = ipHELIUM;
00398 for( ipHi=1; ipHi<iso.numLevels_max[ipHE_LIKE][ipHELIUM]; ++ipHi )
00399 {
00400 for( ipLo=0; ipLo<ipHi; ++ipLo )
00401 {
00402 fprintf(ioQQQ,"Z %li Lo %li n %li s %li l %li \t ",
00403 nelem+1, ipLo,
00404 N_(ipLo),
00405 S_(ipLo),
00406 L_(ipLo) );
00407 fprintf(ioQQQ," Hi %li n %li s %li l %li \t",
00408 ipHi,
00409 N_(ipHi),
00410 S_(ipHi),
00411 L_(ipHi) );
00412 fprintf(ioQQQ,"%.4e\t%.4e\tcs\t%.4e\n",
00413 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
00414 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul ,
00415 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs);
00416 }
00417 }
00418 exit(1);
00419 }
00420 }
00421
00422
00423 free( ColIonizPerN );
00424
00425 DEBUG_EXIT( "HeCollid()" );
00426 }
00427
00428
00429
00430 float HeCSInterp(long int nelem,
00431 long int ipHi,
00432 long int ipLo,
00433 long int Collider )
00434 {
00435 float cs, factor1;
00436
00437
00438
00439 const char *where = " ";
00440
00441 if( iso.lgColl_excite[ipHE_LIKE] == false )
00442 return (float)1E-10;
00443
00444 if( nelem == ipHELIUM )
00445 {
00446
00447 cs = AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
00448 }
00449 else
00450 {
00451
00452 cs = IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
00453 }
00454
00455
00456 putError(nelem, ipHi, ipLo, IPCOLLIS, 0.15f );
00457
00458 ASSERT( cs >= 0.f );
00459
00460
00461
00462
00463 ASSERT( factor1 >= 0.f || nelem!=ipHELIUM );
00464 if( factor1 < 0.f )
00465 {
00466 ASSERT( helike.lgCS_Vriens );
00467
00468 factor1 = 1.f;
00469 }
00470
00471
00472
00473 cs *= factor1;
00474
00475 {
00476
00477 enum {DEBUG_LOC=false};
00478
00479
00480 if( DEBUG_LOC && ( nelem==ipOXYGEN ) && (cs > 0.f) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == 2)
00481 && ( iso.quant_desig[ipHE_LIKE][nelem][ipLo].n <= 2 ) )
00482 fprintf(ioQQQ,"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n",
00483 iso.quant_desig[ipHE_LIKE][nelem][ipLo].n , iso.quant_desig[ipHE_LIKE][nelem][ipLo].s ,
00484 iso.quant_desig[ipHE_LIKE][nelem][ipLo].l , iso.quant_desig[ipHE_LIKE][nelem][ipHi].n ,
00485 iso.quant_desig[ipHE_LIKE][nelem][ipHi].s , iso.quant_desig[ipHE_LIKE][nelem][ipHi].l , cs,where);
00486 }
00487
00488 return MAX2(cs, 1.e-10f);
00489 }
00490
00491 float AtomCSInterp(long int nelem,
00492 long int ipHi,
00493 long int ipLo,
00494 float *factor1,
00495 const char **where,
00496 long int Collider )
00497 {
00498 long ipArray;
00499 float cs;
00500
00501 DEBUG_ENTRY( "AtomCSInterp()" );
00502
00503 ASSERT( nelem == ipHELIUM );
00504
00505
00506 cs = -1.f;
00507
00508
00509
00510
00511 *factor1 = -1.f;
00512
00513
00514
00515
00516
00517
00518 if( ipLo >= ipHe2p3P0 && ipHi <= ipHe2p3P2 && Collider==ipELECTRON )
00519 {
00520 *factor1 = 1.f;
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533 if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P1 )
00534 {
00535 cs = 1.2f;
00536 }
00537 else if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P2 )
00538 {
00539 cs = 2.1f;
00540 }
00541 else if( ipLo == ipHe2p3P1 && ipHi == ipHe2p3P2 )
00542 {
00543 cs = 6.0f;
00544 }
00545 else
00546 {
00547 cs = 1.0f;
00548 TotalInsanity();
00549 }
00550
00551 *where = "Berr ";
00552
00553 }
00554
00555
00556
00557 else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n <= 5 &&
00558 ( ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) &&
00559 helike.HeCS[nelem][ipHi][ipLo][0] >= 0.f && Collider== ipELECTRON )
00560 {
00561 ASSERT( *factor1 == -1.f );
00562 ASSERT( ipLo < ipHi );
00563 ASSERT( ipHe2p3P0 == 3 );
00564
00565
00566 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
00567 {
00568
00569
00570
00571 *factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
00572
00573
00574 ASSERT( ipHi > ipHe2p3P2 );
00575 }
00576
00577 else if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00578 {
00579 ASSERT( ipLo < ipHe2p3P0 );
00580
00581 *factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
00582 }
00583
00584 else
00585 {
00586 *factor1 = 1.f;
00587 }
00588
00589
00590
00591
00592
00593
00594
00595
00596 if( phycon.alogte <= helike.CSTemp[0] )
00597 {
00598 cs = helike.HeCS[nelem][ipHi][ipLo][0];
00599 }
00600 else if( phycon.alogte >= helike.CSTemp[helike.nCS-1] )
00601 {
00602 cs = helike.HeCS[nelem][ipHi][ipLo][helike.nCS-1];
00603 }
00604 else
00605 {
00606 float flow;
00607
00608
00609 ipArray = (long)((phycon.alogte - helike.CSTemp[0])/(helike.CSTemp[1]-helike.CSTemp[0]));
00610 ASSERT( ipArray < helike.nCS );
00611 ASSERT( ipArray >= 0 );
00612
00613 flow = (float)( (phycon.alogte - helike.CSTemp[ipArray])/
00614 (helike.CSTemp[ipArray+1]-helike.CSTemp[ipArray]));
00615 ASSERT( (flow >= 0.f) && (flow <= 1.f) );
00616
00617 cs = helike.HeCS[nelem][ipHi][ipLo][ipArray] * (1.f-flow) +
00618 helike.HeCS[nelem][ipHi][ipLo][ipArray+1] * flow;
00619 }
00620
00621 *where = "Bray ";
00622
00623
00624 if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
00625
00626 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
00627 else
00628 {
00629
00630 cs *= (float)iso.lgColl_excite[ipHE_LIKE];
00631
00632 if( ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >= 5 ) && helike.lgSetBenjamin )
00633 cs = 0.f;
00634 }
00635
00636 ASSERT( cs >= 0.f );
00637
00638 }
00639
00640 else if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) &&
00641 (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == iso.quant_desig[ipHE_LIKE][nelem][ipLo].s ) )
00642 {
00643 ASSERT( *factor1 == -1.f );
00644 *factor1 = 1.f;
00645
00646
00647 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] );
00648
00649 if( (iso.quant_desig[ipHE_LIKE][nelem][ipLo].l <=2) &&
00650 abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1 )
00651 {
00652
00653
00654
00655 cs = (float)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
00656 }
00657 else if( helike.lgCS_Vrinceanu == true )
00658 {
00659 if( iso.quant_desig[ipHE_LIKE][nelem][ipLo].l >=3 &&
00660 iso.quant_desig[ipHE_LIKE][nelem][ipHi].l >=3 )
00661 {
00662
00663
00664
00665 cs = (float)CS_l_mixing_VF01( nelem,
00666 iso.quant_desig[ipHE_LIKE][nelem][ipLo].n,
00667 iso.quant_desig[ipHE_LIKE][nelem][ipLo].l,
00668 iso.quant_desig[ipHE_LIKE][nelem][ipHi].l,
00669 iso.quant_desig[ipHE_LIKE][nelem][ipLo].s,
00670 (double)phycon.te,
00671 Collider );
00672 }
00673 else
00674 {
00675 cs = 0.f;
00676 }
00677 }
00678
00679 else if( abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1)
00680 {
00681
00682
00683 cs = (float)CS_l_mixing_PS64( nelem, ipLo, ipHi, Collider);
00684 }
00685 else
00686 {
00687
00688 cs = 0.f;
00689 }
00690
00691
00692 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
00693 {
00694 *factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
00695 }
00696
00697
00698 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00699 {
00700 *factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
00701 }
00702
00703 *where = "lmix ";
00704 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
00705 }
00706
00707
00708 else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n != iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
00709 {
00710 ASSERT( *factor1 == -1.f );
00711
00712
00713
00714
00715 if( helike.lgCS_Vriens )
00716 {
00717
00718
00719 cs = (float)CS_VS80(
00720 ipHE_LIKE,
00721 nelem,
00722 ipHi,
00723 ipLo,
00724 phycon.te,
00725 Collider );
00726
00727 *factor1 = 1.f;
00728 *where = "Vriens";
00729 }
00730 else if( helike.lgCS_None )
00731 {
00732 cs = 0.f;
00733 *factor1 = 1.f;
00734 *where = "no gb";
00735 }
00736 else if( helike.nCS_new )
00737 {
00738 *factor1 = 1.f;
00739
00740
00741
00742
00743
00744 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul > 1. )
00745 {
00746
00747 double x =
00748 log10(MAX2(34.7,EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN));
00749
00750 if( helike.nCS_new == 1 )
00751 {
00752
00753
00754 if( x < 4.5 )
00755 {
00756
00757 cs = (float)pow( 10. , -1.45*x + 6.75);
00758 }
00759 else
00760 {
00761
00762 cs = (float)pow( 10. , -3.33*x+15.15);
00763 }
00764 }
00765 else if( helike.nCS_new == 2 )
00766 {
00767
00768 cs = (float)pow( 10. , -2.3*x+10.3);
00769 }
00770 else
00771 TotalInsanity();
00772 }
00773 else
00774 {
00775
00776 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN < 25119.f )
00777 {
00778 cs = 0.631f;
00779 }
00780 else
00781 {
00782 cs = (float)pow(10.,
00783 -3.*log10(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN)+12.8);
00784 }
00785 }
00786
00787 *where = "newgb";
00788 }
00789 else
00790 TotalInsanity();
00791
00792
00793 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
00794 {
00795 *factor1 = (2.f*((float)ipLo-3.f)+1.f) / 9.f;
00796 }
00797
00798
00799 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00800 {
00801 *factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
00802 }
00803
00804
00805 cs *= (float)iso.lgColl_excite[ipHE_LIKE];
00806
00807
00808 if( ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >= 5 ) && helike.lgSetBenjamin )
00809 cs = 0.f;
00810 }
00811 else
00812 {
00813
00814
00815
00816 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].s != iso.quant_desig[ipHE_LIKE][nelem][ipLo].s );
00817 cs = 0.f;
00818 *factor1 = 1.f;
00819 }
00820
00821 ASSERT( cs >= 0.f );
00822
00823
00824 if( helike.lgSetBenjamin )
00825 {
00826 if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >=5 &&
00827 iso.quant_desig[ipHE_LIKE][nelem][ipLo].n != iso.quant_desig[ipHE_LIKE][nelem][ipHi].n )
00828 {
00829 ASSERT( cs == 0.f );
00830 }
00831 else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n >=5 &&
00832 iso.quant_desig[ipHE_LIKE][nelem][ipLo].s != iso.quant_desig[ipHE_LIKE][nelem][ipHi].s )
00833 {
00834 cs = 0.f;
00835 }
00836 }
00837
00838
00839
00840 DEBUG_EXIT( "AtomCSInterp()" );
00841
00842 return(cs);
00843
00844 }
00845
00846
00847 float IonCSInterp( long nelem , long ipHi , long ipLo, float *factor1, const char **where, long Collider )
00848 {
00849 float cs;
00850
00851 DEBUG_ENTRY( "IonCSInterp()" );
00852
00853 ASSERT( nelem > ipHELIUM );
00854 ASSERT( nelem < LIMELM );
00855
00856
00857 cs = -1.f;
00858
00859
00860
00861
00862 *factor1 = -1.f;
00863
00864
00865
00866
00867
00868
00869
00870
00871 if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n==2
00872 && iso.quant_desig[ipHE_LIKE][nelem][ipLo].n<=2 && Collider==ipELECTRON)
00873 {
00874 *where = "Zhang";
00875 *factor1 = 1.;
00876
00877
00878 if( ipLo == ipHe1s1S )
00879 {
00880 switch( ipHi )
00881 {
00882 case 1:
00883 cs = 0.25f/(float)POW2(nelem+1.);
00884 break;
00885 case 2:
00886 cs = 0.4f/(float)POW2(nelem+1.);
00887 break;
00888 case 3:
00889 cs = 0.15f/(float)POW2(nelem+1.);
00890 break;
00891 case 4:
00892 cs = 0.45f/(float)POW2(nelem+1.);
00893 break;
00894 case 5:
00895 cs = 0.75f/(float)POW2(nelem+1.);
00896 break;
00897 case 6:
00898 cs = 1.3f/(float)POW2(nelem+1.);
00899 break;
00900 default:
00901 TotalInsanity();
00902 break;
00903 }
00904 cs *= (float)iso.lgColl_excite[ipHE_LIKE];
00905 }
00906
00907 else if( ipLo == ipHe2s3S )
00908 {
00909 switch( ipHi )
00910 {
00911 case 2:
00912 cs = 2.75f/(float)POW2(nelem+1.);
00913 break;
00914 case 3:
00915 cs = 60.f/(float)POW2(nelem+1.);
00916 break;
00917 case 4:
00918 cs = 180.f/(float)POW2(nelem+1.);
00919 break;
00920 case 5:
00921 cs = 300.f/(float)POW2(nelem+1.);
00922 break;
00923 case 6:
00924 cs = 5.8f/(float)POW2(nelem+1.);
00925 break;
00926 default:
00927 TotalInsanity();
00928 break;
00929 }
00930 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
00931 }
00932
00933 else if( ipLo == ipHe2s1S )
00934 {
00935 switch( ipHi )
00936 {
00937 case 3:
00938 cs = 0.56f/(float)POW2(nelem+1.);
00939 break;
00940 case 4:
00941 cs = 1.74f/(float)POW2(nelem+1.);
00942 break;
00943 case 5:
00944 cs = 2.81f/(float)POW2(nelem+1.);
00945 break;
00946 case 6:
00947 cs = 190.f/(float)POW2(nelem+1.);
00948 break;
00949 default:
00950 TotalInsanity();
00951 break;
00952 }
00953 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
00954 }
00955
00956 else if( ipLo == ipHe2p3P0 )
00957 {
00958 switch( ipHi )
00959 {
00960 case 4:
00961 cs = 8.1f/(float)POW2(nelem+1.);
00962 break;
00963 case 5:
00964 cs = 8.2f/(float)POW2(nelem+1.);
00965 break;
00966 case 6:
00967 cs = 3.9f/(float)POW2(nelem+1.);
00968 break;
00969 default:
00970 TotalInsanity();
00971 break;
00972 }
00973 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
00974 }
00975
00976 else if( ipLo == ipHe2p3P1 )
00977 {
00978 switch( ipHi )
00979 {
00980 case 5:
00981 cs = 30.f/(float)POW2(nelem+1.);
00982 break;
00983 case 6:
00984 cs = 11.7f/(float)POW2(nelem+1.);
00985 break;
00986 default:
00987 TotalInsanity();
00988 break;
00989 }
00990 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
00991 }
00992
00993 else if( ipLo == ipHe2p3P2 )
00994 {
00995
00996 cs = 19.5f/(float)POW2(nelem+1.) * (float)iso.lgColl_l_mixing[ipHE_LIKE];
00997 }
00998 else
00999 TotalInsanity();
01000
01001
01002 }
01003
01004
01005 else if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) &&
01006 (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == iso.quant_desig[ipHE_LIKE][nelem][ipLo].s ) )
01007 {
01008 ASSERT( *factor1 == -1.f );
01009 *factor1 = 1.f;
01010
01011 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] );
01012
01013 if( (iso.quant_desig[ipHE_LIKE][nelem][ipLo].l <=2) &&
01014 abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1 )
01015 {
01016
01017
01018
01019 cs = (float)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
01020 }
01021 else if( helike.lgCS_Vrinceanu == true )
01022 {
01023 if( iso.quant_desig[ipHE_LIKE][nelem][ipLo].l >=3 &&
01024 iso.quant_desig[ipHE_LIKE][nelem][ipHi].l >=3 )
01025 {
01026
01027
01028
01029 cs = (float)CS_l_mixing_VF01( nelem,
01030 iso.quant_desig[ipHE_LIKE][nelem][ipLo].n,
01031 iso.quant_desig[ipHE_LIKE][nelem][ipLo].l,
01032 iso.quant_desig[ipHE_LIKE][nelem][ipHi].l,
01033 iso.quant_desig[ipHE_LIKE][nelem][ipLo].s,
01034 (double)phycon.te,
01035 Collider );
01036 }
01037 else
01038 {
01039 cs = 0.f;
01040 }
01041 }
01042
01043 else if( abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)== 1)
01044 {
01045
01046
01047 cs = (float)CS_l_mixing_PS64( nelem, ipLo, ipHi, Collider);
01048 }
01049 else
01050 {
01051
01052 cs = 0.f;
01053 }
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
01066 {
01067 *factor1 = (2.f*((float)ipHi-3.f)+1.f) / 9.f;
01068 }
01069
01070 *where = "lmix ";
01071 cs *= (float)iso.lgColl_l_mixing[ipHE_LIKE];
01072 }
01073
01074
01075 else if( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n != iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
01076 {
01077 if( helike.lgCS_Vriens )
01078 {
01079
01080
01081
01082 cs = (float)CS_VS80(
01083 ipHE_LIKE,
01084 nelem,
01085 ipHi,
01086 ipLo,
01087 phycon.te,
01088 Collider );
01089
01090 *factor1 = 1.f;
01091 *where = "Vriens";
01092 }
01093 else
01094 {
01095
01096
01097 long int nlo = iso.quant_desig[ipHE_LIKE][nelem][ipLo].n;
01098
01099 cs = (float)CS_VS80(
01100 ipHE_LIKE,
01101 nelem,
01102 ipHi,
01103 ipLo,
01104 phycon.te,
01105 Collider );
01106
01107
01108
01109 if( nlo == 1)
01110 nlo = 0;
01111
01112
01113
01114
01115 # if 0
01116
01117 cs = (float)Hion_colldeexc_cs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].n, nlo, nelem, ipHE_LIKE );
01118 cs *= iso.lgColl_excite[ipHE_LIKE];
01119 if( nelem == ipCARBON )
01120 # endif
01121
01122 *where = "hydro";
01123
01124 }
01125 }
01126 else
01127 {
01128
01129
01130
01131 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == iso.quant_desig[ipHE_LIKE][nelem][ipLo].n );
01132 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].s != iso.quant_desig[ipHE_LIKE][nelem][ipLo].s );
01133 cs = 0.f;
01134 *where = "spin ";
01135 }
01136
01137 ASSERT( cs >= 0.f );
01138
01139 DEBUG_EXIT( "IonCSInterp()" );
01140
01141
01142
01143 return(cs);
01144 }
01145
01146
01147
01148
01149 double CS_l_mixing_S62(
01150 long ipISO,
01151 long nelem ,
01152 long ipLo ,
01153 long ipHi ,
01154 double temp,
01155 long Collider)
01156 {
01157
01158 double coll_str, deltaE;
01159
01160 DEBUG_ENTRY( "CS_l_mixing_S62()" );
01161
01162 global_temp = temp;
01163 global_stat_weight = iso.stat[ipISO][nelem][ipLo];
01164
01165
01166
01167 global_deltaE = EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg/EN1EV;
01168 deltaE = global_deltaE;
01169 global_I_energy_eV = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][0];
01170 global_Collider = Collider;
01171
01172 ASSERT( TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD) > 0. );
01173
01174 global_osc_str = EmisLines[ipISO][nelem][ipHi][ipLo].Aul/
01175 (TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD));
01176
01177
01178 coll_str = qg32( 0.0, 1.0, S62_Therm_ave_coll_str);
01179 coll_str += qg32( 1.0, 10.0, S62_Therm_ave_coll_str);
01180 ASSERT( coll_str > 0. );
01181
01182 DEBUG_EXIT( "CS_l_mixing_S62()" );
01183
01184
01185
01186 return coll_str;
01187 }
01188
01189
01190 STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT )
01191 {
01192
01193 double integrand, cross_section, osc_strength, coll_str, zOverB2;
01194 double deltaE, betaone, zeta_of_betaone, cs2, temp, stat_weight;
01195 double I_energy_eV, Dubya, proj_energy;
01196 long i, Collider;
01197
01198 DEBUG_ENTRY( "S62_Therm_ave_coll_str()" );
01199
01200
01201 proj_energy = proj_energy_overKT * phycon.te / EVDEGK;
01202
01203 deltaE = global_deltaE;
01204 osc_strength = global_osc_str;
01205 temp = (double)global_temp;
01206 stat_weight = global_stat_weight;
01207 I_energy_eV = global_I_energy_eV;
01208 Collider = global_Collider;
01209
01210
01211
01212 proj_energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
01213
01214
01215 proj_energy += deltaE;
01216 Dubya = 0.5*(2.*proj_energy-deltaE);
01217 ASSERT( Dubya > 0. );
01218
01219
01220
01221 ASSERT( I_energy_eV > 0. );
01222 ASSERT( osc_strength > 0. );
01223
01224
01225 zOverB2 = 0.5*POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
01226
01227 ASSERT( zOverB2 > 0. );
01228
01229 if( zOverB2 > 100. )
01230 {
01231 betaone = sqrt( 1./zOverB2 );
01232 }
01233 else if( zOverB2 < 0.54 )
01234 {
01235
01236 betaone = (1./3.)*(log(PI)-log(zOverB2)+1.28);
01237 if( betaone > 2.38 )
01238 {
01239
01240 betaone += 0.5*(log(PI)-log(zOverB2));
01241 betaone *= 0.5;
01242 }
01243 }
01244 else
01245 {
01246 long ip_zOverB2 = 0;
01247 double zetaOVERbeta2[101] = {
01248 1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
01249 7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
01250 4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
01251 3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
01252 2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
01253 1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
01254 1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
01255 7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
01256 4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
01257 3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
01258 2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
01259 1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
01260 7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
01261
01262 ASSERT( zOverB2 >= zetaOVERbeta2[100] );
01263
01264
01265 for( i=0; i< 100; ++i )
01266 {
01267 if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
01268 {
01269
01270 ip_zOverB2 = i;
01271 break;
01272 }
01273 }
01274
01275 ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
01276
01277 betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) /
01278 (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
01279 (pow(10., ((double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((double)ip_zOverB2)/100. - 1.)) +
01280 pow(10., (double)ip_zOverB2/100. - 1.);
01281
01282 }
01283
01284 zeta_of_betaone = zOverB2 * POW2(betaone);
01285
01286
01287 cs2 = 0.5*zeta_of_betaone + betaone * bessel_k0(betaone) * bessel_k1(betaone);
01288
01289
01290 cross_section = cs2;
01291
01292
01293 cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
01294
01295
01296 coll_str = cross_section * stat_weight * POW2( proj_energy/EVRYD );
01297
01298 integrand = exp( -1.*(proj_energy-deltaE)*EVDEGK/temp ) * coll_str;
01299
01300 DEBUG_EXIT( "S62_Therm_ave_coll_str()" );
01301 return integrand;
01302 }
01303
01304
01305 double CS_l_mixing_PS64(
01306 long nelem ,
01307 long ipLo ,
01308 long ipHi ,
01309 long Collider)
01310 {
01311
01312 double cs, Dul,
01313 TwoLogDebye, TwoLogRc1,
01314 factor1, factor2,
01315
01316 bestfactor, factorpart,
01317 reduced_mass, reduced_mass_2_emass,
01318 rate, tau ;
01319
01320 const double ChargIncoming = 1.;
01321
01322 DEBUG_ENTRY( "CS_l_mixing_PS64()" );
01323
01324 ASSERT( ipHi > ipLo );
01325
01326
01327
01328
01329 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
01330 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
01331
01332
01333 reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;
01334
01335
01336 tau = helike.Lifetime[nelem][ipLo];
01337
01338
01339
01340
01341
01342
01343 TwoLogDebye = 1.68 + log10( phycon.te / MIN2(1e12 , dense.eden ) );
01344
01345
01346 TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass );
01347
01348
01349
01350 Dul = POW2( ChargIncoming / (double)nelem ) * 6. * POW2( (double)iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) *
01351 ( POW2((double)iso.quant_desig[ipHE_LIKE][nelem][ipLo].n) -
01352 POW2((double)iso.quant_desig[ipHE_LIKE][nelem][ipLo].l) - iso.quant_desig[ipHE_LIKE][nelem][ipLo].l - 1);
01353
01354 ASSERT( Dul > 0. );
01355 ASSERT( phycon.te / Dul / reduced_mass_2_emass > 0. );
01356
01357 factorpart = (11.54 + log10( phycon.te / Dul / reduced_mass_2_emass ) );
01358
01359 if( (factor1 = factorpart + TwoLogDebye) <= 0.)
01360 factor1 = BIGDOUBLE;
01361 if( (factor2 = factorpart + TwoLogRc1) <= 0.)
01362 factor2 = BIGDOUBLE;
01363
01364
01365 bestfactor = MIN2(factor1,factor2);
01366
01367
01368 ASSERT( (bestfactor>0.) && (bestfactor<100.) );
01369
01370
01371 rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dul / phycon.sqrte * bestfactor;
01372
01373
01374 rate /= 2.;
01375
01376
01377 cs = rate / COLL_CONST * phycon.sqrte * iso.stat[ipHE_LIKE][nelem][ipLo];
01378
01379 ASSERT( cs > 0. );
01380
01381 DEBUG_EXIT( "CS_l_mixing_PS64()" );
01382
01383
01384
01385 return cs;
01386 }
01387
01388
01389
01390
01391
01392
01393 double CS_l_mixing_VF01(long int nelem,
01394 long int n,
01395 long int l,
01396 long int lp,
01397 long int s,
01398 double temp,
01399 long int Collider )
01400 {
01401
01402 double coll_str;
01403
01404 DEBUG_ENTRY( "CS_l_mixing_VF01()" );
01405
01406 global_z = nelem;
01407 global_n = n;
01408 global_l = l;
01409 global_l_prime = lp;
01410 global_s = s;
01411 global_temp = temp;
01412 global_collider_charge = 1.;
01413 global_Collider = Collider;
01414
01415 ASSERT( l != 0 );
01416 ASSERT( lp != 0 );
01417
01418 kTRyd = temp / TE1RYD;
01419 if( helike.lgCS_therm_ave == false )
01420 {
01421
01422
01423
01424
01425 if( (dense.eden > 10000.) && (dense.eden < 1E10 ) )
01426 {
01427 coll_str = qg32( 0.0, 6.0, Therm_ave_coll_str_int_VF01);
01428 }
01429 else
01430 {
01431
01432 coll_str = collision_strength_VF01( kTRyd, false );
01433 }
01434 }
01435 else
01436 {
01437
01438 coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_VF01);
01439 coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VF01);
01440 }
01441
01442 DEBUG_EXIT( "CS_l_mixing_VF01()" );
01443
01444 return coll_str;
01445 }
01446
01447
01448 STATIC double Therm_ave_coll_str_int_VF01( double EOverKT )
01449 {
01450 double integrand;
01451
01452 DEBUG_ENTRY( "Therm_ave_coll_str_int_VF01()" );
01453
01454 integrand = exp( -1.*EOverKT ) * collision_strength_VF01( EOverKT * kTRyd, false );
01455
01456 DEBUG_EXIT( "Therm_ave_coll_str_int_VF01()" );
01457
01458 return integrand;
01459 }
01460
01461 STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel )
01462 {
01463
01464 double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
01465 double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
01466 double ColliderCharge = global_collider_charge;
01467 double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
01468 double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 ;
01469
01470 long nelem, n, l, lp, s, ipLo, ipHi, Collider = global_Collider;
01471
01472 DEBUG_ENTRY( "collision_strength_VF01()" );
01473
01474 nelem = global_z;
01475 n = global_n;
01476 l = global_l;
01477 lp = global_l_prime;
01478 s = global_s;
01479
01480
01481 ipLo = QuantumNumbers2Index[nelem][n][l][s];
01482 ipHi = QuantumNumbers2Index[nelem][n][lp][s];
01483
01484 stat_weight = iso.stat[ipHE_LIKE][nelem][ipLo];
01485
01486
01487 ASSERT( n > 0 );
01488 ASSERT( l > 0 );
01489 ASSERT( lp > 0 );
01490
01491
01492 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
01493 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
01494 ASSERT( reduced_mass > 0. );
01495
01496
01497
01498 aveRadius = (BOHR_RADIUS_CM/nelem)*POW2((double)n);
01499 ASSERT( aveRadius < 1.e-4 );
01500
01501
01502
01503 ASSERT( aveRadius > 3.9/LIMELM * BOHR_RADIUS_CM );
01504 global_an = aveRadius;
01505
01506 RMSv = nelem*POW2(ELEM_CHARGE_ESU)/n/H_BAR;
01507 ASSERT( RMSv > 0. );
01508
01509 ASSERT( ColliderMass[Collider] > 0. );
01510
01511 if( lgParamIsRedVel == true )
01512 {
01513
01514 reduced_vel = velOrEner;
01515
01516
01517 E_Proj_Ryd = 0.5 * POW2( reduced_vel * RMSv ) * ColliderMass[Collider] *
01518 PROTON_MASS / EN1RYD;
01519 }
01520 else
01521 {
01522
01523 E_Proj_Ryd = velOrEner;
01524 reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/ColliderMass[Collider]/PROTON_MASS )/RMSv;
01525 }
01526
01527
01528 ASSERT( reduced_vel > 1.e-10 );
01529 ASSERT( reduced_vel < 1.e10 );
01530
01531 global_red_vel = reduced_vel;
01532
01533
01534 ConstantFactors = 4.5*PI*POW2(ColliderCharge*aveRadius/reduced_vel);
01535
01536
01537 reduced_b_min = 1.5 * ColliderCharge / reduced_vel;
01538 ASSERT( reduced_b_min > 1.e-10 );
01539 ASSERT( reduced_b_min < 1.e10 );
01540
01541
01542
01543
01544
01545 quantum_defect1 = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]);
01546 quantum_defect2 = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi]);
01547
01548
01549 ASSERT( fabs(quantum_defect1) < 1.0 );
01550 ASSERT( fabs(quantum_defect1) > 0.0 );
01551 ASSERT( fabs(quantum_defect2) < 1.0 );
01552 ASSERT( fabs(quantum_defect2) > 0.0 );
01553
01554
01555 omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
01556
01557 omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)lp/(double)n)) / POW3( (double)n ) / (double)lp );
01558
01559 omega_qd = 0.5*( omega_qd1 + omega_qd2 );
01560
01561 ASSERT( omega_qd > 0. );
01562
01563 reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
01564 reduced_b_max = MAX2( reduced_b_max, reduced_b_min );
01565 ASSERT( reduced_b_max > 0. );
01566 global_bmax = reduced_b_max*aveRadius;
01567 alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max);
01568 alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min);
01569
01570 ASSERT( alphamin > 0. );
01571 ASSERT( alphamax > 0. );
01572
01573 alphamin = MAX2( alphamin, 1.e-30 );
01574 alphamax = MAX2( alphamax, 1.e-20 );
01575
01576 CSIntegral = 0.;
01577
01578 if( alphamax > alphamin )
01579 {
01580
01581 step = (alphamax - alphamin)/5.;
01582 alpha1 = alphamin;
01583 CSIntegral += qg32( alpha1, alpha1+step, L_mix_integrand_VF01);
01584 CSIntegral += qg32( alpha1+step, alpha1+4.*step, L_mix_integrand_VF01);
01585 }
01586
01587
01588 cross_section = ConstantFactors * CSIntegral;
01589
01590
01591 coll_str = cross_section * stat_weight / PI / BOHR_RADIUS_CM / BOHR_RADIUS_CM;
01592 coll_str *= 0.25 * POW2( E_Proj_Ryd );
01593
01594 coll_str = MAX2( SMALLFLOAT, coll_str);
01595
01596 DEBUG_EXIT( "collision_strength_VF01()" );
01597
01598 return coll_str;
01599 }
01600
01601 STATIC double L_mix_integrand_VF01( double alpha )
01602 {
01603 double integrand, deltaPhi, b;
01604
01605 DEBUG_ENTRY( "L_mix_integrand_VF01()" );
01606
01607 ASSERT( alpha >= 1.e-30 );
01608 ASSERT( global_bmax > 0. );
01609 ASSERT( global_red_vel > 0. );
01610
01611
01612 b = 1.5*global_collider_charge*global_an/(global_red_vel * alpha);
01613
01614 if( b < global_bmax )
01615 {
01616 deltaPhi = -1.*PI + 2.*asin(b/global_bmax);
01617 }
01618 else
01619 {
01620 deltaPhi = 0.;
01621 }
01622 integrand = 1./alpha/alpha/alpha;
01623 integrand *= StarkCollTransProb_VF01( global_n, global_l, global_l_prime, alpha, deltaPhi );
01624
01625 DEBUG_EXIT( "L_mix_integrand_VF01()" );
01626
01627 return integrand;
01628 }
01629
01630 STATIC double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi)
01631 {
01632 double probability;
01633 double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
01634
01635 DEBUG_ENTRY( "StarkCollTransProb_VF01()" );
01636
01637 ASSERT( alpha > 0. );
01638
01639
01640 cosU1 = 2.*POW2((double)l/(double)n) - 1.;
01641 cosU2 = 2.*POW2((double)lp/(double)n) - 1.;
01642
01643 sinU1 = sqrt( 1. - cosU1*cosU1 );
01644 sinU2 = sqrt( 1. - cosU2*cosU2 );
01645
01646
01647 cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
01648 sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
01649 cosChi = 2. * POW2( cosChiOver2 ) - 1.;
01650
01651 if( l == 0 )
01652 {
01653 if( -1.*cosU2 - cosChi < 0. )
01654 {
01655 probability = 0.;
01656 }
01657 else
01658 {
01659
01660 ASSERT( sinChiOver2 > 0. );
01661 ASSERT( sinChiOver2*sinChiOver2 > POW2((double)lp/(double)n) );
01662
01663
01664 probability = (double)lp/(POW2((double)n)*sinChiOver2*sqrt( POW2(sinChiOver2) - POW2((double)lp/(double)n) ) );
01665 }
01666 }
01667 else
01668 {
01669 double OneMinusCosChi = 1. - cosChi;
01670
01671 if( OneMinusCosChi == 0. )
01672 {
01673 double hold = sin( deltaPhi / 2. );
01674
01675 OneMinusCosChi = 8. * alpha * alpha * POW2( hold );
01676 }
01677
01678 if( OneMinusCosChi == 0. )
01679 {
01680 probability = 0.;
01681 }
01682 else
01683 {
01684
01685 A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
01686 B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
01687
01688 ASSERT( B > A );
01689
01690
01691 if( B <= 0. )
01692 {
01693 probability = 0.;
01694 }
01695 else
01696 {
01697 ASSERT( POW2( sinChiOver2 ) > 0. );
01698
01699 probability = 2.*lp/(PI* n*n*POW2( sinChiOver2 ));
01700
01701 if( A < 0. )
01702 {
01703 probability *= ellpk( -A/(B-A) );
01704 probability /= sqrt( B - A );
01705 }
01706 else
01707 {
01708 probability *= ellpk( A/B );
01709 probability /= sqrt( B );
01710 }
01711 }
01712 }
01713
01714 }
01715
01716 DEBUG_EXIT( "StarkCollTransProb_VF01()" );
01717
01718 return probability;
01719
01720 }
01721