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