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 "dense.h"
00017 #include "trace.h"
00018 #include "hydro_bauman.h"
00019 #include "iso.h"
00020 #include "helike.h"
00021 #include "helike_einsta.h"
00022 #include "hydroeinsta.h"
00023
00024
00025 static double ***TransProbs;
00026
00027
00028
00029 STATIC double ritoa( long li, long lf, long nelem, double k, double RI2 );
00030
00031
00032 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem );
00033
00034
00035
00036
00037
00038
00039
00040 static double vJint , zJint;
00041
00042
00043 STATIC double Jint( double theta )
00044 {
00045
00046 double
00047 d0 = ( 1.0 / PI ),
00048 d1 = vJint * theta,
00049 d2 = zJint * sin(theta),
00050 d3 = (d1 - d2),
00051 d4 = cos(d3),
00052 d5 = (d0 * d4);
00053
00054 return( d5 );
00055 }
00056
00057
00058 STATIC double AngerJ( double vv, double zz )
00059 {
00060 long int rep = 0, ddiv, divsor;
00061
00062 double y = 0.0;
00063
00064 DEBUG_ENTRY( "AngerJ()" );
00065
00066
00067
00068
00069 if( (fabs(vv)) - (int)(fabs(vv)) > 0.5 )
00070 ddiv = (int)(fabs(vv)) + 1;
00071 else
00072 ddiv = (int)(fabs(vv));
00073
00074 divsor = ((ddiv == 0) ? 1 : ddiv);
00075 vJint = vv;
00076 zJint = zz;
00077
00078 for( rep = 0; rep < divsor; rep++ )
00079 {
00080 double
00081 rl = (((double) rep)/((double) divsor)),
00082 ru = (((double) (rep+1))/((double) divsor)),
00083 x_low = (PI * rl),
00084 x_up = (PI * ru);
00085
00086 y += qg32( x_low, x_up, Jint );
00087 }
00088
00089 return( y );
00090 }
00091
00092
00093
00094
00095
00096
00097
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 double scqdri(
00147 double nstar, long int l,
00148 double npstar, long int lp,
00149 double iz
00150 )
00151 {
00152 double n_c = ((2.0 * nstar * npstar ) / ( nstar + npstar ));
00153
00154 double D_n = (nstar - npstar);
00155 double D_l = (double) ( l - lp );
00156 double lg = (double) ( (lp > l) ? lp : l);
00157
00158 double h = (lg/n_c);
00159 double g = h*h;
00160 double f = ( 1.0 - g );
00161 double e = (( f >= 0.0) ? sqrt( f ) : 0.0 );
00162
00163 double x = (e * D_n);
00164 double z = (-1.0 * x);
00165 double v1 = (D_n + 1.0);
00166 double v2 = (D_n - 1.0);
00167
00168 double d1,d2,d7,d8,d9,d34,d56,d6_1;
00169
00170 DEBUG_ENTRY( "scqdri()" );
00171
00172 if( iz == 0.0 )
00173 iz += 1.0;
00174
00175 if( D_n == 0.0 )
00176 {
00177 return( -1.0 );
00178 }
00179
00180 if( D_n < 0.0 )
00181 {
00182 return( -1.0 );
00183 }
00184
00185 if( f < 0.0 )
00186 {
00187
00188
00189
00190
00191 return( -1.0 );
00192 }
00193
00194 d1 = ( 1.0 / iz );
00195
00196 d2 = (n_c * n_c)/(2.0 * D_n);
00197
00198 d34 = (1.0 - ((D_l * lg)/n_c)) * AngerJ( v1, z );
00199
00200 d56 = (1.0 + ((D_l * lg)/n_c)) * AngerJ( v2, z );
00201
00202 d6_1 = PI * D_n;
00203
00204 d7 = (2./PI) * sin( d6_1 ) * (1.0 - e);
00205
00206 d8 = d1 * d2 * ( (d34) - (d56) + d7 );
00207
00208 d9 = d8 * d8;
00209
00210 ASSERT( D_n > 0.0 );
00211 ASSERT( l >= 0 );
00212 ASSERT( lp >= 0 );
00213 ASSERT( (l == lp + 1) || ( l == lp - 1) );
00214 ASSERT( n_c != 0.0 );
00215 ASSERT( f >= 0.0 );
00216 ASSERT( d9 > 0.0 );
00217
00218 return( d9 );
00219 }
00220
00221 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem )
00222 {
00223 double A;
00224
00225
00226
00227 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,
00228 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,
00229 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};
00230
00231
00232
00233
00234
00235 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,
00236 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,
00237 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};
00238
00239 DEBUG_ENTRY( "ForbiddenAuls()" );
00240
00241 int ipISO = ipHE_LIKE;
00242
00243 if( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
00244 {
00245 if( nelem == ipHELIUM )
00246 {
00247
00248
00249
00250
00251 double ForbiddenHe[5] = { 1.272E-4, 51.02, 1E-20, 177.58, 0.327 };
00252
00253 fixit();
00254 A = ForbiddenHe[ipHi - 1];
00255 iso_put_error(ipHE_LIKE,nelem,ipHe2s3S ,ipHe1s1S,IPRAD, 0.01f, 0.20f);
00256 iso_put_error(ipHE_LIKE,nelem,ipHe2s1S ,ipHe1s1S,IPRAD, 0.01f, 0.05f);
00257 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P0,ipHe1s1S,IPRAD, 0.00f, 0.00f);
00258 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P1,ipHe1s1S,IPRAD, 0.01f, 0.05f);
00259 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P2,ipHe1s1S,IPRAD, 0.01f, 0.01f);
00260 }
00261 else
00262 {
00263 switch ( (int)ipHi )
00264 {
00265 case 1:
00266
00267
00268 A = (3.9061E-7) * pow( (double)nelem+1., 10.419 ) + As2nuFrom3S[nelem-2];
00269 break;
00270 case 2:
00271 A = As2nuFrom1S[nelem-2];
00272 break;
00273 case 3:
00274 A = iso.SmallA;
00275 break;
00276 case 4:
00277
00278
00279 if( nelem <= ipARGON )
00280 {
00281 A = ( 11.431 * pow((double)nelem, 9.091) );
00282 }
00283 else
00284 {
00285
00286 A = ( 383.42 * pow((double)nelem, 7.8901) );
00287 }
00288 break;
00289 case 5:
00290
00291
00292
00293
00294 A = ( 0.11012 * pow((double)nelem, 7.6954) );
00295 break;
00296 default:
00297 TotalInsanity();
00298 }
00299 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00300 }
00301
00302 return A;
00303 }
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313 else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==3 &&
00314 L_(ipLo)==0 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
00315 {
00316 A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
00317 pow((double)N_(ipHi),2.877);
00318 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00319 }
00320
00321
00322 else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==1 &&
00323 L_(ipLo)==1 && S_(ipLo)==3 && N_(ipLo) < N_(ipHi) )
00324 {
00325 A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
00326
00327 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
00328 {
00329
00330
00331 A *= (2.*(ipLo-3)+1.0)/3.0;
00332 }
00333 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00334 }
00335
00336 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe3s3S ) )
00337 {
00338 double As_3TripS_to_2TripS[9] = {
00339 7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
00340 1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
00341
00342
00343
00344 if( nelem <= ipNEON )
00345 {
00346 A = As_3TripS_to_2TripS[nelem-1];
00347
00348
00349 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f,0.2f);
00350 }
00351 else
00352 {
00353
00354
00355 A= 7.22E-9*pow((double)nelem, 9.33);
00356 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f,0.3f);
00357 }
00358 }
00359
00360 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
00361 {
00362
00363 if( nelem == ipHELIUM )
00364 {
00365 A = 1.549;
00366 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00367 }
00368 else
00369 {
00370
00371
00372
00373 A= 0.1834*pow((double)nelem, 6.5735);
00374 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00375 }
00376 }
00377
00378 else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 )
00379 {
00380
00381 fixit();
00382
00387 A = 5.78E-12;
00388 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00389
00390 }
00391
00392 else if( nelem==ipHELIUM && ipHi==5 && ipLo==4 )
00393 {
00394 fixit();
00395
00400 A = 3.61E-15;
00401 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00402 }
00403
00404 else
00405 {
00406
00407 A = iso.SmallA;
00408 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00409 }
00410
00411 ASSERT( A > 0.);
00412
00413 return A;
00414 }
00415
00416
00417 double he_1trans(
00418
00419 long nelem , double Enerwn ,
00420
00421 double Eff_nupper, long lHi, long sHi, long jHi,
00422
00423 double Eff_nlower, long lLo, long sLo, long jLo )
00424
00425
00426 {
00427 int ipISO = ipHE_LIKE;
00428 double RI2, Aul;
00429 long nHi, nLo, ipHi, ipLo;
00430
00431 DEBUG_ENTRY( "he_1trans()" );
00432
00433 ASSERT(nelem > ipHYDROGEN);
00434
00435
00436
00437 nHi = (int)(Eff_nupper + 0.4);
00438 nLo = (int)(Eff_nlower + 0.4);
00439
00440
00441 ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
00442 ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );
00443
00444 ipHi = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nHi][lHi][sHi];
00445 if( (nHi==2) && (lHi==1) && (sHi==3) )
00446 {
00447 ASSERT( (jHi>=0) && (jHi<=2) );
00448 ipHi -= (2 - jHi);
00449 }
00450
00451 ipLo = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nLo][lLo][sLo];
00452 if( (nLo==2) && (lLo==1) && (sLo==3) )
00453 {
00454 ASSERT( (jLo>=0) && (jLo<=2) );
00455 ipLo -= (2 - jLo);
00456 }
00457
00458 ASSERT( StatesElemNEW[nelem][nelem-ipHE_LIKE][ipHi].n == nHi );
00459 if( nHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00460 {
00461 ASSERT( StatesElemNEW[nelem][nelem-ipHE_LIKE][ipHi].l == lHi );
00462 ASSERT( StatesElemNEW[nelem][nelem-ipHE_LIKE][ipHi].S == sHi );
00463 }
00464 ASSERT( StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLo].n == nLo );
00465 if( nLo <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00466 {
00467 ASSERT( StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLo].l == lLo );
00468 ASSERT( StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLo].S == sLo );
00469 }
00470
00471
00472 if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
00473 {
00474 Aul = -2.;
00475
00476
00477 if( nelem == ipHELIUM )
00478 {
00479
00480
00481 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
00482 {
00483
00484 ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ) );
00485 ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ) );
00486 ASSERT( ipHi > 2 );
00487
00488 Aul = TransProbs[nelem][ipHi][ipLo];
00489
00490 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0001f,0.002f);
00491 }
00492
00493 if( Aul < 0. )
00494 {
00495
00496 if( ipLo == ipHe1s1S )
00497 {
00498 ASSERT( (lHi == 1) && (sHi == 1) );
00499
00500
00501 if( nLo == 1 )
00502 Aul = (1.59208e10) / pow(Eff_nupper,3.0);
00503 ASSERT( Aul > 0.);
00504 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f);
00505 }
00506
00507
00508
00509 else if( lHi>=2 && lLo>=2 && nHi>nLo )
00510 {
00511 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
00512 ASSERT( Aul > 0.);
00513 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.006f,0.04f);
00514 }
00515
00516
00517
00518
00519
00520
00521 else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
00522 {
00523 int paramSet=0;
00524 double emisOscStr, x, a, b, c;
00525 double extrapol_Params[2][4][4][3] = {
00526
00527 {
00528 {
00529 { 0.8267396 , 1.4837624 , -0.4615955 },
00530 { 1.2738405 , 1.5841806 , -0.3022984 },
00531 { 1.6128996 , 1.6842538 , -0.2393057 },
00532 { 1.8855491 , 1.7709125 , -0.2115213 }},
00533 {
00534 { -1.4293664 , 2.3294080 , -0.0890470 },
00535 { -0.3608082 , 2.3337636 , -0.0712380 },
00536 { 0.3027974 , 2.3326252 , -0.0579008 },
00537 { 0.7841193 , 2.3320138 , -0.0497094 }},
00538 {
00539 { 1.1341403 , 3.1702435 , -0.2085843 },
00540 { 1.7915926 , 2.4942946 , -0.2266493 },
00541 { 2.1979400 , 2.2785377 , -0.1518743 },
00542 { 2.5018229 , 2.1925720 , -0.1081966 }},
00543 {
00544 { 0.0000000 , 0.0000000 , 0.0000000 },
00545 { -2.6737396 , 2.9379143 , -0.3805367 },
00546 { -1.4380124 , 2.7756396 , -0.2754625 },
00547 { -0.6630196 , 2.6887253 , -0.2216493 }},
00548 },
00549
00550 {
00551 {
00552 { 0.3075287 , 0.9087130 , -1.0387207 },
00553 { 0.687069 , 1.1485864 , -0.6627317 },
00554 { 0.9776064 , 1.3382024 , -0.5331906 },
00555 { 1.2107725 , 1.4943721 , -0.4779232 }},
00556 {
00557 { -1.3659605 , 2.3262253 , -0.0306439 },
00558 { -0.2899490 , 2.3279391 , -0.0298695 },
00559 { 0.3678878 , 2.3266603 , -0.0240021 },
00560 { 0.8427457 , 2.3249540 , -0.0194091 }},
00561 {
00562 { 1.3108281 , 2.8446367 , -0.1649923 },
00563 { 1.8437692 , 2.2399326 , -0.2583398 },
00564 { 2.1820792 , 2.0693762 , -0.1864091 },
00565 { 2.4414052 , 2.0168255 , -0.1426083 }},
00566 {
00567 { 0.0000000 , 0.0000000 , 0.0000000 },
00568 { -1.9219877 , 2.7689624 , -0.2536072 },
00569 { -0.7818065 , 2.6595150 , -0.1895313 },
00570 { -0.0665624 , 2.5955623 , -0.1522616 }},
00571 }
00572 };
00573
00574 if( lLo==0 )
00575 {
00576 paramSet = 0;
00577 }
00578 else if( lLo==1 && lHi==0 )
00579 {
00580 paramSet = 1;
00581 }
00582 else if( lLo==1 && lHi==2 )
00583 {
00584 paramSet = 2;
00585 }
00586 else if( lLo==2 )
00587 {
00588 paramSet = 3;
00589 ASSERT( lHi==1 );
00590 }
00591
00592 ASSERT( (int)((sHi-1)/2) == 0 || (int)((sHi-1)/2) == 1 );
00593 a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
00594 b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
00595 c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
00596 ASSERT( Enerwn > 0. );
00597 x = log( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]*RYD_INF/Enerwn );
00598
00599 emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
00600 (2.*lLo+1)/(2.*lHi+1);
00601
00602 Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;
00603
00604 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
00605 {
00606 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00607 }
00608
00609 ASSERT( Aul > 0. );
00610 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f);
00611 }
00612 else
00613 {
00614
00615 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
00616 ASSERT( RI2 > 0. );
00617
00618 Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
00619
00620
00621 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
00622 {
00623 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00624 }
00625
00626 ASSERT( Aul > 0. );
00627 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.07f);
00628 }
00629 }
00630 }
00631
00632
00633 else
00634 {
00635
00636
00637
00638 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
00639 {
00640
00641 ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) );
00642 ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) );
00643 ASSERT( ipHi > 2 );
00644
00645 Aul = TransProbs[nelem][ipHi][ipLo];
00646 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00647 }
00648
00649 if( Aul < 0. )
00650 {
00651
00652 if( nLo==nHi && lHi<=2 && lLo<=2 )
00653 {
00654
00655
00656 if( ipLo == ipHe2s3S )
00657 {
00658 if(ipHi == ipHe2p3P0)
00659 Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
00660 else if(ipHi == ipHe2p3P1)
00661 Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
00662 else if(ipHi == ipHe2p3P2)
00663 Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
00664 else
00665 {
00666 fprintf(ioQQQ," PROBLEM Impossible value for ipHi in he_1trans.\n");
00667 TotalInsanity();
00668 }
00669 }
00670
00671
00672 else if( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) )
00673 {
00674 Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
00675 }
00676
00677 else
00678 {
00679
00680 ASSERT( nLo > 2 );
00681
00682
00683 if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
00684 {
00685 Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
00686 }
00687
00688 else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
00689 {
00690 Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
00691 }
00692
00693 else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
00694 {
00695 Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077);
00696 }
00697 else
00698 {
00699 ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
00700 Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
00701 if( (lHi >2) || (lLo > 2) )
00702 Aul *= (lHi/2.);
00703 if(lLo > 2)
00704 Aul *= (1./9.);
00705 }
00706 }
00707 ASSERT( Aul > 0.);
00708 }
00709
00710
00711 else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
00712 {
00713 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
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 == 1) );
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 == 1) );
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 == 3) );
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 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00779 }
00780 }
00781 }
00782
00783
00784
00785
00786
00787
00788
00789
00790 else
00791 {
00792 ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
00793 Aul = ForbiddenAuls(ipHi, ipLo, nelem);
00794 ASSERT( Aul > 0. );
00795 }
00796
00797 Aul = MAX2( Aul, iso.SmallA );
00798 ASSERT( Aul >= iso.SmallA );
00799
00800
00801
00802 if( Enerwn < 0. && Aul > iso.SmallA )
00803 {
00804 fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
00805 }
00806
00807 return Aul;
00808 }
00809
00810 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
00811 {
00812
00813
00814
00815
00816
00817
00818
00819 long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
00820 double Ass, Att, Ast, Ats;
00821 double SinHi, SinLo, CosHi, CosLo;
00822 double HiMixingAngle, LoMixingAngle , error;
00823 double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
00824
00825 DEBUG_ENTRY( "DoFSMixing()" );
00826
00827 nHi = StatesElemNEW[nelem][nelem-ipHE_LIKE][ipHiSing].n;
00828 lHi = StatesElemNEW[nelem][nelem-ipHE_LIKE][ipHiSing].l;
00829 sHi = StatesElemNEW[nelem][nelem-ipHE_LIKE][ipHiSing].S;
00830 nLo = StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLoSing].n;
00831 lLo = StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLoSing].l;
00832 sLo = StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLoSing].S;
00833
00834 if( ( sHi == 3 || sLo == 3 ) ||
00835 ( abs(lHi - lLo) != 1 ) ||
00836 ( nLo < 2 ) ||
00837 ( lHi <= 1 || lLo <= 1 ) ||
00838 ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
00839 ( nHi > nLo && lHi != 1 && lLo != 1 ) )
00840 {
00841 return;
00842 }
00843
00844 ASSERT( lHi > 0 );
00845
00846
00847 ipHiTrip = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nHi][lHi][3];
00848 ipLoTrip = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nLo][lLo][3];
00849
00850 if( lHi == 2 )
00851 {
00852 HiMixingAngle = 0.01;
00853 }
00854 else if( lHi == 3 )
00855 {
00856 HiMixingAngle = 0.5;
00857 }
00858 else
00859 {
00860 HiMixingAngle = PI/4.;
00861 }
00862
00863 if( lLo == 2 )
00864 {
00865 LoMixingAngle = 0.01;
00866 }
00867 else if( lLo == 3 )
00868 {
00869 LoMixingAngle = 0.5;
00870 }
00871 else
00872 {
00873 LoMixingAngle = PI/4.;
00874 }
00875
00876
00877 ASSERT( ipHiTrip > ipLoTrip );
00878 ASSERT( ipHiTrip > ipLoSing );
00879 ASSERT( ipHiSing > ipLoTrip );
00880 ASSERT( ipHiSing > ipLoSing );
00881
00882 SinHi = sin( HiMixingAngle );
00883 SinLo = sin( LoMixingAngle );
00884 CosHi = cos( HiMixingAngle );
00885 CosLo = cos( LoMixingAngle );
00886
00887 Kss = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].EnergyWN;
00888 Ktt = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].EnergyWN;
00889 Kst = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].EnergyWN;
00890 Kts = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].EnergyWN;
00891
00892 fss = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul/TRANS_PROB_CONST/POW2(Kss)/
00893 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Lo->g*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Hi->g;
00894
00895 ftt = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul/TRANS_PROB_CONST/POW2(Ktt)/
00896 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Lo->g*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Hi->g;
00897
00898 temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
00899 fssNew = Kss*POW2( temp );
00900 temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
00901 fttNew = Ktt*POW2( temp );
00902 temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
00903 fstNew = Kst*POW2( temp );
00904 temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
00905 ftsNew = Kts*POW2( temp );
00906
00907 Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Lo->g/
00908 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Hi->g;
00909
00910 Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Lo->g/
00911 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Hi->g;
00912
00913 Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Lo->g/
00914 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Hi->g;
00915
00916 Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Lo->g/
00917 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Hi->g;
00918
00919 error = fabs( ( Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul+
00920 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul )/
00921 (Ass+Ast+Ats+Att) - 1.f );
00922
00923 if( error > 0.001 )
00924 {
00925 fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
00926 ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
00927 Ass/Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul,
00928 Att/Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul,
00929 Ast/Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Emis->Aul,
00930 Ats/Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Emis->Aul );
00931 }
00932
00933 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul = (realnum)Ass;
00934 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul = (realnum)Att;
00935 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Emis->Aul = (realnum)Ast;
00936 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Emis->Aul = (realnum)Ats;
00937
00938 return;
00939 }
00940
00941
00942
00943 STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
00944 {
00945
00946
00947
00948
00949
00950
00951
00952
00953 long lg;
00954 double fmean,mu,EinsteinA,w,RI2_cm;
00955
00956 DEBUG_ENTRY( "ritoa()" );
00957
00958 mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT));
00959
00960 w = 2. * PI * k * SPEEDLIGHT;
00961
00962 RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM;
00963
00964 lg = (lf > li) ? lf : li;
00965
00966 fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));
00967
00968 EinsteinA = TRANS_PROB_CONST*k*k*fmean;
00969
00970
00971
00972 return EinsteinA;
00973 }
00974
00975 realnum helike_transprob( long nelem, long ipHi, long ipLo )
00976 {
00977 double Aul, Aul1;
00978 double Enerwn, n_eff_hi, n_eff_lo;
00979 long ipISO = ipHE_LIKE;
00980
00981 double z4 = POW2((double)nelem);
00982 z4 *= z4;
00983
00984 DEBUG_ENTRY( "helike_transprob()" );
00985
00986 Enerwn = Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN;
00987 n_eff_hi = N_(ipHi) - helike_quantum_defect(nelem,ipHi);
00988 n_eff_lo = N_(ipLo) - helike_quantum_defect(nelem,ipLo);
00989
00990 if( ipHi >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] )
00991 {
00992 if( ipLo >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] )
00993 {
00994
00995 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
00996 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f,0.001f);
00997
00998 ASSERT( Aul > 0.);
00999 }
01000 else
01001 {
01002
01003
01004 Aul = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)+1,
01005 S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
01006
01007 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][0] = (realnum)Aul;
01008
01009 Aul *= (2.*L_(ipLo)+3.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
01010
01011 if( L_(ipLo) != 0 )
01012 {
01013
01014 Aul1 = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)-1,
01015 S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
01016
01017 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = (realnum)Aul1;
01018
01019
01020 Aul += Aul1*(2.*L_(ipLo)-1.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
01021 }
01022 else
01023 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = 0.f;
01024
01025 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.01f);
01026 ASSERT( Aul > 0.);
01027 }
01028 }
01029 else
01030 {
01031
01032 if( Enerwn < 0. )
01033 {
01034 Aul = he_1trans( nelem, -1.*Enerwn, n_eff_lo,
01035 L_(ipLo), S_(ipLo), ipLo-3, n_eff_hi, L_(ipHi), S_(ipHi), ipHi-3);
01036 }
01037 else
01038 {
01039 Aul = he_1trans( nelem, Enerwn, n_eff_hi,
01040 L_(ipHi), S_(ipHi), ipHi-3, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
01041 }
01042 }
01043
01044 return (realnum)Aul;
01045 }
01046
01047
01048 void HelikeTransProbSetup( void )
01049 {
01050
01051 const int chLine_LENGTH = 1000;
01052 char chLine[chLine_LENGTH];
01053
01054 FILE *ioDATA;
01055 bool lgEOL;
01056
01057 long nelem, ipLo, ipHi, i, i1, i2, i3;
01058
01059 DEBUG_ENTRY( "HelikeTransProbSetup()" );
01060
01061 TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
01062
01063 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
01064 {
01065
01066 TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) );
01067
01068 for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo )
01069 {
01070 TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX );
01071 }
01072 }
01073
01074
01075
01076 if( trace.lgTrace )
01077 fprintf( ioQQQ," HelikeTransProbSetup opening he_transprob.dat:");
01078
01079 ioDATA = open_data( "he_transprob.dat", "r" );
01080
01081
01082 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01083 {
01084 fprintf( ioQQQ, " HelikeTransProbSetup could not read first line of he_transprob.dat.\n");
01085 cdEXIT(EXIT_FAILURE);
01086 }
01087 i = 1;
01088 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01089 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01090 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
01091 {
01092 fprintf( ioQQQ,
01093 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
01094 fprintf( ioQQQ,
01095 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
01096 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
01097 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01098 cdEXIT(EXIT_FAILURE);
01099 }
01100
01101
01102 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
01103 {
01104 for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
01105 {
01106 for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
01107 {
01108 TransProbs[nelem][ipHi][ipLo] = -1.;
01109 }
01110 }
01111 }
01112
01113 for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
01114 {
01115 char *chTemp;
01116
01117
01118 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01119 BadRead();
01120
01121 while( chLine[0]=='#' )
01122 {
01123 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01124 BadRead();
01125 }
01126
01127 i3 = 1;
01128 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
01129 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
01130
01131 if( i1<0 || i2<=i1 )
01132 {
01133 fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
01134 cdEXIT(EXIT_FAILURE);
01135 }
01136
01137 chTemp = chLine;
01138
01139
01140 for( i=0; i<1; ++i )
01141 {
01142 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
01143 {
01144 fprintf( ioQQQ, " HelikeTransProbSetup could not init he_transprob\n" );
01145 cdEXIT(EXIT_FAILURE);
01146 }
01147 ++chTemp;
01148 }
01149
01150
01151 for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
01152 {
01153 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
01154 {
01155 fprintf( ioQQQ, " HelikeTransProbSetup could not scan he_transprob\n" );
01156 cdEXIT(EXIT_FAILURE);
01157 }
01158 ++chTemp;
01159
01160 sscanf( chTemp , "%le" , &TransProbs[nelem][i2][i1] );
01161
01163 if( lgEOL )
01164 {
01165 fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
01166 cdEXIT(EXIT_FAILURE);
01167 }
01168 }
01169 }
01170
01171
01172 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01173 {
01174 fprintf( ioQQQ, " HelikeTransProbSetup could not read last line of he_transprob.dat.\n");
01175 cdEXIT(EXIT_FAILURE);
01176 }
01177 i = 1;
01178 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01179 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01180 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
01181 {
01182 fprintf( ioQQQ,
01183 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
01184 fprintf( ioQQQ,
01185 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
01186 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
01187 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01188 cdEXIT(EXIT_FAILURE);
01189 }
01190
01191
01192 fclose( ioDATA );
01193 return;
01194 }