00001
00002
00003
00004
00005
00006
00007
00008 #include "cddefines.h"
00009 #include "atmdat.h"
00010 #include "conv.h"
00011 #include "dense.h"
00012 #include "helike.h"
00013 #include "helike_cs.h"
00014 #include "hydro_vs_rates.h"
00015 #include "iso.h"
00016 #include "lines_service.h"
00017 #include "opacity.h"
00018 #include "phycon.h"
00019 #include "physconst.h"
00020 #include "rfield.h"
00021 #include "taulines.h"
00022 #include "thirdparty.h"
00023 #include "trace.h"
00024
00026 vector<double> CSTemp;
00028 multi_arr<realnum,3> HeCS;
00029
00030
00031 STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT, long nelem, long Collider, double deltaE, double osc_strength, double temp,
00032 double stat_weight, double I_energy_eV );
00033
00034
00035
00036 STATIC double collision_strength_VF01( long ipISO, long nelem, long n, long l, long lp, long s, long Collider,
00037 double ColliderCharge, double temp, double velOrEner, bool lgParamIsRedVel );
00038 STATIC double L_mix_integrand_VF01( long n, long l, long lp, double bmax, double red_vel, double an, double ColliderCharge, double alpha );
00039 STATIC double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi);
00040
00041
00042 class my_Integrand_S62
00043 {
00044 public:
00045 long nelem, Collider;
00046 double deltaE, osc_strength, temp, stat_weight, I_energy_eV;
00047
00048 double operator() (double proj_energy_overKT)
00049 {
00050 double col_str = S62_Therm_ave_coll_str( proj_energy_overKT, nelem, Collider, deltaE, osc_strength,
00051 temp, stat_weight, I_energy_eV );
00052 return col_str;
00053 }
00054 };
00055
00056 class my_Integrand_VF01_E
00057 {
00058 public:
00059 long ipISO, nelem, n, l, lp, s, Collider;
00060 double ColliderCharge, temp, velOrEner;
00061 bool lgParamIsRedVel;
00062
00063 double operator() (double EOverKT)
00064 {
00065 double col_str = collision_strength_VF01( ipISO, nelem, n, l, lp, s, Collider,
00066 ColliderCharge, temp, EOverKT * temp / TE1RYD, lgParamIsRedVel );
00067 return exp( -1.*EOverKT ) * col_str;
00068 }
00069 };
00070
00071 class my_Integrand_VF01_alpha
00072 {
00073 public:
00074 long n, l, lp;
00075 double bmax, red_vel, an, ColliderCharge;
00076
00077 double operator() (double alpha)
00078 {
00079 double integrand = L_mix_integrand_VF01( n, l, lp,
00080 bmax, red_vel, an, ColliderCharge, alpha );
00081 return integrand;
00082 }
00083 };
00084
00085
00086
00087 static const double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
00088 static const double ColliderCharge[4] = {1.0, 1.0, 1.0, 2.0};
00089
00090 void HeCollidSetup( void )
00091 {
00092
00093 long i, i1, j, nelem, ipHi, ipLo;
00094 bool lgEOL, lgHIT;
00095 FILE *ioDATA;
00096
00097 # define chLine_LENGTH 1000
00098 char chLine[chLine_LENGTH];
00099
00100 DEBUG_ENTRY( "HeCollidSetup()" );
00101
00102
00103 if( trace.lgTrace )
00104 fprintf( ioQQQ," HeCollidSetup opening he1_cs.dat:");
00105
00106 ioDATA = open_data( "he1_cs.dat", "r" );
00107
00108
00109 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00110 {
00111 fprintf( ioQQQ, " HeCollidSetup could not read first line of he1_cs.dat.\n");
00112 cdEXIT(EXIT_FAILURE);
00113 }
00114 i = 1;
00115 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00116
00117
00118
00119
00120 if( i1 !=COLLISMAGIC )
00121 {
00122 fprintf( ioQQQ,
00123 " HeCollidSetup: the version of he1_cs.dat is not the current version.\n" );
00124 fprintf( ioQQQ,
00125 " HeCollidSetup: I expected to find the number %i and got %li instead.\n" ,
00126 COLLISMAGIC, i1 );
00127 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00128 cdEXIT(EXIT_FAILURE);
00129 }
00130
00131
00132 lgHIT = false;
00133 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00134 {
00135
00136 if( chLine[0] == '#')
00137 continue;
00138
00139 lgHIT = true;
00140 lgEOL = false;
00141 char *chTemp = strtok(chLine," \t\n");
00142 while( chTemp != NULL )
00143 {
00144 CSTemp.push_back( atof(chTemp) );
00145 chTemp = strtok(NULL," \t\n");
00146 }
00147 break;
00148 }
00149 if( !lgHIT )
00150 {
00151 fprintf( ioQQQ, " HeCollidSetup could not find line in CS temperatures in he1_cs.dat.\n");
00152 cdEXIT(EXIT_FAILURE);
00153 }
00154 ASSERT( CSTemp.size() == 9U );
00155
00156
00157 {
00158 long nelem = ipHELIUM;
00159 long numLevs = iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max;
00160 HeCS.reserve( numLevs );
00161 for( long ipHi=1; ipHi < numLevs; ++ipHi )
00162 {
00163 HeCS.reserve( ipHi, ipHi );
00164 for( long ipLo=0; ipLo<ipHi; ++ipLo )
00165 HeCS.reserve( ipHi, ipLo, CSTemp.size() );
00166 }
00167 HeCS.alloc();
00168 }
00169
00170
00171 lgHIT = false;
00172 nelem = ipHELIUM;
00173 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00174 {
00175 char *chTemp;
00176
00177 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
00178 break;
00179 if( chLine[0] != '#')
00180 {
00181 lgHIT = true;
00182
00183
00184 j = 1;
00185 ipLo = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
00186 ipHi = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
00187 ASSERT( ipLo < ipHi );
00188 if( ipHi >= iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
00189 continue;
00190 else
00191 {
00192 chTemp = chLine;
00193
00194 for( long i=0; i<3; ++i )
00195 {
00196 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
00197 {
00198 fprintf( ioQQQ, " HeCollidSetup could not init cs\n" );
00199 cdEXIT(EXIT_FAILURE);
00200 }
00201 ++chTemp;
00202 }
00203
00204 for( unsigned i=0; i< CSTemp.size(); ++i )
00205 {
00206 double a;
00207 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
00208 {
00209 fprintf( ioQQQ, " HeCollidSetup could not scan cs, current indices: %li %li\n", ipHi, ipLo );
00210 cdEXIT(EXIT_FAILURE);
00211 }
00212 ++chTemp;
00213 sscanf( chTemp , "%le" , &a );
00214 HeCS[ipHi][ipLo][i] = (realnum)a;
00215 }
00216 }
00217 }
00218 }
00219
00220
00221 fclose( ioDATA );
00222
00223 return;
00224 }
00225
00226
00227 realnum HeCSInterp(long int nelem,
00228 long int ipHi,
00229 long int ipLo,
00230 long int Collider )
00231 {
00232 realnum cs, factor1;
00233
00234
00235
00236 const char *where = " ";
00237
00238 DEBUG_ENTRY( "HeCSInterp()" );
00239
00240 if( !iso_ctrl.lgColl_excite[ipHE_LIKE] )
00241 {
00242 return (realnum)1E-10;
00243 }
00244
00245 if( nelem == ipHELIUM )
00246 {
00247
00248 cs = AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
00249 }
00250 else
00251 {
00252
00253 cs = IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
00254 }
00255
00256 ASSERT( cs >= 0.f );
00257
00258
00259
00260
00261 ASSERT( factor1 >= 0.f || nelem!=ipHELIUM );
00262 if( factor1 < 0.f )
00263 {
00264 ASSERT( iso_ctrl.lgCS_Vriens[ipHE_LIKE] );
00265
00266 factor1 = 1.f;
00267 }
00268
00269
00270
00271 cs *= factor1;
00272
00273 {
00274
00275 enum {DEBUG_LOC=false};
00276
00277
00278 if( DEBUG_LOC && ( nelem==ipOXYGEN ) && (cs > 0.f) && (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == 2)
00279 && ( iso_sp[ipHE_LIKE][nelem].st[ipLo].n() <= 2 ) )
00280 fprintf(ioQQQ,"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n",
00281 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(), iso_sp[ipHE_LIKE][nelem].st[ipLo].S() ,
00282 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(), iso_sp[ipHE_LIKE][nelem].st[ipHi].n() ,
00283 iso_sp[ipHE_LIKE][nelem].st[ipHi].S(), iso_sp[ipHE_LIKE][nelem].st[ipHi].l() , cs,where);
00284 }
00285
00286 return MAX2(cs, 1.e-10f);
00287 }
00288
00289 realnum AtomCSInterp(long int nelem,
00290 long int ipHi,
00291 long int ipLo,
00292 realnum *factor1,
00293 const char **where,
00294 long int Collider )
00295 {
00296 long ipArray;
00297 realnum cs;
00298
00299 DEBUG_ENTRY( "AtomCSInterp()" );
00300
00301 ASSERT( nelem == ipHELIUM );
00302
00303
00304 cs = -1.f;
00305
00306
00307
00308
00309 *factor1 = -1.f;
00310
00311
00312
00313
00314
00315
00316 if( ipLo >= ipHe2p3P0 && ipHi <= ipHe2p3P2 && Collider==ipELECTRON )
00317 {
00318 *factor1 = 1.f;
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P1 )
00332 {
00333 cs = 1.2f;
00334 }
00335 else if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P2 )
00336 {
00337 cs = 2.1f;
00338 }
00339 else if( ipLo == ipHe2p3P1 && ipHi == ipHe2p3P2 )
00340 {
00341 cs = 6.0f;
00342 }
00343 else
00344 {
00345 cs = 1.0f;
00346 TotalInsanity();
00347 }
00348
00349 *where = "Berr ";
00350
00351 }
00352
00353
00354
00355 else if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() <= 5 &&
00356 ( ipHi < iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) &&
00357 nelem==ipHELIUM && HeCS[ipHi][ipLo][0] >= 0.f && Collider== ipELECTRON )
00358 {
00359 ASSERT( *factor1 == -1.f );
00360 ASSERT( ipLo < ipHi );
00361 ASSERT( ipHe2p3P0 == 3 );
00362
00363
00364 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
00365 {
00366
00367
00368
00369 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
00370
00371
00372 ASSERT( ipHi > ipHe2p3P2 );
00373 }
00374
00375 else if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00376 {
00377 ASSERT( ipLo < ipHe2p3P0 );
00378
00379 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
00380 }
00381
00382 else
00383 {
00384 *factor1 = 1.f;
00385 }
00386
00387
00388
00389
00390
00391
00392
00393
00394 if( phycon.alogte <= CSTemp[0] )
00395 {
00396 cs = HeCS[ipHi][ipLo][0];
00397 }
00398 else if( phycon.alogte >= CSTemp.back() )
00399 {
00400 cs = HeCS[ipHi][ipLo][CSTemp.size()-1];
00401 }
00402 else
00403 {
00404 realnum flow;
00405
00406
00407 ipArray = (long)((phycon.alogte - CSTemp[0])/(CSTemp[1]-CSTemp[0]));
00408 ASSERT( (unsigned)ipArray < CSTemp.size() );
00409 ASSERT( ipArray >= 0 );
00410
00411 flow = (realnum)( (phycon.alogte - CSTemp[ipArray])/
00412 (CSTemp[ipArray+1]-CSTemp[ipArray]));
00413 ASSERT( (flow >= 0.f) && (flow <= 1.f) );
00414
00415 cs = HeCS[ipHi][ipLo][ipArray] * (1.f-flow) +
00416 HeCS[ipHi][ipLo][ipArray+1] * flow;
00417 }
00418
00419 *where = "Bray ";
00420
00421
00422 if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() )
00423
00424 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00425 else
00426 {
00427
00428 cs *= (realnum)iso_ctrl.lgColl_excite[ipHE_LIKE];
00429 }
00430
00431 ASSERT( cs >= 0.f );
00432
00433 }
00434
00435 else if( (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() ) &&
00436 (iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == iso_sp[ipHE_LIKE][nelem].st[ipLo].S() ) )
00437 {
00438 ASSERT( *factor1 == -1.f );
00439 *factor1 = 1.f;
00440
00441
00442 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max );
00443
00444 if( (iso_sp[ipHE_LIKE][nelem].st[ipLo].l() <=2) &&
00445 abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1 )
00446 {
00447
00448
00449
00450 cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
00451 }
00452 else if( iso_ctrl.lgCS_Vrinceanu[ipHE_LIKE] )
00453 {
00454 if( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() >=3 &&
00455 iso_sp[ipHE_LIKE][nelem].st[ipHi].l() >=3 )
00456 {
00457
00458
00459
00460 cs = (realnum)CS_l_mixing_VF01( ipHE_LIKE,
00461 nelem,
00462 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
00463 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
00464 iso_sp[ipHE_LIKE][nelem].st[ipHi].l(),
00465 iso_sp[ipHE_LIKE][nelem].st[ipLo].S(),
00466 (double)phycon.te,
00467 Collider );
00468 }
00469 else
00470 {
00471 cs = 0.f;
00472 }
00473 }
00474
00475 else if( abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1)
00476 {
00477
00478
00479 cs = (realnum)CS_l_mixing_PS64(
00480 nelem,
00481 iso_sp[ipHE_LIKE][nelem].st[ipLo].lifetime(),
00482 nelem+1.-ipHE_LIKE,
00483 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
00484 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
00485 iso_sp[ipHE_LIKE][nelem].st[ipHi].g(),
00486 Collider);
00487 }
00488 else
00489 {
00490
00491 cs = 0.f;
00492 }
00493
00494
00495 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
00496 {
00497 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
00498 }
00499
00500
00501 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00502 {
00503 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
00504 }
00505
00506 *where = "lmix ";
00507 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00508 }
00509
00510
00511 else if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() != iso_sp[ipHE_LIKE][nelem].st[ipLo].n() )
00512 {
00513 ASSERT( *factor1 == -1.f );
00514
00515
00516
00517
00518 if( iso_ctrl.lgCS_Vriens[ipHE_LIKE] )
00519 {
00520
00521
00522 cs = (realnum)CS_VS80(
00523 ipHE_LIKE,
00524 nelem,
00525 ipHi,
00526 ipLo,
00527 iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul(),
00528 phycon.te,
00529 Collider );
00530
00531 *factor1 = 1.f;
00532 *where = "Vriens";
00533 }
00534 else if( iso_ctrl.lgCS_None[ipHE_LIKE] )
00535 {
00536 cs = 0.f;
00537 *factor1 = 1.f;
00538 *where = "no gb";
00539 }
00540 else if( iso_ctrl.nCS_new[ipHE_LIKE] )
00541 {
00542 *factor1 = 1.f;
00543
00544
00545
00546
00547
00548 if( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul() > 1. )
00549 {
00550
00551 double x =
00552 log10(MAX2(34.7,iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN()));
00553
00554 if( iso_ctrl.nCS_new[ipHE_LIKE] == 1 )
00555 {
00556
00557
00558 if( x < 4.5 )
00559 {
00560
00561 cs = (realnum)pow( 10. , -1.45*x + 6.75);
00562 }
00563 else
00564 {
00565
00566 cs = (realnum)pow( 10. , -3.33*x+15.15);
00567 }
00568 }
00569 else if( iso_ctrl.nCS_new[ipHE_LIKE] == 2 )
00570 {
00571
00572 cs = (realnum)pow( 10. , -2.3*x+10.3);
00573 }
00574 else
00575 TotalInsanity();
00576 }
00577 else
00578 {
00579
00580 if( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN() < 25119.f )
00581 {
00582 cs = 0.631f;
00583 }
00584 else
00585 {
00586 cs = (realnum)pow(10.,
00587 -3.*log10(iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN())+12.8);
00588 }
00589 }
00590
00591 *where = "newgb";
00592 }
00593 else
00594 TotalInsanity();
00595
00596
00597 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
00598 {
00599 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
00600 }
00601
00602
00603 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00604 {
00605 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
00606 }
00607
00608
00609 cs *= (realnum)iso_ctrl.lgColl_excite[ipHE_LIKE];
00610
00611 }
00612 else
00613 {
00614
00615
00616
00617 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() != iso_sp[ipHE_LIKE][nelem].st[ipLo].S() );
00618 cs = 0.f;
00619 *factor1 = 1.f;
00620 }
00621
00622 ASSERT( cs >= 0.f );
00623
00624 return(cs);
00625
00626 }
00627
00628
00629 realnum IonCSInterp( long nelem , long ipHi , long ipLo, realnum *factor1, const char **where, long Collider )
00630 {
00631 realnum cs;
00632
00633 DEBUG_ENTRY( "IonCSInterp()" );
00634
00635 ASSERT( nelem > ipHELIUM );
00636 ASSERT( nelem < LIMELM );
00637
00638
00639 cs = -1.f;
00640
00641
00642
00643
00644 *factor1 = -1.f;
00645
00646
00647
00648
00649
00650
00651
00652
00653 if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n()==2
00654 && iso_sp[ipHE_LIKE][nelem].st[ipLo].n()<=2 && Collider==ipELECTRON)
00655 {
00656 *where = "Zhang";
00657 *factor1 = 1.;
00658
00659
00660 if( ipLo == ipHe1s1S )
00661 {
00662 switch( ipHi )
00663 {
00664 case 1:
00665 cs = 0.25f/POW2(nelem+1.f);
00666 break;
00667 case 2:
00668 cs = 0.4f/POW2(nelem+1.f);
00669 break;
00670 case 3:
00671 cs = 0.15f/POW2(nelem+1.f);
00672 break;
00673 case 4:
00674 cs = 0.45f/POW2(nelem+1.f);
00675 break;
00676 case 5:
00677 cs = 0.75f/POW2(nelem+1.f);
00678 break;
00679 case 6:
00680 cs = 1.3f/POW2(nelem+1.f);
00681 break;
00682 default:
00683 TotalInsanity();
00684 break;
00685 }
00686 cs *= (realnum)iso_ctrl.lgColl_excite[ipHE_LIKE];
00687 }
00688
00689 else if( ipLo == ipHe2s3S )
00690 {
00691 switch( ipHi )
00692 {
00693 case 2:
00694 cs = 2.75f/POW2(nelem+1.f);
00695 break;
00696 case 3:
00697 cs = 60.f/POW2(nelem+1.f);
00698 break;
00699 case 4:
00700 cs = 180.f/POW2(nelem+1.f);
00701 break;
00702 case 5:
00703 cs = 300.f/POW2(nelem+1.f);
00704 break;
00705 case 6:
00706 cs = 5.8f/POW2(nelem+1.f);
00707 break;
00708 default:
00709 TotalInsanity();
00710 break;
00711 }
00712 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00713 }
00714
00715 else if( ipLo == ipHe2s1S )
00716 {
00717 switch( ipHi )
00718 {
00719 case 3:
00720 cs = 0.56f/POW2(nelem+1.f);
00721 break;
00722 case 4:
00723 cs = 1.74f/POW2(nelem+1.f);
00724 break;
00725 case 5:
00726 cs = 2.81f/POW2(nelem+1.f);
00727 break;
00728 case 6:
00729 cs = 190.f/POW2(nelem+1.f);
00730 break;
00731 default:
00732 TotalInsanity();
00733 break;
00734 }
00735 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00736 }
00737
00738 else if( ipLo == ipHe2p3P0 )
00739 {
00740 switch( ipHi )
00741 {
00742 case 4:
00743 cs = 8.1f/POW2(nelem+1.f);
00744 break;
00745 case 5:
00746 cs = 8.2f/POW2(nelem+1.f);
00747 break;
00748 case 6:
00749 cs = 3.9f/POW2(nelem+1.f);
00750 break;
00751 default:
00752 TotalInsanity();
00753 break;
00754 }
00755 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00756 }
00757
00758 else if( ipLo == ipHe2p3P1 )
00759 {
00760 switch( ipHi )
00761 {
00762 case 5:
00763 cs = 30.f/POW2(nelem+1.f);
00764 break;
00765 case 6:
00766 cs = 11.7f/POW2(nelem+1.f);
00767 break;
00768 default:
00769 TotalInsanity();
00770 break;
00771 }
00772 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00773 }
00774
00775 else if( ipLo == ipHe2p3P2 )
00776 {
00777
00778 cs = 19.5f/POW2(nelem+1.f) * (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00779 }
00780 else
00781 TotalInsanity();
00782
00783
00784 }
00785
00786
00787 else if( (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() ) &&
00788 (iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == iso_sp[ipHE_LIKE][nelem].st[ipLo].S() ) )
00789 {
00790 ASSERT( *factor1 == -1.f );
00791 *factor1 = 1.f;
00792
00793 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max );
00794
00795 if( (iso_sp[ipHE_LIKE][nelem].st[ipLo].l() <=2) &&
00796 abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1 )
00797 {
00798
00799
00800
00801 cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
00802 }
00803 else if( iso_ctrl.lgCS_Vrinceanu[ipHE_LIKE] )
00804 {
00805 if( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() >=3 &&
00806 iso_sp[ipHE_LIKE][nelem].st[ipHi].l() >=3 )
00807 {
00808
00809
00810
00811 cs = (realnum)CS_l_mixing_VF01( ipHE_LIKE,
00812 nelem,
00813 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
00814 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
00815 iso_sp[ipHE_LIKE][nelem].st[ipHi].l(),
00816 iso_sp[ipHE_LIKE][nelem].st[ipLo].S(),
00817 (double)phycon.te,
00818 Collider );
00819 }
00820 else
00821 {
00822 cs = 0.f;
00823 }
00824 }
00825
00826 else if( abs(iso_sp[ipHE_LIKE][nelem].st[ipHi].l() - iso_sp[ipHE_LIKE][nelem].st[ipLo].l())== 1)
00827 {
00828
00829
00830 cs = (realnum)CS_l_mixing_PS64(
00831 nelem,
00832 iso_sp[ipHE_LIKE][nelem].st[ipLo].lifetime(),
00833 nelem+1.-ipHE_LIKE,
00834 iso_sp[ipHE_LIKE][nelem].st[ipLo].n(),
00835 iso_sp[ipHE_LIKE][nelem].st[ipLo].l(),
00836 iso_sp[ipHE_LIKE][nelem].st[ipHi].g(),
00837 Collider);
00838 }
00839 else
00840 {
00841
00842 cs = 0.f;
00843 }
00844
00845
00846 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
00847 {
00848 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
00849 }
00850
00851 *where = "lmix ";
00852 cs *= (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
00853 }
00854
00855
00856 else if( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() != iso_sp[ipHE_LIKE][nelem].st[ipLo].n() )
00857 {
00858 if( iso_ctrl.lgCS_Vriens[ipHE_LIKE] )
00859 {
00860
00861
00862
00863 cs = (realnum)CS_VS80(
00864 ipHE_LIKE,
00865 nelem,
00866 ipHi,
00867 ipLo,
00868 iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul(),
00869 phycon.te,
00870 Collider );
00871
00872 *factor1 = 1.f;
00873 *where = "Vriens";
00874 }
00875 else
00876 {
00877
00878
00879 fixit();
00880
00881 cs = 0.f;
00882 *where = "hydro";
00883 }
00884 }
00885 else
00886 {
00887
00888
00889
00890 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == iso_sp[ipHE_LIKE][nelem].st[ipLo].n() );
00891 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() != iso_sp[ipHE_LIKE][nelem].st[ipLo].S() );
00892 cs = 0.f;
00893 *where = "spin ";
00894 }
00895
00896 ASSERT( cs >= 0.f );
00897
00898 return(cs);
00899 }
00900
00901
00902
00903
00904 double CS_l_mixing_S62(
00905 long ipISO,
00906 long nelem ,
00907 long ipLo ,
00908 long ipHi ,
00909 double temp,
00910 long Collider)
00911 {
00912
00913 double coll_str;
00914
00915 DEBUG_ENTRY( "CS_l_mixing_S62()" );
00916
00917 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
00918 {
00919 return 0.;
00920 }
00921
00922 my_Integrand_S62 func;
00923
00924 func.temp = temp;
00925 func.stat_weight = iso_sp[ipISO][nelem].st[ipLo].g();
00926
00927
00928
00929 func.deltaE = iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg()/EN1EV;
00930 func.I_energy_eV = EVRYD*iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd;
00931 func.Collider = Collider;
00932 func.nelem = nelem;
00933
00934 ASSERT( TRANS_PROB_CONST*POW2(func.deltaE/WAVNRYD/EVRYD) > 0. );
00935
00936 func.osc_strength = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul()/
00937 (TRANS_PROB_CONST*POW2(func.deltaE/WAVNRYD/EVRYD));
00938
00939 Integrator<my_Integrand_S62,Gaussian32> S62;
00940
00941 coll_str = S62.sum( 0., 1., func );
00942 coll_str += S62.sum( 1., 10., func );
00943 ASSERT( coll_str > 0. );
00944
00945 return coll_str;
00946 }
00947
00948
00949 STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT, long nelem, long Collider, double deltaE, double osc_strength, double temp,
00950 double stat_weight, double I_energy_eV )
00951 {
00952
00953 double integrand, cross_section, coll_str, zOverB2;
00954 double betaone, zeta_of_betaone, cs2;
00955 double Dubya, proj_energy;
00956 double reduced_mass;
00957
00958 DEBUG_ENTRY( "S62_Therm_ave_coll_str()" );
00959
00960
00961 proj_energy = proj_energy_overKT * phycon.te / EVDEGK;
00962
00963 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
00964 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
00965
00966
00967
00968 proj_energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
00969
00970
00971 proj_energy += deltaE;
00972 Dubya = 0.5*(2.*proj_energy-deltaE);
00973 ASSERT( Dubya > 0. );
00974
00975
00976
00977 ASSERT( I_energy_eV > 0. );
00978 ASSERT( osc_strength > 0. );
00979
00980
00981 zOverB2 = 0.5*POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
00982
00983 ASSERT( zOverB2 > 0. );
00984
00985 if( zOverB2 > 100. )
00986 {
00987 betaone = sqrt( 1./zOverB2 );
00988 }
00989 else if( zOverB2 < 0.54 )
00990 {
00991
00992 betaone = (1./3.)*(log(PI)-log(zOverB2)+1.28);
00993 if( betaone > 2.38 )
00994 {
00995
00996 betaone += 0.5*(log(PI)-log(zOverB2));
00997 betaone *= 0.5;
00998 }
00999 }
01000 else
01001 {
01002 long ip_zOverB2 = 0;
01003 double zetaOVERbeta2[101] = {
01004 1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
01005 7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
01006 4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
01007 3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
01008 2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
01009 1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
01010 1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
01011 7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
01012 4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
01013 3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
01014 2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
01015 1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
01016 7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
01017
01018 ASSERT( zOverB2 >= zetaOVERbeta2[100] );
01019
01020
01021 for( long i=0; i< 100; ++i )
01022 {
01023 if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
01024 {
01025
01026 ip_zOverB2 = i;
01027 break;
01028 }
01029 }
01030
01031 ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
01032
01033 betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) /
01034 (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
01035 (pow(10., ((double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((double)ip_zOverB2)/100. - 1.)) +
01036 pow(10., (double)ip_zOverB2/100. - 1.);
01037
01038 }
01039
01040 zeta_of_betaone = zOverB2 * POW2(betaone);
01041
01042
01043 cs2 = 0.5*zeta_of_betaone + betaone * bessel_k0(betaone) * bessel_k1(betaone);
01044
01045
01046 cross_section = cs2;
01047
01048
01049 cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
01050
01051
01052 coll_str = ConvCrossSect2CollStr( cross_section * PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM, stat_weight, proj_energy/EVRYD, reduced_mass );
01053
01054 integrand = exp( -1.*(proj_energy-deltaE)*EVDEGK/temp ) * coll_str;
01055 return integrand;
01056 }
01057
01058
01059 double CS_l_mixing_PS64(
01060 long nelem,
01061 double tau,
01062 double target_charge,
01063 long n,
01064 long l,
01065 double gHi,
01066 long Collider)
01067 {
01068
01069
01070 double cs, Dnl,
01071 TwoLogDebye, TwoLogRc1,
01072 factor1, factor2,
01073 bestfactor, factorpart,
01074 reduced_mass, reduced_mass_2_emass,
01075 rate;
01076 const double ChargIncoming = ColliderCharge[Collider];
01077
01078 DEBUG_ENTRY( "CS_l_mixing_PS64()" );
01079
01080
01081
01082
01083
01084 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
01085 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
01086
01087 reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;
01088
01089
01090
01091
01092
01093
01094
01095
01096 TwoLogDebye = 1.181 + log10( phycon.te / MIN2(1e10 , dense.eden ) );
01097
01098
01099
01100 TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass );
01101
01102 #if 0
01103
01104
01105 if( deltaE > 0. )
01106 TwoLogRc2 = 2. * log10( 1.12 * H_BAR * v / deltaE );
01107 else
01108 TwoLogRc2 = BIGDOUBLE;
01109 #endif
01110
01111
01112 Dnl = POW2( ChargIncoming / target_charge) * 6. * POW2( (double)n) *
01113 ( POW2((double)n) - POW2((double)l) - l - 1);
01114
01115 ASSERT( Dnl > 0. );
01116 ASSERT( phycon.te / Dnl / reduced_mass_2_emass > 0. );
01117
01118 factorpart = (11.54 + log10( phycon.te / Dnl / reduced_mass_2_emass ) );
01119
01120 if( (factor1 = factorpart + TwoLogDebye) <= 0.)
01121 factor1 = BIGDOUBLE;
01122 if( (factor2 = factorpart + TwoLogRc1) <= 0.)
01123 factor2 = BIGDOUBLE;
01124
01125
01126 bestfactor = MIN2(factor1,factor2);
01127
01128 ASSERT( bestfactor > 0. );
01129
01130
01131 if( bestfactor > 100. )
01132 {
01133 return SMALLFLOAT;
01134 }
01135
01136
01137 rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnl / phycon.sqrte * bestfactor;
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147 if( l > 0 )
01148 rate /=2;
01149
01150
01151
01152
01153 cs = rate / ( COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) *
01154 phycon.sqrte * gHi;
01155
01156 ASSERT( cs > 0. );
01157
01158 return cs;
01159 }
01160
01161
01162
01163
01164
01165
01166 double CS_l_mixing_VF01(long int ipISO,
01167 long int nelem,
01168 long int n,
01169 long int l,
01170 long int lp,
01171 long int s,
01172 double temp,
01173 long int Collider )
01174 {
01175
01176 double coll_str;
01177
01178 DEBUG_ENTRY( "CS_l_mixing_VF01()" );
01179
01180 my_Integrand_VF01_E func;
01181 func.ipISO = ipISO;
01182 func.nelem = nelem;
01183 func.n = n;
01184 func.l = l;
01185 func.lp = lp;
01186 func.s = s;
01187 func.Collider = Collider;
01188 func.temp = temp;
01189 func.ColliderCharge = ColliderCharge[Collider];
01190 func.lgParamIsRedVel = false;
01191 ASSERT( func.ColliderCharge > 0. );
01192
01193 Integrator<my_Integrand_VF01_E, Gaussian32> VF01_E;
01194
01195
01196 if( ipISO > ipH_LIKE )
01197 {
01198 ASSERT( l != 0 );
01199 ASSERT( lp != 0 );
01200 }
01201
01202 if( !iso_ctrl.lgCS_therm_ave[ipISO] )
01203 {
01204
01205
01206
01207
01208 if( (dense.eden > 10000.) && (dense.eden < 1E10 ) )
01209 {
01210 coll_str = VF01_E.sum( 0.0, 6.0, func );
01211 }
01212 else
01213 {
01214
01215 coll_str = collision_strength_VF01( ipISO, nelem, n, l, lp, s, Collider,
01216 ColliderCharge[Collider], temp, temp/TE1RYD, false );
01217 }
01218 }
01219 else
01220 {
01221
01222 coll_str = VF01_E.sum( 0.0, 1.0, func );
01223 coll_str += VF01_E.sum( 1.0, 10.0, func );
01224 }
01225
01226 return coll_str;
01227 }
01228
01229 STATIC double collision_strength_VF01( long ipISO, long nelem, long n, long l, long lp, long s, long Collider, double ColliderCharge, double temp, double velOrEner, bool lgParamIsRedVel )
01230 {
01231 double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
01232 double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
01233 double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
01234 double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 ;
01235 long ipLo, ipHi;
01236
01237 DEBUG_ENTRY( "collision_strength_VF01()" );
01238
01239 ASSERT( n > 0 );
01240
01241
01242 ipLo = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][s];
01243 ipHi = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][lp][s];
01244 stat_weight = iso_sp[ipISO][nelem].st[ipLo].g();
01245
01246
01247 if( ipISO > ipH_LIKE )
01248 {
01249
01250 ASSERT( l > 0 );
01251 ASSERT( lp > 0 );
01252 }
01253
01254
01255 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
01256 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
01257 ASSERT( reduced_mass > 0. );
01258
01259
01260
01261 aveRadius = (BOHR_RADIUS_CM/((double)nelem+1.-(double)ipISO))*POW2((double)n);
01262 ASSERT( aveRadius < 1.e-4 );
01263
01264
01265
01266 ASSERT( aveRadius > 3.9/LIMELM * BOHR_RADIUS_CM );
01267
01268
01269
01270 RMSv = ((double)nelem+1.-(double)ipISO)*POW2(ELEM_CHARGE_ESU)/(double)n/H_BAR;
01271 ASSERT( RMSv > 0. );
01272
01273 ASSERT( ColliderMass[Collider] > 0. );
01274
01275 if( lgParamIsRedVel )
01276 {
01277
01278 reduced_vel = velOrEner;
01279
01280
01281 E_Proj_Ryd = 0.5 * POW2( reduced_vel * RMSv ) * ColliderMass[Collider] *
01282 PROTON_MASS / EN1RYD;
01283 }
01284 else
01285 {
01286
01287 E_Proj_Ryd = velOrEner;
01288 reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/ColliderMass[Collider]/PROTON_MASS )/RMSv;
01289 }
01290
01291
01292 ASSERT( reduced_vel > 1.e-10 );
01293 ASSERT( reduced_vel < 1.e10 );
01294
01295
01296 ConstantFactors = 4.5*PI*POW2(ColliderCharge*aveRadius/reduced_vel);
01297
01298
01299 reduced_b_min = 1.5 * ColliderCharge / reduced_vel;
01300 ASSERT( reduced_b_min > 1.e-10 );
01301 ASSERT( reduced_b_min < 1.e10 );
01302
01303 if( ipISO == ipH_LIKE )
01304 {
01305
01306 reduced_b_max = sqrt( BOLTZMANN*temp/ColliderCharge/dense.eden )
01307 / (PI2*ELEM_CHARGE_ESU)/aveRadius;
01308 }
01309 else if( ipISO == ipHE_LIKE )
01310 {
01311 quantum_defect1 = (double)n- (double)nelem/sqrt(iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd);
01312 quantum_defect2 = (double)n- (double)nelem/sqrt(iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
01313
01314
01315 ASSERT( fabs(quantum_defect1) < 1.0 );
01316 ASSERT( fabs(quantum_defect1) > 0.0 );
01317 ASSERT( fabs(quantum_defect2) < 1.0 );
01318 ASSERT( fabs(quantum_defect2) > 0.0 );
01319
01320
01321 omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
01322
01323 omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)lp/(double)n)) / POW3( (double)n ) / (double)lp );
01324
01325 omega_qd = 0.5*( omega_qd1 + omega_qd2 );
01326
01327 ASSERT( omega_qd > 0. );
01328
01329 reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
01330 }
01331 else
01332
01333 TotalInsanity();
01334
01335 reduced_b_max = MAX2( reduced_b_max, reduced_b_min );
01336 ASSERT( reduced_b_max > 0. );
01337
01338
01339 my_Integrand_VF01_alpha func;
01340 func.n = n;
01341 func.l = l;
01342 func.lp = lp;
01343 func.red_vel = reduced_vel;
01344 func.an = aveRadius;
01345 func.ColliderCharge = ColliderCharge;
01346 func.bmax = reduced_b_max*aveRadius;
01347 Integrator<my_Integrand_VF01_alpha, Gaussian32> VF01_alpha;
01348
01349 alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max);
01350 alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min);
01351
01352 ASSERT( alphamin > 0. );
01353 ASSERT( alphamax > 0. );
01354
01355 alphamin = MAX2( alphamin, 1.e-30 );
01356 alphamax = MAX2( alphamax, 1.e-20 );
01357
01358 CSIntegral = 0.;
01359
01360 if( alphamax > alphamin )
01361 {
01362
01363 step = (alphamax - alphamin)/5.;
01364 alpha1 = alphamin;
01365 CSIntegral += VF01_alpha.sum( alpha1, alpha1+step, func );
01366 CSIntegral += VF01_alpha.sum( alpha1+step, alpha1+4.*step, func );
01367 }
01368
01369
01370 cross_section = ConstantFactors * CSIntegral;
01371
01372
01373 coll_str = ConvCrossSect2CollStr( cross_section, stat_weight, E_Proj_Ryd, reduced_mass );
01374
01375 coll_str = MAX2( SMALLFLOAT, coll_str);
01376
01377 return coll_str;
01378 }
01379
01380 STATIC double L_mix_integrand_VF01( long n, long l, long lp, double bmax, double red_vel, double an, double ColliderCharge, double alpha )
01381 {
01382 double integrand, deltaPhi, b;
01383
01384 DEBUG_ENTRY( "L_mix_integrand_VF01()" );
01385
01386 ASSERT( alpha >= 1.e-30 );
01387 ASSERT( bmax > 0. );
01388 ASSERT( red_vel > 0. );
01389
01390
01391 b = 1.5*ColliderCharge*an/(red_vel * alpha);
01392
01393 if( b < bmax )
01394 {
01395 deltaPhi = -1.*PI + 2.*asin(b/bmax);
01396 }
01397 else
01398 {
01399 deltaPhi = 0.;
01400 }
01401 integrand = 1./alpha/alpha/alpha;
01402 integrand *= StarkCollTransProb_VF01( n, l, lp, alpha, deltaPhi );
01403
01404 return integrand;
01405 }
01406
01407 STATIC double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi)
01408 {
01409 double probability;
01410 double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
01411
01412 DEBUG_ENTRY( "StarkCollTransProb_VF01()" );
01413
01414 ASSERT( alpha > 0. );
01415
01416
01417 cosU1 = 2.*POW2((double)l/(double)n) - 1.;
01418 cosU2 = 2.*POW2((double)lp/(double)n) - 1.;
01419
01420 sinU1 = sqrt( 1. - cosU1*cosU1 );
01421 sinU2 = sqrt( 1. - cosU2*cosU2 );
01422
01423
01424 cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
01425 sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
01426 cosChi = 2. * POW2( cosChiOver2 ) - 1.;
01427
01428 if( l == 0 )
01429 {
01430 if( -1.*cosU2 - cosChi < 0. )
01431 {
01432 probability = 0.;
01433 }
01434 else
01435 {
01436
01437 ASSERT( sinChiOver2 > 0. );
01438 ASSERT( sinChiOver2*sinChiOver2 > POW2((double)lp/(double)n) );
01439
01440
01441 probability = (double)lp/(POW2((double)n)*sinChiOver2*sqrt( POW2(sinChiOver2) - POW2((double)lp/(double)n) ) );
01442 }
01443 }
01444 else
01445 {
01446 double OneMinusCosChi = 1. - cosChi;
01447
01448 if( OneMinusCosChi == 0. )
01449 {
01450 double hold = sin( deltaPhi / 2. );
01451
01452 OneMinusCosChi = 8. * alpha * alpha * POW2( hold );
01453 }
01454
01455 if( OneMinusCosChi == 0. )
01456 {
01457 probability = 0.;
01458 }
01459 else
01460 {
01461
01462 A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
01463 B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
01464
01465 ASSERT( B > A );
01466
01467
01468 if( B <= 0. )
01469 {
01470 probability = 0.;
01471 }
01472 else
01473 {
01474 ASSERT( POW2( sinChiOver2 ) > 0. );
01475
01476 probability = 2.*lp/(PI* n*n*POW2( sinChiOver2 ));
01477
01478 if( A < 0. )
01479 {
01480 probability *= ellpk( -A/(B-A) );
01481 probability /= sqrt( B - A );
01482 }
01483 else
01484 {
01485 probability *= ellpk( A/B );
01486 probability /= sqrt( B );
01487 }
01488 }
01489 }
01490
01491 }
01492
01493 return probability;
01494
01495 }