00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "cddefines.h"
00014 #include "physconst.h"
00015 #include "taulines.h"
00016 #include "path.h"
00017 #include "dense.h"
00018 #include "trace.h"
00019 #include "hydro_bauman.h"
00020 #include "iso.h"
00021 #include "helike.h"
00022 #include "helike_einsta.h"
00023
00024
00025
00026
00027
00028 static double ***TransProbs;
00029
00030
00031
00032 static double ritoa( long li, long lf, long nelem, double k, double RI2 );
00033
00034
00035 static double ForbiddenAuls( long ipHi, long ipLo, long nelem );
00036
00037
00038 static double CalculateRadialIntegral( long nelem, long ipLo, long ipHi );
00039
00040 static double RadialIntegrand( double r12 );
00041
00042 static double WaveFunction( double r12 );
00043
00044 static long RI_ipHi, RI_ipLo, RI_ipLev;
00045 static long globalZ;
00046
00047
00048
00049
00050 static double vJint , zJint;
00051
00052
00053 static double Jint( double theta )
00054 {
00055
00056 double
00057 d0 = ( 1.0 / PI ),
00058 d1 = vJint * theta,
00059 d2 = zJint * sin(theta),
00060 d3 = (d1 - d2),
00061 d4 = cos(d3),
00062 d5 = (d0 * d4);
00063
00064 return( d5 );
00065 }
00066
00067
00068 static double AngerJ( double vv, double zz )
00069 {
00070 long int rep = 0, ddiv, divsor;
00071
00072 double y = 0.0;
00073
00074
00075
00076
00077 if( (fabs(vv)) - (int)(fabs(vv)) > 0.5 )
00078 ddiv = (int)(fabs(vv)) + 1;
00079 else
00080 ddiv = (int)(fabs(vv));
00081
00082 divsor = ((ddiv == 0) ? 1 : ddiv);
00083 vJint = vv;
00084 zJint = zz;
00085
00086 for( rep = 0; rep < divsor; rep++ )
00087 {
00088 double
00089 rl = (((double) rep)/((double) divsor)),
00090 ru = (((double) (rep+1))/((double) divsor)),
00091 x_low = (PI * rl),
00092 x_up = (PI * ru);
00093
00094 y += qg32( x_low, x_up, Jint );
00095 }
00096
00097 return( y );
00098 }
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 double scqdri(
00155 double nstar, long int l,
00156 double npstar, long int lp,
00157 double iz
00158 )
00159 {
00160 double n_c = ((2.0 * nstar * npstar ) / ( nstar + npstar ));
00161
00162
00163 double D_n = (nstar - npstar);
00164 double D_l = (double) ( l - lp );
00165 double lg = (double) ( (lp > l) ? lp : l);
00166
00167 double h = (lg/n_c);
00168 double g = h*h;
00169 double f = ( 1.0 - g );
00170 double e = (( f >= 0.0) ? sqrt( f ) : 0.0 );
00171
00172 double x = (e * D_n);
00173 double z = (-1.0 * x);
00174 double v1 = (D_n + 1.0);
00175 double v2 = (D_n - 1.0);
00176
00177 double d1,d2,d7,d8,d9,d34,d56,d6_1;
00178
00179 if( iz == 0.0 )
00180 iz += 1.0;
00181
00182 if( D_n == 0.0 )
00183 {
00184 return( -1.0 );
00185 }
00186
00187 if( D_n < 0.0 )
00188 {
00189 return( -1.0 );
00190 }
00191
00192 if( f < 0.0 )
00193 {
00194
00195
00196
00197
00198 return( -1.0 );
00199 }
00200
00201 d1 = ( 1.0 / iz );
00202
00203 d2 = (n_c * n_c)/(2.0 * D_n);
00204
00205 d34 = (1.0 - ((D_l * lg)/n_c)) * AngerJ( v1, z );
00206
00207 d56 = (1.0 + ((D_l * lg)/n_c)) * AngerJ( v2, z );
00208
00209 d6_1 = PI * D_n;
00210
00211 d7 = (2./PI) * sin( d6_1 ) * (1.0 - e);
00212
00213 d8 = d1 * d2 * ( (d34) - (d56) + d7 );
00214
00215 d9 = d8 * d8;
00216
00217 ASSERT( D_n > 0.0 );
00218 ASSERT( l >= 0 );
00219 ASSERT( lp >= 0 );
00220 ASSERT( (l == lp + 1) || ( l == lp - 1) );
00221 ASSERT( n_c != 0.0 );
00222 ASSERT( f >= 0.0 );
00223 ASSERT( d9 > 0.0 );
00224
00225 return( d9 );
00226 }
00227
00228 static double ForbiddenAuls( long ipHi, long ipLo, long nelem )
00229 {
00230 double A;
00231
00232
00233
00234 double As2nuFrom1S[28] = {1940.,1.82E+04,9.21E+04,3.30E+05,9.44E+05,2.31E+06,5.03E+06,1.00E+07,
00235 1.86E+07,3.25E+07,5.42E+07,8.69E+07,1.34E+08,2.02E+08,2.96E+08,4.23E+08,5.93E+08,8.16E+08,
00236 1.08E+09,1.43E+09,1.88E+09,2.43E+09,3.25E+09,3.95E+09,4.96E+09,6.52E+09,7.62E+09,9.94E+09};
00237
00238
00239
00240
00241
00242 double As2nuFrom3S[28] = {1.25E-06,5.53E-05,8.93E-04,8.05E-03,4.95E-02,2.33E-01,8.94E-01,2.95E+00,
00243 8.59E+00,2.26E+01,5.49E+01,1.24E+02,2.64E+02,5.33E+02,1.03E+03,1.91E+03,3.41E+03,5.91E+03,
00244 9.20E+03,1.50E+04,2.39E+04,3.72E+04,6.27E+04,8.57E+04,1.27E+05,2.04E+05,2.66E+05,4.17E+05};
00245
00246 if( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
00247 {
00248 if( nelem == ipHELIUM )
00249 {
00250
00251
00252
00253
00254 double ForbiddenHe[5] = { 1.272E-4, 51.02, 1E-20, 177.58, 0.327 };
00255
00256 A = ForbiddenHe[ipHi - 1];
00257 putError(nelem,ipHi,ipLo,IPRAD,-1);
00258 }
00259 else
00260 {
00261 switch ( (int)ipHi )
00262 {
00263 case 1:
00264
00265
00266 A = (3.9061E-7) * pow( (double)nelem+1., 10.419 ) + As2nuFrom3S[nelem-2];
00267 break;
00268 case 2:
00269 A = As2nuFrom1S[nelem-2];
00270 break;
00271 case 3:
00272 A = iso.SmallA;
00273 break;
00274 case 4:
00275
00276
00277 if( nelem <= ipARGON )
00278 {
00279 A = ( 11.431 * pow((double)nelem, 9.091) );
00280 }
00281 else
00282 {
00283
00284 A = ( 383.42 * pow((double)nelem, 7.8901) );
00285 }
00286 break;
00287 case 5:
00288
00289
00290
00291
00292
00293
00294 A = ( 0.11012 * pow((double)nelem, 7.6954) );
00295
00296 break;
00297 default:
00298 TotalInsanity();
00299 }
00300 putError(nelem,ipHi,ipLo,IPRAD,-1);
00301 }
00302
00303 putError(nelem,ipHi,ipLo,IPRAD,.01f);
00304 return A;
00305 }
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==1 &&
00316 L_(ipLo)==0 && S_(ipLo)==0 && N_(ipLo) < N_(ipHi) )
00317 {
00318 A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
00319 pow((double)N_(ipHi),2.877);
00320 putError(nelem,ipHi,ipLo,IPRAD,-1);
00321 }
00322
00323
00324 else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==0 &&
00325 L_(ipLo)==1 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
00326 {
00327 A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
00328
00329 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
00330 {
00331
00332
00333 A *= (2.*(ipLo-3)+1.0)/3.0;
00334 }
00335 putError(nelem,ipHi,ipLo,IPRAD,-1);
00336 }
00337
00338 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
00339 {
00340
00341 if( nelem == ipHELIUM )
00342 {
00343 A = 1.549;
00344 putError(nelem,ipHi,ipLo,IPRAD,-1);
00345 }
00346 else
00347 {
00348
00349
00350
00351 A= 0.1834*pow((double)nelem, 6.5735);
00352 putError(nelem,ipHi,ipLo,IPRAD,-1);
00353 }
00354 }
00355
00356 else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 )
00357 {
00358
00363 A = iso.SmallA;
00364 putError(nelem,ipHi,ipLo,IPRAD,-1);
00365 }
00366
00367 else
00368 {
00369
00370 A = iso.SmallA;
00371 putError(nelem,ipHi,ipLo,IPRAD,-1);
00372 }
00373
00374
00375 putError(nelem,ipHi,ipLo,IPRAD,.01f);
00376
00377 ASSERT( A > 0.);
00378 return A;
00379 }
00380
00381
00382 double he_1trans(
00383
00384 long nelem , double Enerwn ,
00385
00386 double Eff_nupper, long lHi, long sHi, long jHi,
00387
00388 double Eff_nlower, long lLo, long sLo, long jLo )
00389
00390
00391 {
00392
00393
00394
00395
00396 double RI2, MyRI, Aul;
00397 long nHi, nLo, ipHi, ipLo;
00398
00399
00400 ASSERT(nelem > ipHYDROGEN);
00401
00402
00403
00404 nHi = (int)(Eff_nupper + 0.4);
00405 nLo = (int)(Eff_nlower + 0.4);
00406
00407
00408 ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
00409 ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );
00410
00411 ipHi = QuantumNumbers2Index[nelem][nHi][lHi][sHi];
00412 if( (nHi==2) && (lHi==1) && (sHi==1) )
00413 {
00414 ASSERT( (jHi>=0) && (jHi<=2) );
00415 ipHi -= (2 - jHi);
00416 }
00417
00418 ipLo = QuantumNumbers2Index[nelem][nLo][lLo][sLo];
00419 if( (nLo==2) && (lLo==1) && (sLo==1) )
00420 {
00421 ASSERT( (jLo>=0) && (jLo<=2) );
00422 ipLo -= (2 - jLo);
00423 }
00424
00425 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n == nHi );
00426 if( nHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00427 {
00428 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].l == lHi );
00429 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == sHi );
00430 }
00431 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipLo].n == nLo );
00432 if( nLo <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00433 {
00434 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipLo].l == lLo );
00435 ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ipLo].s == sLo );
00436 }
00437
00438
00439 if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
00440 {
00441 Aul = -2.;
00442
00443
00444 if( nelem == ipHELIUM )
00445 {
00446 {
00447
00448
00449 enum {DEBUG_LOC=false};
00450
00451
00452 if( ( DEBUG_LOC ) && (nLo<=2) && (nHi==2) )
00453 {
00454 MyRI = CalculateRadialIntegral( ipHELIUM, ipLo, ipHi );
00455 fprintf(ioQQQ,"ipHi\t%li\tipLo\t%li\tMyRI^2\t%.2e\n",
00456 ipHi,
00457 ipLo,
00458 MyRI*MyRI);
00459 }
00460 }
00461
00462
00463
00464 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
00465 {
00466
00467 ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ) );
00468 ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ) );
00469 ASSERT( ipHi > 2 );
00470
00471 Aul = TransProbs[nelem][ipHi][ipLo];
00472
00473 putError(nelem,ipHi,ipLo,IPRAD,0.0005f);
00474 }
00475
00476 if( Aul < 0. )
00477 {
00478
00479 if( ipLo == ipHe1s1S )
00480 {
00481 ASSERT( (lHi == 1) && (sHi == 0) );
00482
00483
00484 if( nLo == 1 )
00485 Aul = (1.59208e10) / pow(Eff_nupper,3.0);
00486 ASSERT( Aul > 0.);
00487 putError(nelem,ipHi,ipLo,IPRAD,0.005f);
00488 }
00489
00490
00491
00492 else if( lHi>=2 && lLo>=2 && nHi>nLo )
00493 {
00494
00495 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
00496
00497 ASSERT( Aul > 0.);
00498
00499 if( lHi + lLo >= 7 )
00500 {
00501 putError(nelem,ipHi,ipLo,IPRAD,0.001f);
00502 }
00503 else
00504 {
00505 putError(nelem,ipHi,ipLo,IPRAD,0.01f);
00506 }
00507 }
00508
00509
00510
00511
00512
00513
00514 else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
00515 {
00516 int paramSet=0;
00517 double emisOscStr, x, a, b, c;
00518 double extrapol_Params[2][4][4][3] = {
00519
00520 {
00521 {
00522 { 0.8267396 , 1.4837624 , -0.4615955 },
00523 { 1.2738405 , 1.5841806 , -0.3022984 },
00524 { 1.6128996 , 1.6842538 , -0.2393057 },
00525 { 1.8855491 , 1.7709125 , -0.2115213 }},
00526 {
00527 { -1.4293664 , 2.3294080 , -0.0890470 },
00528 { -0.3608082 , 2.3337636 , -0.0712380 },
00529 { 0.3027974 , 2.3326252 , -0.0579008 },
00530 { 0.7841193 , 2.3320138 , -0.0497094 }},
00531 {
00532 { 1.1341403 , 3.1702435 , -0.2085843 },
00533 { 1.7915926 , 2.4942946 , -0.2266493 },
00534 { 2.1979400 , 2.2785377 , -0.1518743 },
00535 { 2.5018229 , 2.1925720 , -0.1081966 }},
00536 {
00537 { 0.0000000 , 0.0000000 , 0.0000000 },
00538 { -2.6737396 , 2.9379143 , -0.3805367 },
00539 { -1.4380124 , 2.7756396 , -0.2754625 },
00540 { -0.6630196 , 2.6887253 , -0.2216493 }},
00541 },
00542
00543 {
00544 {
00545 { 0.3075287 , 0.9087130 , -1.0387207 },
00546 { 0.687069 , 1.1485864 , -0.6627317 },
00547 { 0.9776064 , 1.3382024 , -0.5331906 },
00548 { 1.2107725 , 1.4943721 , -0.4779232 }},
00549 {
00550 { -1.3659605 , 2.3262253 , -0.0306439 },
00551 { -0.2899490 , 2.3279391 , -0.0298695 },
00552 { 0.3678878 , 2.3266603 , -0.0240021 },
00553 { 0.8427457 , 2.3249540 , -0.0194091 }},
00554 {
00555 { 1.3108281 , 2.8446367 , -0.1649923 },
00556 { 1.8437692 , 2.2399326 , -0.2583398 },
00557 { 2.1820792 , 2.0693762 , -0.1864091 },
00558 { 2.4414052 , 2.0168255 , -0.1426083 }},
00559 {
00560 { 0.0000000 , 0.0000000 , 0.0000000 },
00561 { -1.9219877 , 2.7689624 , -0.2536072 },
00562 { -0.7818065 , 2.6595150 , -0.1895313 },
00563 { -0.0665624 , 2.5955623 , -0.1522616 }},
00564 }
00565 };
00566
00567 if( lLo==0 )
00568 {
00569 paramSet = 0;
00570 }
00571 else if( lLo==1 && lHi==0 )
00572 {
00573 paramSet = 1;
00574 }
00575 else if( lLo==1 && lHi==2 )
00576 {
00577 paramSet = 2;
00578 }
00579 else if( lLo==2 )
00580 {
00581 paramSet = 3;
00582 ASSERT( lHi==1 );
00583 }
00584
00585 a = extrapol_Params[sHi][paramSet][nLo-2][0];
00586 b = extrapol_Params[sHi][paramSet][nLo-2][1];
00587 c = extrapol_Params[sHi][paramSet][nLo-2][2];
00588 ASSERT( Enerwn > 0. );
00589 x = log( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]*RYD_INF/Enerwn );
00590
00591 emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
00592 (2.*lLo+1)/(2.*lHi+1);
00593
00594 Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;
00595
00596 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
00597 {
00598 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00599 }
00600
00601 ASSERT( Aul > 0. );
00602 putError(nelem,ipHi,ipLo,IPRAD,0.01f);
00603 }
00604 else
00605 {
00606
00607 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
00608 ASSERT( RI2 > 0. );
00609
00610 Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
00611
00612
00613 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
00614 {
00615 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00616 }
00617
00618 ASSERT( Aul > 0. );
00619 putError(nelem,ipHi,ipLo,IPRAD,0.03f);
00620 }
00621 }
00622 }
00623
00624
00625 else
00626 {
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00637 {
00638
00639 ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) );
00640 ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) );
00641 ASSERT( ipHi > 2 );
00642
00643 Aul = TransProbs[nelem][ipHi][ipLo];
00644 putError(nelem,ipHi,ipLo,IPRAD,-1);
00645 }
00646
00647 if( Aul < 0. )
00648 {
00649
00650 if( nLo == nHi )
00651 {
00652
00653
00654 if( ipLo == ipHe2s3S )
00655 {
00656 if(ipHi == ipHe2p3P0)
00657 Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
00658 else if(ipHi == ipHe2p3P1)
00659 Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
00660 else if(ipHi == ipHe2p3P2)
00661 Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
00662 else
00663 {
00664 fprintf(ioQQQ,"WHOA! Stop in Helike.c");
00665 cdEXIT(EXIT_FAILURE);
00666 }
00667 }
00668
00669
00670 else if( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) )
00671 {
00672 Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
00673 }
00674
00675 else
00676 {
00677
00678 ASSERT( nLo > 2 );
00679
00680
00681 if( (lHi == 1) && (sHi == 1) && (lLo == 0) && (sLo == 1))
00682 {
00683 Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
00684 }
00685
00686 else if( (lHi == 1) && (sHi == 0) && (lLo == 2) && (sLo == 0))
00687 {
00688 Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
00689 }
00690
00691 else if( (lHi == 1) && (sHi == 0) && (lLo == 0) )
00692 {
00693 Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077);
00694 }
00695 else
00696 {
00697 Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
00698 if( (lHi >2) || (lLo > 2) )
00699 Aul *= (lHi/2.);
00700 if(lLo > 2)
00701 Aul *= (1./9.);
00702 }
00703 }
00704 ASSERT( Aul > 0.);
00705 }
00706
00707
00708 else if( ((lHi > 2) || (lLo > 2)) )
00709 {
00710 ASSERT( nHi > nLo );
00711
00712 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
00713
00714 ASSERT( Aul > 0.);
00715 }
00716
00717
00718
00719
00720
00721 else if( ( ipLo == 0 ) || ( ipLo == ipHe2s1S ) || ( ipLo == ipHe2s3S )
00722 || ( ipLo == ipHe3s3S ) || ( ipLo == ipHe4s3S ) )
00723 {
00724
00725 if( ipLo == 0 )
00726 {
00727 ASSERT( (lHi == 1) && (sHi == 0) );
00728
00729
00730
00731
00732
00733
00734 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
00735 }
00736
00737
00738 else if( ipLo == ipHe2s1S )
00739 {
00740 ASSERT( (lHi == 1) && (sHi == 0) );
00741
00742
00743 Aul = 5.0e8 * pow((double)nelem,4.) / pow((double)nHi, 2.889);
00744 }
00745
00746
00747 else
00748 {
00749 ASSERT( (lHi == 1) && (sHi == 1) );
00750
00751
00752 if( nLo == 2 )
00753 Aul = 1.5 * 3.405E8 * pow((double)nelem,4.) / pow((double)nHi, 2.883);
00754 else if( nLo == 3 )
00755 Aul = 2.5 * 4.613E7 * pow((double)nelem,4.) / pow((double)nHi, 2.672);
00756 else
00757 Aul = 3.0 * 1.436E7 * pow((double)nelem,4.) / pow((double)nHi, 2.617);
00758 }
00759
00760 ASSERT( Aul > 0.);
00761 }
00762
00763
00764 else
00765 {
00766
00767 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem));
00768
00769 Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2);
00770
00771
00772 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) && (Aul > iso.SmallA) )
00773 {
00774 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00775 }
00776
00777 }
00778 putError(nelem,ipHi,ipLo,IPRAD,-1);
00779 }
00780
00781
00782 putError(nelem,ipHi,ipLo,IPRAD,0.05f);
00783 }
00784 }
00785
00786
00787
00788
00789
00790
00791
00792
00793 else
00794 {
00795 ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
00796 Aul = ForbiddenAuls(ipHi, ipLo, nelem);
00797 ASSERT( Aul > 0. );
00798 }
00799
00800 Aul = MAX2( Aul, iso.SmallA );
00801 ASSERT( Aul >= iso.SmallA );
00802
00803
00804
00805 if( Enerwn < 0. && Aul > iso.SmallA )
00806 {
00807 fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
00808 }
00809
00810 return Aul;
00811 }
00812
00813 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
00814 {
00815
00816
00817
00818
00819
00820
00821
00822 long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
00823 double Ass, Att, Ast, Ats;
00824 double SinHi, SinLo, CosHi, CosLo;
00825 double HiMixingAngle, LoMixingAngle , error;
00826 double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
00827
00828 nHi = iso.quant_desig[ipHE_LIKE][nelem][ipHiSing].n;
00829 lHi = iso.quant_desig[ipHE_LIKE][nelem][ipHiSing].l;
00830 sHi = iso.quant_desig[ipHE_LIKE][nelem][ipHiSing].s;
00831 nLo = iso.quant_desig[ipHE_LIKE][nelem][ipLoSing].n;
00832 lLo = iso.quant_desig[ipHE_LIKE][nelem][ipLoSing].l;
00833 sLo = iso.quant_desig[ipHE_LIKE][nelem][ipLoSing].s;
00834
00835 if( sHi==1 || sLo ==1 )
00836 return;
00837 if( abs(lHi - lLo)!=1 )
00838 return;
00839 if( nLo < 2 )
00840 return;
00841 if( (lHi<=1) || (lLo<=1) )
00842 return;
00843 if( (nHi==nLo) && (lHi==1) && (lLo==2) )
00844 return;
00845 if( (nHi > nLo ) && ( lHi != 1 ) && ( lLo != 1) )
00846 return;
00847
00848 ASSERT( lHi > 0 );
00849
00850
00851 ipHiTrip = QuantumNumbers2Index[nelem][nHi][lHi][1];
00852 ipLoTrip = QuantumNumbers2Index[nelem][nLo][lLo][1];
00853
00854 if( lHi == 2 )
00855 {
00856 HiMixingAngle = 0.01;
00857 }
00858 else if( lHi == 3 )
00859 {
00860 HiMixingAngle = 0.5;
00861 }
00862 else
00863 {
00864 HiMixingAngle = PI/4.;
00865 }
00866
00867 if( lLo == 2 )
00868 {
00869 LoMixingAngle = 0.01;
00870 }
00871 else if( lLo == 3 )
00872 {
00873 LoMixingAngle = 0.5;
00874 }
00875 else
00876 {
00877 LoMixingAngle = PI/4.;
00878 }
00879
00880
00881 ASSERT( ipHiTrip > ipLoTrip );
00882 ASSERT( ipHiTrip > ipLoSing );
00883 ASSERT( ipHiSing > ipLoTrip );
00884 ASSERT( ipHiSing > ipLoSing );
00885
00886 SinHi = sin( HiMixingAngle );
00887 SinLo = sin( LoMixingAngle );
00888 CosHi = cos( HiMixingAngle );
00889 CosLo = cos( LoMixingAngle );
00890
00891 Kss = EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].EnergyWN;
00892 Ktt = EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].EnergyWN;
00893 Kst = EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].EnergyWN;
00894 Kts = EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].EnergyWN;
00895
00896 fss = EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul/TRANS_PROB_CONST/POW2(Kss)/
00897 EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gLo*EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gHi;
00898
00899 ftt = EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul/TRANS_PROB_CONST/POW2(Ktt)/
00900 EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gLo*EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gHi;
00901
00902 temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
00903 fssNew = Kss*POW2( temp );
00904 temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
00905 fttNew = Ktt*POW2( temp );
00906 temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
00907 fstNew = Kst*POW2( temp );
00908 temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
00909 ftsNew = Kts*POW2( temp );
00910
00911 Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gLo/
00912 EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].gHi;
00913
00914 Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gLo/
00915 EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].gHi;
00916
00917 Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].gLo/
00918 EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].gHi;
00919
00920 Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].gLo/
00921 EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].gHi;
00922
00923 error = fabs( ( EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul+
00924 EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul
00925
00926 )/
00927 (Ass+Ast+Ats+Att) - 1.f );
00928
00929 if( error > 0.001 )
00930 {
00931 fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
00932 ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
00933 Ass/EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul,
00934 Att/EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul,
00935 Ast/EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Aul,
00936 Ats/EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Aul );
00937 }
00938
00939 EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Aul = (float)Ass;
00940 EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Aul = (float)Att;
00941 EmisLines[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Aul = (float)Ast;
00942 EmisLines[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Aul = (float)Ats;
00943
00944 return;
00945 }
00946
00947
00948
00949 static double ritoa(long li, long lf, long nelem, double k, double RI2)
00950 {
00951
00952
00953
00954
00955
00956
00957
00958
00959 long lg;
00960 double fmean,mu,EinsteinA,w,RI2_cm;
00961
00962 mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT));
00963
00964 w = 2. * PI * k * SPEEDLIGHT;
00965
00966 RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM;
00967
00968 lg = (lf > li) ? lf : li;
00969
00970 fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));
00971
00972 EinsteinA = TRANS_PROB_CONST*k*k*fmean;
00973
00974
00975
00976 return EinsteinA;
00977 }
00978
00979
00980 void HelikeTransProbSetup( void )
00981 {
00982
00983 const int chLine_LENGTH = 1000;
00984 char chLine[chLine_LENGTH] ,
00985
00986 chFilename[FILENAME_PATH_LENGTH_2];
00987
00988 FILE *ioDATA;
00989 bool lgEOL;
00990
00991 long nelem, ipLo, ipHi, i, i1, i2, i3;
00992
00993 TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
00994
00995 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
00996 {
00997
00998 TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) );
00999
01000 for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo )
01001 {
01002 TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX );
01003 }
01004 }
01005
01006
01007
01008
01009
01010
01011 if( lgDataPathSet == true )
01012 {
01013
01014 strcpy( chFilename , chDataPath );
01015 strcat( chFilename , "he_transprob.dat" );
01016 }
01017 else
01018 {
01019
01020 strcpy( chFilename , "he_transprob.dat" );
01021 }
01022
01023 if( trace.lgTrace )
01024 fprintf( ioQQQ," HeCreate opening he_transprob.dat:");
01025
01026 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
01027 {
01028 fprintf( ioQQQ, " HeCreate could not open he_transprob.dat\n" );
01029 if( lgDataPathSet == true )
01030 fprintf( ioQQQ, " even tried path\n" );
01031
01032 if( lgDataPathSet == true )
01033 {
01034 fprintf( ioQQQ, " HeCreate could not open he_transprob.dat\n");
01035 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
01036 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
01037 fprintf( ioQQQ, " Computation can not continue without he_transprob.dat.\n");
01038
01039 path_not_set();
01040 cdEXIT(EXIT_FAILURE);
01041 }
01042 }
01043 else
01044 {
01045
01046 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01047 {
01048 fprintf( ioQQQ, " HeCreate could not read first line of he_transprob.dat.\n");
01049 puts( "[Stop in HeCreate]" );
01050 cdEXIT(EXIT_FAILURE);
01051 }
01052 i = 1;
01053 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01054 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01055 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
01056 {
01057 fprintf( ioQQQ,
01058 " HeCreate: the version of he_transprob.dat is not the current version.\n" );
01059 fprintf( ioQQQ,
01060 " HeCreate: I expected to find the number %i %i and got %li %li instead.\n" ,
01061 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
01062 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01063 puts( "[Stop in HeCreate]" );
01064 cdEXIT(EXIT_FAILURE);
01065 }
01066
01067
01068 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
01069 {
01070 for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
01071 {
01072 for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
01073 {
01074 TransProbs[nelem][ipHi][ipLo] = -1.;
01075 }
01076 }
01077 }
01078
01079 for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
01080 {
01081 char *chTemp;
01082
01083
01084 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01085 BadRead();
01086
01087 while( chLine[0]=='#' )
01088 {
01089 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01090 BadRead();
01091 }
01092
01093 i3 = 1;
01094 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
01095 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
01096
01097 if( i1<0 || i2<=i1 )
01098 {
01099 fprintf( ioQQQ, " HeCreate detected insanity in he_transprob.dat.\n");
01100 puts( "[Stop in HeCreate]" );
01101 cdEXIT(EXIT_FAILURE);
01102 }
01103
01104 chTemp = chLine;
01105
01106
01107 for( i=0; i<1; ++i )
01108 {
01109 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
01110 {
01111 fprintf( ioQQQ, " HeCreate could not init he_transprob\n" );
01112 puts( "[Stop in HeCreate]" );
01113 cdEXIT(EXIT_FAILURE);
01114 }
01115 ++chTemp;
01116 }
01117
01118
01119 for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
01120 {
01121 float a;
01122
01123 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
01124 {
01125 fprintf( ioQQQ, " HeCreate could not scan he_transprob\n" );
01126 puts( "[Stop in HeCreate]" );
01127 cdEXIT(EXIT_FAILURE);
01128 }
01129 ++chTemp;
01130
01131 sscanf( chTemp , "%e" , &a );
01132 TransProbs[nelem][i2][i1] = a;
01133
01134
01135 if( lgEOL )
01136 {
01137 fprintf( ioQQQ, " HeCreate detected insanity in he_transprob.dat.\n");
01138 puts( "[Stop in HeCreate]" );
01139 cdEXIT(EXIT_FAILURE);
01140 }
01141 }
01142 }
01143
01144
01145 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01146 {
01147 fprintf( ioQQQ, " HeCreate could not read last line of he_transprob.dat.\n");
01148 puts( "[Stop in HeCreate]" );
01149 cdEXIT(EXIT_FAILURE);
01150 }
01151 i = 1;
01152 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01153 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01154 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
01155 {
01156 fprintf( ioQQQ,
01157 " HeCreate: the version of he_transprob.dat is not the current version.\n" );
01158 fprintf( ioQQQ,
01159 " HeCreate: I expected to find the number %i %i and got %li %li instead.\n" ,
01160 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
01161 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01162 puts( "[Stop in HeCreate]" );
01163 cdEXIT(EXIT_FAILURE);
01164 }
01165
01166
01167 fclose( ioDATA );
01168 }
01169
01170 return;
01171 }
01172
01173
01174
01175
01176
01177
01178 static double CalculateRadialIntegral( long nelem, long ipLo, long ipHi )
01179 {
01180 double RadInt = 0.;
01181
01182 RI_ipHi = ipHi;
01183 RI_ipLo = ipLo;
01184
01185
01186 globalZ = nelem;
01187
01188 if( ( nelem == ipHELIUM) )
01189 {
01190
01191 RadInt = qg32( 0., 0.2, RadialIntegrand );
01192 RadInt += qg32( 0.2, 0.4, RadialIntegrand );
01193 RadInt += qg32( 0.4, 0.6, RadialIntegrand );
01194 RadInt += qg32( 0.8, 1.0, RadialIntegrand );
01195 RadInt += qg32( 1.0, 2.0, RadialIntegrand );
01196 RadInt += qg32( 2.0, 4.0, RadialIntegrand );
01197 RadInt += qg32( 4.0, 10., RadialIntegrand );
01198 }
01199 else
01200 fprintf(ioQQQ, "No can do in CalculateRadialIntegral, nelem %li, iplo %li, ipHi %li\n",
01201 nelem, ipLo, ipHi);
01202
01203
01204 RadInt *= 1.;
01205
01206 return RadInt;
01207 }
01208
01209 static double RadialIntegrand( double r12 )
01210 {
01211 double result;
01212
01213 RI_ipLev = RI_ipHi;
01214
01215 result = r12 * WaveFunction( r12 );
01216
01217 RI_ipLev = RI_ipLo;
01218
01219 result *= WaveFunction( r12 );
01220
01221 return result;
01222 }
01223
01224 static double WaveFunction( double r12 )
01225 {
01226 double value,PsiSlow2,alphaOfR2, PsiFastOneTwo;
01227 long Z = globalZ;
01228
01229
01230 double R1 = POW2((double)Z) * 5.3E-5;
01231 double R2 = r12 + R1;
01232 double a=0.,b=0.,mu=0.,nu=0.,Beta=0.;
01233 long l = iso.quant_desig[ipHE_LIKE][Z][RI_ipLev].l;
01234
01235
01236
01237
01238
01239
01240
01241
01242 if( l == 0 )
01243 {
01244 Beta = 0.411;
01245 PsiSlow2 = sqrt( POW3(Beta) / PI ) * exp( -1. * Beta * R2 );
01246
01247 a = -0.1;
01248 b = 0.14;
01249 mu = 1.34;
01250 nu = 2.86;
01251 }
01252 else if( l == 1 )
01253 {
01254 Beta = 0.5;
01255
01256 PsiSlow2 = 1.1547 * pow( Beta, 2.5 ) * R2 * exp( -1. * Beta * R2 );
01257
01258 a = -1;
01259 b = 0.2;
01260 mu = 0.1;
01261 nu = 4.0;
01262 }
01263 else
01264 {
01265
01266 PsiSlow2 = 0.;
01267 }
01268
01269 alphaOfR2 = 2 - (1. + a * pow( R2, mu ) )/ ( 1. + b * pow( R2, nu ) );
01270
01271
01272 PsiFastOneTwo = sqrt( POW3(alphaOfR2) / PI ) * exp( -1. * alphaOfR2 * R1 );
01273
01274 value = PsiFastOneTwo * PsiSlow2;
01275
01276 return value;
01277 }
01278
01279