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 DEBUG_ENTRY( "ForbiddenAuls()" );
00226
00227 int ipISO = ipHE_LIKE;
00228
00229 if( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
00230 {
00231 if( nelem == ipHELIUM )
00232 {
00233
00234
00235
00236 double ForbiddenHe[5] = { 1.272E-4, 1E-20, 1E-20, 177.58, 0.327 };
00237
00238 A = ForbiddenHe[ipHi - 1];
00239 iso_put_error(ipHE_LIKE,nelem,ipHe2s3S ,ipHe1s1S,IPRAD, 0.01f, 0.20f);
00240 iso_put_error(ipHE_LIKE,nelem,ipHe2s1S ,ipHe1s1S,IPRAD, 0.01f, 0.05f);
00241 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P0,ipHe1s1S,IPRAD, 0.00f, 0.00f);
00242 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P1,ipHe1s1S,IPRAD, 0.01f, 0.05f);
00243 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P2,ipHe1s1S,IPRAD, 0.01f, 0.01f);
00244 }
00245 else
00246 {
00247 switch ( (int)ipHi )
00248 {
00249 case 1:
00250
00251
00252
00253 A = (3.9061E-7) * pow( (double)nelem+1., 10.419 );
00254 break;
00255 case 2:
00256 A = iso_ctrl.SmallA;
00257 break;
00258 case 3:
00259 A = iso_ctrl.SmallA;
00260 break;
00261 case 4:
00262
00263
00264 if( nelem <= ipARGON )
00265 {
00266 A = ( 11.431 * pow((double)nelem, 9.091) );
00267 }
00268 else
00269 {
00270
00271 A = ( 383.42 * pow((double)nelem, 7.8901) );
00272 }
00273 break;
00274 case 5:
00275
00276
00277
00278
00279 A = ( 0.11012 * pow((double)nelem, 7.6954) );
00280 break;
00281 default:
00282 TotalInsanity();
00283 }
00284 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00285 }
00286
00287 return A;
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==3 &&
00299 L_(ipLo)==0 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
00300 {
00301 A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
00302 pow((double)N_(ipHi),2.877);
00303 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00304 }
00305
00306
00307 else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==1 &&
00308 L_(ipLo)==1 && S_(ipLo)==3 && N_(ipLo) < N_(ipHi) )
00309 {
00310 A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
00311
00312 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
00313 {
00314
00315
00316 A *= (2.*(ipLo-3)+1.0)/3.0;
00317 }
00318 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00319 }
00320
00321 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe3s3S ) )
00322 {
00323 double As_3TripS_to_2TripS[9] = {
00324 7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
00325 1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
00326
00327
00328
00329 if( nelem <= ipNEON )
00330 {
00331 A = As_3TripS_to_2TripS[nelem-1];
00332
00333
00334 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f,0.2f);
00335 }
00336 else
00337 {
00338
00339
00340 A= 7.22E-9*pow((double)nelem, 9.33);
00341 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f,0.3f);
00342 }
00343 }
00344
00345 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
00346 {
00347
00348 if( nelem == ipHELIUM )
00349 {
00350 A = 1.549;
00351 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00352 }
00353 else
00354 {
00355
00356
00357
00358 A= 0.1834*pow((double)nelem, 6.5735);
00359 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00360 }
00361 }
00362
00363 else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 )
00364 {
00365
00366 fixit();
00367
00372 A = 5.78E-12;
00373 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00374
00375 }
00376
00377 else if( nelem==ipHELIUM && ipHi==5 && ipLo==4 )
00378 {
00379 fixit();
00380
00385 A = 3.61E-15;
00386 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00387 }
00388
00389 else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe1s1S )
00390 {
00391
00392 A = 44.326 * (1./3.) + 0.1146547 * (5./9.);
00393 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00394 }
00395
00396 else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe2s1S )
00397 {
00398
00399 A = 1.459495 * (1./3.) + 3.6558e-5 * (5./9.);
00400 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00401 }
00402
00403 else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe3s1S )
00404 {
00405
00406 A = 2.2297e-3 * (1./3.);
00407 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00408 }
00409
00410 else if( nelem==ipHELIUM && ipHi==ipHe3p1P && ipLo==ipHe2s3S )
00411 {
00412
00413 A = 0.1815436;
00414 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00415 }
00416
00417 else if( nelem==ipHELIUM && ipHi==ipHe3p1P && ipLo==ipHe3s3S )
00418 {
00419
00420 A = 0.14895852;
00421 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00422 }
00423
00424 else if( ( ipLo==0 && L_(ipHi)==2 && S_(ipHi)==1 ) ||
00425 ( ipLo==ipHe2s3S && L_(ipHi)==2 && S_(ipHi)==3 ) ||
00426 ( ipLo==ipHe2s1S && L_(ipHi)==2 && S_(ipHi)==1 ) ||
00427 ( N_(ipLo)==2 && L_(ipLo)==1 && S_(ipLo)==3 && N_(ipHi)>=3 && L_(ipHi)==1 && S_(ipHi)==3 ) ||
00428 ( ipLo==ipHe2p1P && L_(ipHi)==1 && S_(ipHi)==1 ) )
00429 {
00430
00431
00432 double f_params[5][4][3] = {
00433 {
00434 {9.360591E-007, -3.1937E-006, 3.5186E-006},
00435 {4.631435E-007, -1.4973E-006, 1.4848E-006},
00436 {2.493912E-007, -7.8658E-007, 7.3994E-007},
00437 {1.476742E-007, -4.5953E-007, 4.1932E-007}},
00438 {
00439 {1.646733E-006, -2.0028E-006, -1.3552E-006},
00440 {9.120593E-008, 3.1301E-007, -3.2190E-007},
00441 {1.360965E-008, 1.1467E-007, 8.6977E-008},
00442 {3.199421E-009, 4.5485E-008, 1.1016E-007}},
00443 {
00444 {1.646733E-006, -2.9720E-006, 1.5367E-006},
00445 {9.120593E-008, -3.9037E-008, 3.9156E-008},
00446 {1.360965E-008, 1.4671E-008, 1.5626E-008},
00447 {3.199421E-009, 8.9806E-009, 1.2436E-008}},
00448 {
00449 {1.543812E-007, -2.8220E-007, 3.0261E-008},
00450 {3.648237E-008, -6.6824E-008, 4.5879E-009},
00451 {1.488556E-008, -2.7304E-008, 1.6628E-009},
00452 {7.678610E-009, -1.4112E-008, 6.8453E-010}},
00453 {
00454 {1.543812E-007, -2.8546E-007, 4.6605E-008},
00455 {3.648237E-008, -6.8422E-008, 1.7038E-008},
00456 {1.488556E-008, -2.8170E-008, 8.5012E-009},
00457 {7.678610E-009, -1.4578E-008, 4.6686E-009}}
00458 };
00459
00460 ASSERT( ipLo <= 6 );
00461 long index_lo;
00462 if( ipLo==ipHe2p1P )
00463 index_lo = 4;
00464 else if( ipLo==ipHe2p3P1 || ipLo==ipHe2p3P2 )
00465 index_lo = 3;
00466 else
00467 index_lo = ipLo;
00468
00469 ASSERT( N_(ipHi) >= 3 );
00470 long index_hi = MIN2( N_(ipHi), 6 ) - 3;
00471 double f_lu = POW2(nelem+1) * (
00472 f_params[index_lo][index_hi][0] +
00473 f_params[index_lo][index_hi][1]/(nelem+1) +
00474 f_params[index_lo][index_hi][2]/POW2(nelem+1) );
00475
00476
00477 if( N_(ipHi) > 6 )
00478 f_lu *= pow( 6. / N_(ipHi), 3. );
00479
00480 double gLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].g();
00481 double gHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].g();
00482 double eWN = iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN();
00483
00484 A = eina( gLo * f_lu, eWN, gHi );
00485 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00486 }
00487
00488 else
00489 {
00490
00491 A = iso_ctrl.SmallA;
00492 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00493 }
00494
00495 ASSERT( A > 0.);
00496
00497 return A;
00498 }
00499
00500
00501 double he_1trans(
00502
00503 long nelem , double Enerwn ,
00504
00505 double Eff_nupper, long lHi, long sHi, long jHi,
00506
00507 double Eff_nlower, long lLo, long sLo, long jLo )
00508
00509
00510 {
00511 int ipISO = ipHE_LIKE;
00512 double RI2, Aul;
00513 long nHi, nLo, ipHi, ipLo;
00514
00515 DEBUG_ENTRY( "he_1trans()" );
00516
00517 ASSERT(nelem > ipHYDROGEN);
00518
00519
00520
00521 nHi = (int)(Eff_nupper + 0.4);
00522 nLo = (int)(Eff_nlower + 0.4);
00523
00524
00525 ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
00526 ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );
00527
00528 ipHi = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][sHi];
00529 if( (nHi==2) && (lHi==1) && (sHi==3) )
00530 {
00531 ASSERT( (jHi>=0) && (jHi<=2) );
00532 ipHi -= (2 - jHi);
00533 }
00534
00535 ipLo = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][sLo];
00536 if( (nLo==2) && (lLo==1) && (sLo==3) )
00537 {
00538 ASSERT( (jLo>=0) && (jLo<=2) );
00539 ipLo -= (2 - jLo);
00540 }
00541
00542 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == nHi );
00543 if( nHi <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
00544 {
00545 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].l() == lHi );
00546 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == sHi );
00547 }
00548 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].n() == nLo );
00549 if( nLo <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
00550 {
00551 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() == lLo );
00552 ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].S() == sLo );
00553 }
00554
00555
00556 if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
00557 {
00558 Aul = -2.;
00559
00560
00561 if( nelem == ipHELIUM )
00562 {
00563
00564
00565 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max )
00566 {
00567
00568 ASSERT( ipHi < ( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max ) );
00569 ASSERT( ipLo < ( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max ) );
00570 ASSERT( ipHi > 2 );
00571
00572 Aul = TransProbs[nelem][ipHi][ipLo];
00573
00574 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0001f,0.002f);
00575 }
00576
00577 if( Aul < 0. )
00578 {
00579
00580 if( ipLo == ipHe1s1S )
00581 {
00582 ASSERT( (lHi == 1) && (sHi == 1) );
00583
00584
00585 if( nLo == 1 )
00586 Aul = (1.59208e10) / pow(Eff_nupper,3.0);
00587 ASSERT( Aul > 0.);
00588 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f);
00589 }
00590
00591
00592
00593 else if( lHi>=2 && lLo>=2 && nHi>nLo )
00594 {
00595 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
00596 ASSERT( Aul > 0.);
00597 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.006f,0.04f);
00598 }
00599
00600
00601
00602
00603
00604
00605 else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
00606 {
00607 int paramSet=0;
00608 double emisOscStr, x, a, b, c;
00609 double extrapol_Params[2][4][4][3] = {
00610
00611 {
00612 {
00613 { 0.8267396 , 1.4837624 , -0.4615955 },
00614 { 1.2738405 , 1.5841806 , -0.3022984 },
00615 { 1.6128996 , 1.6842538 , -0.2393057 },
00616 { 1.8855491 , 1.7709125 , -0.2115213 }},
00617 {
00618 { -1.4293664 , 2.3294080 , -0.0890470 },
00619 { -0.3608082 , 2.3337636 , -0.0712380 },
00620 { 0.3027974 , 2.3326252 , -0.0579008 },
00621 { 0.7841193 , 2.3320138 , -0.0497094 }},
00622 {
00623 { 1.1341403 , 3.1702435 , -0.2085843 },
00624 { 1.7915926 , 2.4942946 , -0.2266493 },
00625 { 2.1979400 , 2.2785377 , -0.1518743 },
00626 { 2.5018229 , 2.1925720 , -0.1081966 }},
00627 {
00628 { 0.0000000 , 0.0000000 , 0.0000000 },
00629 { -2.6737396 , 2.9379143 , -0.3805367 },
00630 { -1.4380124 , 2.7756396 , -0.2754625 },
00631 { -0.6630196 , 2.6887253 , -0.2216493 }},
00632 },
00633
00634 {
00635 {
00636 { 0.3075287 , 0.9087130 , -1.0387207 },
00637 { 0.687069 , 1.1485864 , -0.6627317 },
00638 { 0.9776064 , 1.3382024 , -0.5331906 },
00639 { 1.2107725 , 1.4943721 , -0.4779232 }},
00640 {
00641 { -1.3659605 , 2.3262253 , -0.0306439 },
00642 { -0.2899490 , 2.3279391 , -0.0298695 },
00643 { 0.3678878 , 2.3266603 , -0.0240021 },
00644 { 0.8427457 , 2.3249540 , -0.0194091 }},
00645 {
00646 { 1.3108281 , 2.8446367 , -0.1649923 },
00647 { 1.8437692 , 2.2399326 , -0.2583398 },
00648 { 2.1820792 , 2.0693762 , -0.1864091 },
00649 { 2.4414052 , 2.0168255 , -0.1426083 }},
00650 {
00651 { 0.0000000 , 0.0000000 , 0.0000000 },
00652 { -1.9219877 , 2.7689624 , -0.2536072 },
00653 { -0.7818065 , 2.6595150 , -0.1895313 },
00654 { -0.0665624 , 2.5955623 , -0.1522616 }},
00655 }
00656 };
00657
00658 if( lLo==0 )
00659 {
00660 paramSet = 0;
00661 }
00662 else if( lLo==1 && lHi==0 )
00663 {
00664 paramSet = 1;
00665 }
00666 else if( lLo==1 && lHi==2 )
00667 {
00668 paramSet = 2;
00669 }
00670 else if( lLo==2 )
00671 {
00672 paramSet = 3;
00673 ASSERT( lHi==1 );
00674 }
00675
00676 ASSERT( (int)((sHi-1)/2) == 0 || (int)((sHi-1)/2) == 1 );
00677 a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
00678 b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
00679 c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
00680 ASSERT( Enerwn > 0. );
00681 x = log( iso_sp[ipHE_LIKE][nelem].fb[ipLo].xIsoLevNIonRyd*RYD_INF/Enerwn );
00682
00683 emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
00684 (2.*lLo+1)/(2.*lHi+1);
00685
00686 Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;
00687
00688 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
00689 {
00690 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00691 }
00692
00693 ASSERT( Aul > 0. );
00694 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f);
00695 }
00696 else
00697 {
00698
00699 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
00700 ASSERT( RI2 > 0. );
00701
00702 Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
00703
00704
00705 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
00706 {
00707 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00708 }
00709
00710 ASSERT( Aul > 0. );
00711 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.07f);
00712 }
00713 }
00714 }
00715
00716
00717 else
00718 {
00719
00720
00721
00722 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
00723 {
00724
00725 ASSERT( ipHi < ( iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) );
00726 ASSERT( ipLo < ( iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) );
00727 ASSERT( ipHi > 2 );
00728
00729 Aul = TransProbs[nelem][ipHi][ipLo];
00730 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00731 }
00732
00733 if( Aul < 0. )
00734 {
00735
00736 if( nLo==nHi && lHi<=2 && lLo<=2 )
00737 {
00738
00739
00740 if( ipLo == ipHe2s3S )
00741 {
00742 if(ipHi == ipHe2p3P0)
00743 Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
00744 else if(ipHi == ipHe2p3P1)
00745 Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
00746 else if(ipHi == ipHe2p3P2)
00747 Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
00748 else
00749 {
00750 fprintf(ioQQQ," PROBLEM Impossible value for ipHi in he_1trans.\n");
00751 TotalInsanity();
00752 }
00753 }
00754
00755
00756 else if( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) )
00757 {
00758 Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
00759 }
00760
00761 else
00762 {
00763
00764 ASSERT( nLo > 2 );
00765
00766
00767 if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
00768 {
00769 Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
00770 }
00771
00772 else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
00773 {
00774 Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
00775 }
00776
00777 else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
00778 {
00779 Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077);
00780 }
00781 else
00782 {
00783 ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
00784 Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
00785 if( (lHi >2) || (lLo > 2) )
00786 Aul *= (lHi/2.);
00787 if(lLo > 2)
00788 Aul *= (1./9.);
00789 }
00790 }
00791 ASSERT( Aul > 0.);
00792 }
00793
00794
00795 else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
00796 {
00797 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
00798 ASSERT( Aul > 0.);
00799 }
00800
00801
00802
00803
00804
00805 else if( ( ipLo == 0 ) || ( ipLo == ipHe2s1S ) || ( ipLo == ipHe2s3S )
00806 || ( ipLo == ipHe3s3S ) || ( nLo==4 && lLo==0 && sLo==3 ) )
00807 {
00808
00809 if( ipLo == 0 )
00810 {
00811 ASSERT( (lHi == 1) && (sHi == 1) );
00812
00813
00814
00815
00816
00817
00818 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
00819 }
00820
00821
00822 else if( ipLo == ipHe2s1S )
00823 {
00824 ASSERT( (lHi == 1) && (sHi == 1) );
00825
00826
00827 Aul = 5.0e8 * pow((double)nelem,4.) / pow((double)nHi, 2.889);
00828 }
00829
00830
00831 else
00832 {
00833 ASSERT( (lHi == 1) && (sHi == 3) );
00834
00835
00836 if( nLo == 2 )
00837 Aul = 1.5 * 3.405E8 * pow((double)nelem,4.) / pow((double)nHi, 2.883);
00838 else if( nLo == 3 )
00839 Aul = 2.5 * 4.613E7 * pow((double)nelem,4.) / pow((double)nHi, 2.672);
00840 else
00841 Aul = 3.0 * 1.436E7 * pow((double)nelem,4.) / pow((double)nHi, 2.617);
00842 }
00843
00844 ASSERT( Aul > 0.);
00845 }
00846
00847
00848 else
00849 {
00850
00851 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem));
00852
00853 Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2);
00854
00855
00856 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) && (Aul > iso_ctrl.SmallA) )
00857 {
00858 Aul *= (2.*(ipLo-3)+1.0)/9.0;
00859 }
00860
00861 }
00862 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
00863 }
00864 }
00865 }
00866
00867
00868
00869
00870
00871
00872
00873
00874 else
00875 {
00876 ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
00877 Aul = ForbiddenAuls(ipHi, ipLo, nelem);
00878 ASSERT( Aul > 0. );
00879 }
00880
00881 Aul = MAX2( Aul, iso_ctrl.SmallA );
00882 ASSERT( Aul >= iso_ctrl.SmallA );
00883
00884
00885
00886 if( Enerwn < 0. && Aul > iso_ctrl.SmallA )
00887 {
00888 fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
00889 }
00890
00891 return Aul;
00892 }
00893
00894 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
00895 {
00896
00897
00898
00899
00900
00901
00902
00903 long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
00904 double Ass, Att, Ast, Ats;
00905 double SinHi, SinLo, CosHi, CosLo;
00906 double HiMixingAngle, LoMixingAngle , error;
00907 double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
00908
00909 DEBUG_ENTRY( "DoFSMixing()" );
00910
00911 nHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].n();
00912 lHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].l();
00913 sHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].S();
00914 nLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].n();
00915 lLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].l();
00916 sLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].S();
00917
00918 if( ( sHi == 3 || sLo == 3 ) ||
00919 ( abs(lHi - lLo) != 1 ) ||
00920 ( nLo < 2 ) ||
00921 ( lHi <= 1 || lLo <= 1 ) ||
00922 ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
00923 ( nHi > nLo && lHi != 1 && lLo != 1 ) )
00924 {
00925 return;
00926 }
00927
00928 ASSERT( lHi > 0 );
00929
00930
00931 ipHiTrip = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][3];
00932 ipLoTrip = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][3];
00933
00934 if( lHi == 2 )
00935 {
00936 HiMixingAngle = 0.01;
00937 }
00938 else if( lHi == 3 )
00939 {
00940 HiMixingAngle = 0.5;
00941 }
00942 else
00943 {
00944 HiMixingAngle = PI/4.;
00945 }
00946
00947 if( lLo == 2 )
00948 {
00949 LoMixingAngle = 0.01;
00950 }
00951 else if( lLo == 3 )
00952 {
00953 LoMixingAngle = 0.5;
00954 }
00955 else
00956 {
00957 LoMixingAngle = PI/4.;
00958 }
00959
00960
00961 ASSERT( ipHiTrip > ipLoTrip );
00962 ASSERT( ipHiTrip > ipLoSing );
00963 ASSERT( ipHiSing > ipLoTrip );
00964 ASSERT( ipHiSing > ipLoSing );
00965
00966 SinHi = sin( HiMixingAngle );
00967 SinLo = sin( LoMixingAngle );
00968 CosHi = cos( HiMixingAngle );
00969 CosLo = cos( LoMixingAngle );
00970
00971 Kss = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).EnergyWN();
00972 Ktt = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).EnergyWN();
00973 Kst = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).EnergyWN();
00974 Kts = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).EnergyWN();
00975
00976 fss = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()/TRANS_PROB_CONST/POW2(Kss)/
00977 (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Lo()).g()*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Hi()).g();
00978
00979 ftt = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul()/TRANS_PROB_CONST/POW2(Ktt)/
00980 (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Lo()).g()*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Hi()).g();
00981
00982 temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
00983 fssNew = Kss*POW2( temp );
00984 temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
00985 fttNew = Ktt*POW2( temp );
00986 temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
00987 fstNew = Kst*POW2( temp );
00988 temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
00989 ftsNew = Kts*POW2( temp );
00990
00991 Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Lo()).g()/
00992 (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Hi()).g();
00993
00994 Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Lo()).g()/
00995 (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Hi()).g();
00996
00997 Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Lo()).g()/
00998 (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Hi()).g();
00999
01000 Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Lo()).g()/
01001 (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Hi()).g();
01002
01003 error = fabs( ( iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()+
01004 iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul() )/
01005 (Ass+Ast+Ats+Att) - 1.f );
01006
01007 if( error > 0.001 )
01008 {
01009 fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
01010 ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
01011 Ass/iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul(),
01012 Att/iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul(),
01013 Ast/iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Emis().Aul(),
01014 Ats/iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Emis().Aul() );
01015 }
01016
01017 iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul() = (realnum)Ass;
01018 iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul() = (realnum)Att;
01019 iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Emis().Aul() = (realnum)Ast;
01020 iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Emis().Aul() = (realnum)Ats;
01021
01022 return;
01023 }
01024
01025
01026
01027 STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
01028 {
01029
01030
01031
01032
01033
01034
01035
01036
01037 long lg;
01038 double fmean,mu,EinsteinA,w,RI2_cm;
01039
01040 DEBUG_ENTRY( "ritoa()" );
01041
01042 mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT));
01043
01044 w = 2. * PI * k * SPEEDLIGHT;
01045
01046 RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM;
01047
01048 lg = (lf > li) ? lf : li;
01049
01050 fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));
01051
01052 EinsteinA = TRANS_PROB_CONST*k*k*fmean;
01053
01054
01055
01056 return EinsteinA;
01057 }
01058
01059 realnum helike_transprob( long nelem, long ipHi, long ipLo )
01060 {
01061 double Aul, Aul1;
01062 double Enerwn, n_eff_hi, n_eff_lo;
01063 long ipISO = ipHE_LIKE;
01064
01065 double z4 = POW2((double)nelem);
01066 z4 *= z4;
01067
01068 DEBUG_ENTRY( "helike_transprob()" );
01069
01070 Enerwn = iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN();
01071 n_eff_hi = N_(ipHi) - helike_quantum_defect(nelem,ipHi);
01072 n_eff_lo = N_(ipLo) - helike_quantum_defect(nelem,ipLo);
01073
01074 if( ipHi >= iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max )
01075 {
01076 if( ipLo >= iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max )
01077 {
01078
01079 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
01080 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f,0.001f);
01081
01082 ASSERT( Aul > 0.);
01083 }
01084 else
01085 {
01086
01087
01088 Aul = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)+1,
01089 S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
01090
01091 iso_sp[ipISO][nelem].CachedAs[ N_(ipHi)-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][0] = (realnum)Aul;
01092
01093 Aul *= (2.*L_(ipLo)+3.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
01094
01095 if( L_(ipLo) != 0 )
01096 {
01097
01098 Aul1 = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)-1,
01099 S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
01100
01101 iso_sp[ipISO][nelem].CachedAs[ N_(ipHi)-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] = (realnum)Aul1;
01102
01103
01104 Aul += Aul1*(2.*L_(ipLo)-1.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
01105 }
01106 else
01107 iso_sp[ipISO][nelem].CachedAs[ N_(ipHi)-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] = 0.f;
01108
01109 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.01f);
01110 ASSERT( Aul > 0.);
01111 }
01112 }
01113 else
01114 {
01115
01116 if( Enerwn < 0. )
01117 {
01118 Aul = he_1trans( nelem, -1.*Enerwn, n_eff_lo,
01119 L_(ipLo), S_(ipLo), ipLo-3, n_eff_hi, L_(ipHi), S_(ipHi), ipHi-3);
01120 }
01121 else
01122 {
01123 Aul = he_1trans( nelem, Enerwn, n_eff_hi,
01124 L_(ipHi), S_(ipHi), ipHi-3, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
01125 }
01126 }
01127
01128 return (realnum)Aul;
01129 }
01130
01131
01132 void HelikeTransProbSetup( void )
01133 {
01134
01135 const int chLine_LENGTH = 1000;
01136 char chLine[chLine_LENGTH];
01137
01138 FILE *ioDATA;
01139 bool lgEOL;
01140
01141 long nelem, ipLo, ipHi, i, i1, i2, i3;
01142
01143 DEBUG_ENTRY( "HelikeTransProbSetup()" );
01144
01145 TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
01146
01147 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
01148 {
01149
01150 TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) );
01151
01152 for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo )
01153 {
01154 TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX );
01155 }
01156 }
01157
01158
01159
01160 if( trace.lgTrace )
01161 fprintf( ioQQQ," HelikeTransProbSetup opening he_transprob.dat:");
01162
01163 ioDATA = open_data( "he_transprob.dat", "r" );
01164
01165
01166 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01167 {
01168 fprintf( ioQQQ, " HelikeTransProbSetup could not read first line of he_transprob.dat.\n");
01169 cdEXIT(EXIT_FAILURE);
01170 }
01171 i = 1;
01172 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01173 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01174 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
01175 {
01176 fprintf( ioQQQ,
01177 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
01178 fprintf( ioQQQ,
01179 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
01180 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
01181 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01182 cdEXIT(EXIT_FAILURE);
01183 }
01184
01185
01186 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
01187 {
01188 for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
01189 {
01190 for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
01191 {
01192 TransProbs[nelem][ipHi][ipLo] = -1.;
01193 }
01194 }
01195 }
01196
01197 for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
01198 {
01199 char *chTemp;
01200
01201
01202 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01203 BadRead();
01204
01205 while( chLine[0]=='#' )
01206 {
01207 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01208 BadRead();
01209 }
01210
01211 i3 = 1;
01212 i1 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
01213 i2 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
01214
01215 if( i1<0 || i2<=i1 )
01216 {
01217 fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
01218 cdEXIT(EXIT_FAILURE);
01219 }
01220
01221 chTemp = chLine;
01222
01223
01224 for( i=0; i<1; ++i )
01225 {
01226 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
01227 {
01228 fprintf( ioQQQ, " HelikeTransProbSetup could not init he_transprob\n" );
01229 cdEXIT(EXIT_FAILURE);
01230 }
01231 ++chTemp;
01232 }
01233
01234
01235 for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
01236 {
01237 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
01238 {
01239 fprintf( ioQQQ, " HelikeTransProbSetup could not scan he_transprob\n" );
01240 cdEXIT(EXIT_FAILURE);
01241 }
01242 ++chTemp;
01243
01244 sscanf( chTemp , "%le" , &TransProbs[nelem][i2][i1] );
01245
01247 if( lgEOL )
01248 {
01249 fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
01250 cdEXIT(EXIT_FAILURE);
01251 }
01252 }
01253 }
01254
01255
01256 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01257 {
01258 fprintf( ioQQQ, " HelikeTransProbSetup could not read last line of he_transprob.dat.\n");
01259 cdEXIT(EXIT_FAILURE);
01260 }
01261 i = 1;
01262 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01263 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01264 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
01265 {
01266 fprintf( ioQQQ,
01267 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
01268 fprintf( ioQQQ,
01269 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
01270 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
01271 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01272 cdEXIT(EXIT_FAILURE);
01273 }
01274
01275
01276 fclose( ioDATA );
01277 return;
01278 }