cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*he_1trans compute Aul for given line */ 00004 /*ritoa - converts the square of the radial integral for a transition 00005 * (calculated by scqdri) to the transition probability, Aul. */ 00006 /*ForbiddenAuls calculates transition probabilities for forbidden transitions. */ 00007 /*scqdri - stands for Semi-Classical Quantum Defect Radial Integral */ 00008 /*Jint - used by scqdri */ 00009 /*AngerJ - used by scqdri */ 00010 /*DoFSMixing - applies a fine structure mixing approximation to A's. To be replaced by 00011 * method that treats the entire rate matrix. */ 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 /* the array of transitions probabilities read from data file. */ 00025 static double ***TransProbs; 00026 00027 /*ritoa converts the square of the radial integral for a transition 00028 * (calculated by scqdri) to the transition probability, Aul. */ 00029 STATIC double ritoa( long li, long lf, long nelem, double k, double RI2 ); 00030 00031 /* handles all forbidden transition probabilities. */ 00032 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem ); 00033 00034 /* 00035 static long RI_ipHi, RI_ipLo, RI_ipLev; 00036 static long globalZ; 00037 */ 00038 00039 /* used as parameters in qg32 integration */ 00040 static double vJint , zJint; 00041 00042 /* the integrand used in the AngerJ function below. */ 00043 STATIC double Jint( double theta ) 00044 { 00045 /* [ cos[vx - z sin[x]] ] */ 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 /* AngerJ function. */ 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 /* Estimate number of peaks in integrand. */ 00067 /* Divide region of integration by number */ 00068 /* peaks in region. */ 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 /* Semi-Classical Quantum Defect Radial Integral */ 00096 /* */ 00097 /* See for example */ 00098 /* Atomic, Molecular & Optical Physics Handbook */ 00099 /* Gordon W. F. Drake; Editor */ 00100 /* AIP Press */ 00101 /* Woddbury, New York. */ 00102 /* 1996 */ 00103 /* */ 00104 /* NOTE:: we do not include the Bohr Radius a_o in the */ 00105 /* definition of of R(n,L;n'L') as per Drake. */ 00106 /* */ 00107 /* */ 00108 /* 1 (n_c)^2 | { D_l max(L,L') } */ 00109 /* R(n,L;n'L') = --- ------- | { 1 - ------------- } J( D_n-1; -x ) - */ 00110 /* Z 2 D_n | { n_c } */ 00111 /* */ 00112 /* */ 00113 /* { D_L max(L,L') } */ 00114 /* - { 1 + ------------- } J( D_n+1; -x ) */ 00115 /* { n_c } */ 00116 /* */ 00117 /* */ 00118 /* 2 | */ 00119 /* + --- sin(Pi D_n)(1-e) | */ 00120 /* Pi | */ 00121 /* */ 00122 /* where */ 00123 /* n_c = (2n*'n*)/(n*'+n*) */ 00124 /* */ 00125 /* Here is the quantity Drake gives... */ 00126 /* n_c = ( 2.0 * nstar * npstar ) / ( nstar + npstar ); */ 00127 /* */ 00128 /* while V.A. Davidkin uses */ 00129 /* n_c = sqrt( nstar * npstar ); */ 00130 /* */ 00131 /* D_n = n*' - n* */ 00132 /* */ 00133 /* D_L = L*' - L* */ 00134 /* */ 00135 /* x = e D_n */ 00136 /* */ 00137 /* Lmx = max(L',L) */ 00138 /* */ 00139 /* e = sqrt( 1 - {Lmx/n_c}^2 ) */ 00140 /* */ 00141 /* */ 00142 /* Here n* = n - qd where qd is the quantum defect */ 00143 /* */ 00144 /******************************************************************************/ 00145 /******************************************************************************/ 00146 double scqdri(/* upper and lower quantum numbers...n's are effective */ 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 /* This can happen for certain quantum defects */ 00188 /* in the lower n=1:l=0 state. In which case you */ 00189 /* probably should be using some other alogrithm */ 00190 /* or theory to calculate the dipole moment. */ 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 /* >>refer Helike 2pho Derevianko, A., & Johnson, W.R. 1997, Phys. Rev. A 56, 1288 00225 * numbers are not explicitly given in this paper for Z=21-24,26,27,and 29. 00226 * So numbers given here are interpolated. */ 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 /* Important clarification, according to Derevianko & Johnson (see ref above), 2^3S can decay 00231 * to ground in one of two ways: through a two-photon process, or through a single-photon M1 decay, 00232 * but the M1 rates are about 10^4 greater that the two-photon decays throughout the entire 00233 * sequence. Thus these numbers, are much weaker than the effective decay rate, but should probably 00234 * be treated in as a two-photon decay at some point */ 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 /* All of these but the second and third one (values 51.02 and 1E-20) are from 00248 * >>refer HeI As Lach, G, & Pachucki, K, 2001, Phys. Rev. A 64, 042510 00249 ,* 1E-20 is made up 00250 * 51.3 is from the Derevianko & Johnson paper cited above. */ 00251 double ForbiddenHe[5] = { 1.272E-4, 51.02, 1E-20, 177.58, 0.327 }; 00252 00253 fixit(); // adding the 2-photon decay 2^3S - 1^1S may be important in early universe 00254 A = ForbiddenHe[ipHi - 1]; 00255 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00256 } 00257 else 00258 { 00259 switch ( (int)ipHi ) 00260 { 00261 case 1: /* Parameters for 2^3S to ground transition. */ 00262 /* >>refer Helike As Lin, C.D., Johnson, W.R., and Dalgarno, A. 1977, 00263 * >>refercon Phys. Rev. A 15, 1, 010015 */ 00264 A = (3.9061E-7) * pow( (double)nelem+1., 10.419 ) + As2nuFrom3S[nelem-2]; 00265 break; 00266 case 2: /* Parameters for 2^1S to ground transition. */ 00267 A = As2nuFrom1S[nelem-2]; 00268 break; 00269 case 3: /* Parameters for 2^3P0 to ground transition. */ 00270 A = iso.SmallA; 00271 break; 00272 case 4: /* Parameters for 2^3P1 to ground transition. */ 00273 /* >>chng 06 jul 25, only use the fit to Johnson et al. values up to and 00274 * including argon, where there calculation stops. For higher Z use below */ 00275 if( nelem <= ipARGON ) 00276 { 00277 A = ( 11.431 * pow((double)nelem, 9.091) ); 00278 } 00279 else 00280 { 00281 /* a fit to the Lin et al. 1977. values, which go past zinc. */ 00282 A = ( 383.42 * pow((double)nelem, 7.8901) ); 00283 } 00284 break; 00285 case 5: /* Parameters for 2^3P2 to ground transition. */ 00286 /* fit to Lin et al. 1977 values. This fit is good 00287 * to 7% for the range from carbon to iron. The Lin et al. values 00288 * differ from the Hata and Grant (1984) values (only given from 00289 * calcium to copper) by less than 2%. */ 00290 A = ( 0.11012 * pow((double)nelem, 7.6954) ); 00291 break; 00292 default: 00293 TotalInsanity(); 00294 } 00295 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00296 } 00297 /* For now just just put 1% error for forbidden lines. */ 00298 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,.01f); 00299 return A; 00300 } 00301 00302 /* The next two cases are fits to probabilities given in 00303 * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., & 00304 * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */ 00305 /* originally astro.ph. 0201454 */ 00306 /* The involve Triplet P and Singlet S. Rates for Triplet S to Singlet P 00307 * do not seem to be available. */ 00308 00309 /* Triplet P to Singlet S...Delta n not equal to zero! */ 00310 else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==3 && 00311 L_(ipLo)==0 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) ) 00312 { 00313 A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) / 00314 pow((double)N_(ipHi),2.877); 00315 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00316 } 00317 00318 /* Singlet S to Triplet P...Delta n not equal to zero! */ 00319 else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==1 && 00320 L_(ipLo)==1 && S_(ipLo)==3 && N_(ipLo) < N_(ipHi) ) 00321 { 00322 A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111); 00323 00324 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) ) 00325 { 00326 /* This is divided by 3 instead of 9, because value calculated is specifically for 2^3P1. 00327 * Here we assume statistical population of the other two. */ 00328 A *= (2.*(ipLo-3)+1.0)/3.0; 00329 } 00330 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00331 } 00332 00333 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe3s3S ) ) 00334 { 00335 double As_3TripS_to_2TripS[9] = { 00336 7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2, 00337 1.27E-1, 5.56E-1, 2.01E0, 6.26E0 }; 00338 00339 /* These M1 transitions given by 00340 * >>refer He-like As Savukov, I.M., Labzowsky, and Johnson, W.R. 2005, PhysRevA, 72, 012504 */ 00341 if( nelem <= ipNEON ) 00342 { 00343 A = As_3TripS_to_2TripS[nelem-1]; 00344 /* 20% error is based on difference between Savukov, Labzowsky, and Johnson (2005) 00345 * and Lach and Pachucki (2001) for the helium transition. */ 00346 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f); 00347 } 00348 else 00349 { 00350 /* This is an extrapolation to the data given above. The fit reproduces 00351 * the above values to 10% or better. */ 00352 A= 7.22E-9*pow((double)nelem, 9.33); 00353 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f); 00354 } 00355 } 00356 00357 else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) ) 00358 { 00359 /* This transition,1.549 , given by Lach and Pachucki, 2001 for the atom */ 00360 if( nelem == ipHELIUM ) 00361 { 00362 A = 1.549; 00363 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00364 } 00365 else 00366 { 00367 /* This is a fit to data given in 00368 * >>refer He-like As Savukov, I.M., Johnson, W.R., & Safronova, U.I. 00369 * >>refercon astro-ph 0205163 */ 00370 A= 0.1834*pow((double)nelem, 6.5735); 00371 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00372 } 00373 } 00374 00375 else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 ) 00376 { 00377 /* >>refer He As Bulatov, N.N. 1976, Soviet Astronomy, 20, 315 */ 00378 fixit(); 00379 /* This is the 29.6 GHz line that can be seen in radio galaxies. */ 00384 A = 5.78E-12; 00385 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00386 00387 } 00388 00389 else if( nelem==ipHELIUM && ipHi==5 && ipLo==4 ) 00390 { 00391 fixit(); 00392 /* This is the 3 GHz line that can be seen in radio galaxies. */ 00397 A = 3.61E-15; 00398 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00399 } 00400 00401 else 00402 { 00403 /* Current transition is not supported. */ 00404 A = iso.SmallA; 00405 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00406 } 00407 00408 /* For now just put 1% error for forbidden lines. */ 00409 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,.01f); 00410 00411 ASSERT( A > 0.); 00412 00413 return A; 00414 } 00415 00416 /* Calculates Einstein A for a given transition. */ 00417 double he_1trans( 00418 /* charge on c scale, Energy is wavenumbers, Einstein A */ 00419 long nelem , double Enerwn , 00420 /* quantum numbers of upper level: */ 00421 double Eff_nupper, long lHi, long sHi, long jHi, 00422 /* and of lower level: */ 00423 double Eff_nlower, long lLo, long sLo, long jLo ) 00424 /* Note j is only necessary for 2 triplet P...for all other n,l,s, 00425 * j is completely ignored. */ 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 /* Since 0.4 is bigger than any defect, adding that to the effective principle quantum number, 00436 * and truncating to an integer will produce the principal quantum number. */ 00437 nHi = (int)(Eff_nupper + 0.4); 00438 nLo = (int)(Eff_nlower + 0.4); 00439 00440 /* Make sure this worked correctly. */ 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( StatesElem[ipHE_LIKE][nelem][ipHi].n == nHi ); 00459 if( nHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] ) 00460 { 00461 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].l == lHi ); 00462 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].S == sHi ); 00463 } 00464 ASSERT( StatesElem[ipHE_LIKE][nelem][ipLo].n == nLo ); 00465 if( nLo <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] ) 00466 { 00467 ASSERT( StatesElem[ipHE_LIKE][nelem][ipLo].l == lLo ); 00468 ASSERT( StatesElem[ipHE_LIKE][nelem][ipLo].S == sLo ); 00469 } 00470 00471 /* First do allowed transitions */ 00472 if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) ) 00473 { 00474 Aul = -2.; 00475 00476 /* For clarity, let's do this in two separate chunks...one for helium, one for everything else. */ 00477 if( nelem == ipHELIUM ) 00478 { 00479 /* Retrieve transition probabilities for Helium. */ 00480 /* >>refer He As Drake, G.W.F., Atomic, Molecular, and Optical Physics Handbook */ 00481 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] ) 00482 { 00483 /*Must not be accessed by collapsed levels! */ 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.0005f); 00491 } 00492 00493 if( Aul < 0. ) 00494 { 00495 /* Here are the Lyman transitions. */ 00496 if( ipLo == ipHe1s1S ) 00497 { 00498 ASSERT( (lHi == 1) && (sHi == 1) ); 00499 00500 /* these fits calculated from Drake A's (1996) */ 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.005f); 00505 } 00506 00507 /* last resort for transitions involving significant defects, 00508 * except that highest lLo are excluded */ 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 00514 if( lHi + lLo >= 7 ) 00515 { 00516 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f); 00517 } 00518 else 00519 { 00520 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f); 00521 } 00522 } 00523 /* These fits come from extrapolations of Drake's oscillator strengths 00524 * out to the series limit. We also use this method to obtain threshold 00525 * photoionization cross-sections for the lower level of each transition here. */ 00526 /* See these two references : 00527 * >>refer He As Hummer, D. G. \& Storey, P. J. 1998, MNRAS 297, 1073 00528 * >>refer Seaton's Theorem Seaton, M. J. 1958, MNRAS 118, 504 */ 00529 else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 ) 00530 { 00531 int paramSet=0; 00532 double emisOscStr, x, a, b, c; 00533 double extrapol_Params[2][4][4][3] = { 00534 /* these are for singlets */ 00535 { 00536 { /* these are P to S */ 00537 { 0.8267396 , 1.4837624 , -0.4615955 }, 00538 { 1.2738405 , 1.5841806 , -0.3022984 }, 00539 { 1.6128996 , 1.6842538 , -0.2393057 }, 00540 { 1.8855491 , 1.7709125 , -0.2115213 }}, 00541 { /* these are S to P */ 00542 { -1.4293664 , 2.3294080 , -0.0890470 }, 00543 { -0.3608082 , 2.3337636 , -0.0712380 }, 00544 { 0.3027974 , 2.3326252 , -0.0579008 }, 00545 { 0.7841193 , 2.3320138 , -0.0497094 }}, 00546 { /* these are D to P */ 00547 { 1.1341403 , 3.1702435 , -0.2085843 }, 00548 { 1.7915926 , 2.4942946 , -0.2266493 }, 00549 { 2.1979400 , 2.2785377 , -0.1518743 }, 00550 { 2.5018229 , 2.1925720 , -0.1081966 }}, 00551 { /* these are P to D */ 00552 { 0.0000000 , 0.0000000 , 0.0000000 }, 00553 { -2.6737396 , 2.9379143 , -0.3805367 }, 00554 { -1.4380124 , 2.7756396 , -0.2754625 }, 00555 { -0.6630196 , 2.6887253 , -0.2216493 }}, 00556 }, 00557 /* these are for triplets */ 00558 { 00559 { /* these are P to S */ 00560 { 0.3075287 , 0.9087130 , -1.0387207 }, 00561 { 0.687069 , 1.1485864 , -0.6627317 }, 00562 { 0.9776064 , 1.3382024 , -0.5331906 }, 00563 { 1.2107725 , 1.4943721 , -0.4779232 }}, 00564 { /* these are S to P */ 00565 { -1.3659605 , 2.3262253 , -0.0306439 }, 00566 { -0.2899490 , 2.3279391 , -0.0298695 }, 00567 { 0.3678878 , 2.3266603 , -0.0240021 }, 00568 { 0.8427457 , 2.3249540 , -0.0194091 }}, 00569 { /* these are D to P */ 00570 { 1.3108281 , 2.8446367 , -0.1649923 }, 00571 { 1.8437692 , 2.2399326 , -0.2583398 }, 00572 { 2.1820792 , 2.0693762 , -0.1864091 }, 00573 { 2.4414052 , 2.0168255 , -0.1426083 }}, 00574 { /* these are P to D */ 00575 { 0.0000000 , 0.0000000 , 0.0000000 }, 00576 { -1.9219877 , 2.7689624 , -0.2536072 }, 00577 { -0.7818065 , 2.6595150 , -0.1895313 }, 00578 { -0.0665624 , 2.5955623 , -0.1522616 }}, 00579 } 00580 }; 00581 00582 if( lLo==0 ) 00583 { 00584 paramSet = 0; 00585 } 00586 else if( lLo==1 && lHi==0 ) 00587 { 00588 paramSet = 1; 00589 } 00590 else if( lLo==1 && lHi==2 ) 00591 { 00592 paramSet = 2; 00593 } 00594 else if( lLo==2 ) 00595 { 00596 paramSet = 3; 00597 ASSERT( lHi==1 ); 00598 } 00599 00600 ASSERT( (int)((sHi-1)/2) == 0 || (int)((sHi-1)/2) == 1 ); 00601 a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0]; 00602 b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1]; 00603 c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2]; 00604 ASSERT( Enerwn > 0. ); 00605 x = log( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]*RYD_INF/Enerwn ); 00606 00607 emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)* 00608 (2.*lLo+1)/(2.*lHi+1); 00609 00610 Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr; 00611 00612 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) 00613 { 00614 Aul *= (2.*(ipLo-3)+1.0)/9.0; 00615 } 00616 00617 ASSERT( Aul > 0. ); 00618 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f); 00619 } 00620 else 00621 { 00622 /* Calculate the radial integral from the quantum defects. */ 00623 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM)); 00624 ASSERT( RI2 > 0. ); 00625 /* Convert radial integral to Aul. */ 00626 Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2); 00627 /* radial integral routine does not recognize fine structure. 00628 * Here we split 2^3P. */ 00629 if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) 00630 { 00631 Aul *= (2.*(ipLo-3)+1.0)/9.0; 00632 } 00633 00634 ASSERT( Aul > 0. ); 00635 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.03f); 00636 } 00637 } 00638 } 00639 00640 /* Heavier species */ 00641 else 00642 { 00643 /* Retrieve transition probabilities for Helike ions. */ 00644 /* >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., & 00645 * >>refercon Dalgarno, A., 2002, ApJS 141, 543J, originally astro.ph. 0201454 */ 00646 if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] ) 00647 { 00648 /*Must not be accessed by collapsed levels! */ 00649 ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) ); 00650 ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) ); 00651 ASSERT( ipHi > 2 ); 00652 00653 Aul = TransProbs[nelem][ipHi][ipLo]; 00654 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00655 } 00656 00657 if( Aul < 0. ) 00658 { 00659 /* Do same-n transitions. */ 00660 if( nLo==nHi && lHi<=2 && lLo<=2 ) 00661 { 00662 /* These are 2p3Pj to 2s3S fits to (low Z) Porquet & Dubau (2000) & 00663 * (high Z) NIST Atomic Spectra Database. */ 00664 if( ipLo == ipHe2s3S ) 00665 { 00666 if(ipHi == ipHe2p3P0) 00667 Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76); 00668 else if(ipHi == ipHe2p3P1) 00669 Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76); 00670 else if(ipHi == ipHe2p3P2) 00671 Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29); 00672 else 00673 { 00674 fprintf(ioQQQ," PROBLEM Impossible value for ipHi in he_1trans.\n"); 00675 TotalInsanity(); 00676 } 00677 } 00678 00679 /* These are 2p1P to 2s1S fits to data from TOPbase. */ 00680 else if( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) ) 00681 { 00682 Aul = 5.53E6 * exp( 0.171*(nelem+1.) ); 00683 } 00684 00685 else 00686 { 00687 /* This case should only be entered if n > 2. Those cases were done above. */ 00688 ASSERT( nLo > 2 ); 00689 00690 /* Triplet P to triplet S, delta n = 0 */ 00691 if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3)) 00692 { 00693 Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328); 00694 } 00695 /* Singlet P to singlet D, delta n = 0 */ 00696 else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1)) 00697 { 00698 Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269); 00699 } 00700 /* Singlet P to singlet S, delta n = 0 */ 00701 else if( (lHi == 1) && (sHi == 1) && (lLo == 0) ) 00702 { 00703 Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077); 00704 } 00705 else 00706 { 00707 ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) ); 00708 Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9); 00709 if( (lHi >2) || (lLo > 2) ) 00710 Aul *= (lHi/2.); 00711 if(lLo > 2) 00712 Aul *= (1./9.); 00713 } 00714 } 00715 ASSERT( Aul > 0.); 00716 } 00717 00718 /* assume transitions involving F and higher orbitals are hydrogenic. */ 00719 else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) ) 00720 { 00721 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem); 00722 ASSERT( Aul > 0.); 00723 } 00724 00725 /* These transitions are of great importance, but the below radial integral 00726 * routine fails to achieve desirable accuracy, so these are fits as produced 00727 * from He A's for nupper through 9. They are transitions to ground and 00728 * 2, 3, and 4 triplet S. */ 00729 else if( ( ipLo == 0 ) || ( ipLo == ipHe2s1S ) || ( ipLo == ipHe2s3S ) 00730 || ( ipLo == ipHe3s3S ) || ( ipLo == ipHe4s3S ) ) 00731 { 00732 /* Here are the Lyman transitions. */ 00733 if( ipLo == 0 ) 00734 { 00735 ASSERT( (lHi == 1) && (sHi == 1) ); 00736 00737 /* In theory, this Z dependence should be Z^4, but values from TOPbase 00738 * suggest 3.9 is a more accurate exponent. Values from 00739 * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., & 00740 * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */ 00741 /* originally astro.ph. 0201454 */ 00742 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1); 00743 } 00744 00745 /* Here are the Balmer transitions. */ 00746 else if( ipLo == ipHe2s1S ) 00747 { 00748 ASSERT( (lHi == 1) && (sHi == 1) ); 00749 00750 /* More fits, as above. */ 00751 Aul = 5.0e8 * pow((double)nelem,4.) / pow((double)nHi, 2.889); 00752 } 00753 00754 /* Here are transitions down to triplet S */ 00755 else 00756 { 00757 ASSERT( (lHi == 1) && (sHi == 3) ); 00758 00759 /* These fits also as above. */ 00760 if( nLo == 2 ) 00761 Aul = 1.5 * 3.405E8 * pow((double)nelem,4.) / pow((double)nHi, 2.883); 00762 else if( nLo == 3 ) 00763 Aul = 2.5 * 4.613E7 * pow((double)nelem,4.) / pow((double)nHi, 2.672); 00764 else 00765 Aul = 3.0 * 1.436E7 * pow((double)nelem,4.) / pow((double)nHi, 2.617); 00766 } 00767 00768 ASSERT( Aul > 0.); 00769 } 00770 00771 /* Every other allowed transition is calculated as follows. */ 00772 else 00773 { 00774 /* Calculate the radial integral from the quantum defects. */ 00775 RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem)); 00776 /* Convert radial integral to Aul. */ 00777 Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2); 00778 /* radial integral routine does not recognize fine structure. 00779 * Here we split 2^3P. */ 00780 if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) && (Aul > iso.SmallA) ) 00781 { 00782 Aul *= (2.*(ipLo-3)+1.0)/9.0; 00783 } 00784 00785 } 00786 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 00787 } 00788 00789 /* for now just give ions some a 5% error across the board */ 00790 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.05f); 00791 } 00792 } 00793 00794 /* Now do forbidden transitions from 2-1 ... */ 00795 /* and those given by 00796 * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., & 00797 * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */ 00798 /* originally astro.ph. 0201454 00799 * for heavy elements. These are triplet P to singlet S, 00800 * going either up or down...Triplet S to Singlet P are not included, as they are far weaker. */ 00801 else 00802 { 00803 ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) ); 00804 Aul = ForbiddenAuls(ipHi, ipLo, nelem); 00805 ASSERT( Aul > 0. ); 00806 } 00807 00808 Aul = MAX2( Aul, iso.SmallA ); 00809 ASSERT( Aul >= iso.SmallA ); 00810 00811 /* negative energy for a transition with substantial transition probability 00812 * would be major logical error - but is ok for same-n l transitions */ 00813 if( Enerwn < 0. && Aul > iso.SmallA ) 00814 { 00815 fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn ); 00816 } 00817 00818 return Aul; 00819 } 00820 00821 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing ) 00822 { 00823 /* Every bit of this routine is based upon the singlet-triplet mixing formalism given in 00824 * >>refer He FSM Drake, G. W. F. 1996, in Atomic, Molecular, \& Optical Physics Handbook, 00825 * >>refercon ed. G. W. F. Drake (New York: AIP Press). 00826 * That formalism mixes the levels themselves, but since this code is not J-resolved, we simulate 00827 * that by mixing only the transition probabilities. We find results comparable to those calculated 00828 * in the fully J-resolved model spearheaded by Rob Bauman, and described in 00829 * >>refer He FSM Bauman, R. P., Porter, R. L., Ferland, G. J., \& MacAdam, K. B. 2005, ApJ, accepted */ 00830 long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip; 00831 double Ass, Att, Ast, Ats; 00832 double SinHi, SinLo, CosHi, CosLo; 00833 double HiMixingAngle, LoMixingAngle , error; 00834 double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp; 00835 00836 DEBUG_ENTRY( "DoFSMixing()" ); 00837 00838 nHi = StatesElem[ipHE_LIKE][nelem][ipHiSing].n; 00839 lHi = StatesElem[ipHE_LIKE][nelem][ipHiSing].l; 00840 sHi = StatesElem[ipHE_LIKE][nelem][ipHiSing].S; 00841 nLo = StatesElem[ipHE_LIKE][nelem][ipLoSing].n; 00842 lLo = StatesElem[ipHE_LIKE][nelem][ipLoSing].l; 00843 sLo = StatesElem[ipHE_LIKE][nelem][ipLoSing].S; 00844 00845 if( ( sHi == 3 || sLo == 3 ) || 00846 ( abs(lHi - lLo) != 1 ) || 00847 ( nLo < 2 ) || 00848 ( lHi <= 1 || lLo <= 1 ) || 00849 ( nHi == nLo && lHi == 1 && lLo == 2 ) || 00850 ( nHi > nLo && lHi != 1 && lLo != 1 ) ) 00851 { 00852 return; 00853 } 00854 00855 ASSERT( lHi > 0 ); 00856 /*ASSERT( (lHi > 1) && (lLo > 1) );*/ 00857 00858 ipHiTrip = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nHi][lHi][3]; 00859 ipLoTrip = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nLo][lLo][3]; 00860 00861 if( lHi == 2 ) 00862 { 00863 HiMixingAngle = 0.01; 00864 } 00865 else if( lHi == 3 ) 00866 { 00867 HiMixingAngle = 0.5; 00868 } 00869 else 00870 { 00871 HiMixingAngle = PI/4.; 00872 } 00873 00874 if( lLo == 2 ) 00875 { 00876 LoMixingAngle = 0.01; 00877 } 00878 else if( lLo == 3 ) 00879 { 00880 LoMixingAngle = 0.5; 00881 } 00882 else 00883 { 00884 LoMixingAngle = PI/4.; 00885 } 00886 00887 /* These would not work correctly if l<=1 were included in this treatment! */ 00888 ASSERT( ipHiTrip > ipLoTrip ); 00889 ASSERT( ipHiTrip > ipLoSing ); 00890 ASSERT( ipHiSing > ipLoTrip ); 00891 ASSERT( ipHiSing > ipLoSing ); 00892 00893 SinHi = sin( HiMixingAngle ); 00894 SinLo = sin( LoMixingAngle ); 00895 CosHi = cos( HiMixingAngle ); 00896 CosLo = cos( LoMixingAngle ); 00897 00898 Kss = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].EnergyWN; 00899 Ktt = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].EnergyWN; 00900 Kst = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].EnergyWN; 00901 Kts = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].EnergyWN; 00902 00903 fss = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul/TRANS_PROB_CONST/POW2(Kss)/ 00904 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Lo->g*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Hi->g; 00905 00906 ftt = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul/TRANS_PROB_CONST/POW2(Ktt)/ 00907 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Lo->g*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Hi->g; 00908 00909 temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo; 00910 fssNew = Kss*POW2( temp ); 00911 temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo; 00912 fttNew = Ktt*POW2( temp ); 00913 temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo; 00914 fstNew = Kst*POW2( temp ); 00915 temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo; 00916 ftsNew = Kts*POW2( temp ); 00917 00918 Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Lo->g/ 00919 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Hi->g; 00920 00921 Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Lo->g/ 00922 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Hi->g; 00923 00924 Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Lo->g/ 00925 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Hi->g; 00926 00927 Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Lo->g/ 00928 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Hi->g; 00929 00930 error = fabs( ( Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul+ 00931 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul )/ 00932 (Ass+Ast+Ats+Att) - 1.f ); 00933 00934 if( error > 0.001 ) 00935 { 00936 fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error, 00937 ipLoSing, ipHiSing, ipLoTrip, ipHiTrip, 00938 Ass/Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul, 00939 Att/Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul, 00940 Ast/Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Emis->Aul, 00941 Ats/Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Emis->Aul ); 00942 } 00943 00944 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul = (realnum)Ass; 00945 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul = (realnum)Att; 00946 Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Emis->Aul = (realnum)Ast; 00947 Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Emis->Aul = (realnum)Ats; 00948 00949 return; 00950 } 00951 00952 /*ritoa converts the square of the radial integral for a transition 00953 * (calculated by scqdri) to the transition probability, Aul. */ 00954 STATIC double ritoa(long li, long lf, long nelem, double k, double RI2) 00955 { 00956 /* Variables are as follows: */ 00957 /* lg = larger of li and lf */ 00958 /* fmean = mean oscillator strength */ 00959 /* for a given level. */ 00960 /* mu = reduced mass of optical electron. */ 00961 /* EinsteinA = Einstein emission coef. */ 00962 /* w = angular frequency of transition. */ 00963 /* RI2_cm = square of rad. int. in cm^2. */ 00964 long lg; 00965 double fmean,mu,EinsteinA,w,RI2_cm; 00966 00967 DEBUG_ENTRY( "ritoa()" ); 00968 00969 mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT)); 00970 00971 w = 2. * PI * k * SPEEDLIGHT; 00972 00973 RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM; 00974 00975 lg = (lf > li) ? lf : li; 00976 00977 fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0)); 00978 00979 EinsteinA = TRANS_PROB_CONST*k*k*fmean; 00980 00981 /* ASSERT( EinsteinA > SMALLFLOAT ); */ 00982 00983 return EinsteinA; 00984 } 00985 00986 realnum helike_transprob( long nelem, long ipHi, long ipLo ) 00987 { 00988 double Aul, Aul1; 00989 double Enerwn, n_eff_hi, n_eff_lo; 00990 long ipISO = ipHE_LIKE; 00991 /* charge to 4th power, needed for scaling laws for As*/ 00992 double z4 = POW2((double)nelem); 00993 z4 *= z4; 00994 00995 DEBUG_ENTRY( "helike_transprob()" ); 00996 00997 Enerwn = Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN; 00998 n_eff_hi = N_(ipHi) - helike_quantum_defect(nelem,ipHi); 00999 n_eff_lo = N_(ipLo) - helike_quantum_defect(nelem,ipLo); 01000 01001 if( ipHi >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] ) 01002 { 01003 if( ipLo >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] ) 01004 { 01005 /* Neither upper nor lower is resolved...Aul's are easy. */ 01006 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4; 01007 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f); 01008 01009 ASSERT( Aul > 0.); 01010 } 01011 else 01012 { 01013 /* Lower level resolved, upper not. First calculate Aul 01014 * from upper level with ang mom one higher. */ 01015 Aul = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)+1, 01016 S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3); 01017 01018 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][0] = (realnum)Aul; 01019 01020 Aul *= (2.*L_(ipLo)+3.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi)); 01021 01022 if( L_(ipLo) != 0 ) 01023 { 01024 /* For all l>0, add in transitions from upper level with ang mom one lower. */ 01025 Aul1 = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)-1, 01026 S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3); 01027 01028 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = (realnum)Aul1; 01029 01030 /* now add in this part after multiplying by stat weight for lHi = lLo-1. */ 01031 Aul += Aul1*(2.*L_(ipLo)-1.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi)); 01032 } 01033 else 01034 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = 0.f; 01035 01036 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f); 01037 ASSERT( Aul > 0.); 01038 } 01039 } 01040 else 01041 { 01042 /* Both levels are resolved...do the normal bit. */ 01043 if( Enerwn < 0. ) 01044 { 01045 Aul = he_1trans( nelem, -1.*Enerwn, n_eff_lo, 01046 L_(ipLo), S_(ipLo), ipLo-3, n_eff_hi, L_(ipHi), S_(ipHi), ipHi-3); 01047 } 01048 else 01049 { 01050 Aul = he_1trans( nelem, Enerwn, n_eff_hi, 01051 L_(ipHi), S_(ipHi), ipHi-3, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3); 01052 } 01053 } 01054 01055 return (realnum)Aul; 01056 } 01057 01058 /* This routine is pure infrastructure and bookkeeping, no physics. */ 01059 void HelikeTransProbSetup( void ) 01060 { 01061 01062 const int chLine_LENGTH = 1000; 01063 char chLine[chLine_LENGTH]; 01064 01065 FILE *ioDATA; 01066 bool lgEOL; 01067 01068 long nelem, ipLo, ipHi, i, i1, i2, i3; 01069 01070 DEBUG_ENTRY( "HelikeTransProbSetup()" ); 01071 01072 TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM ); 01073 01074 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem ) 01075 { 01076 01077 TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) ); 01078 01079 for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo ) 01080 { 01081 TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX ); 01082 } 01083 } 01084 01085 /********************************************************************/ 01086 /*************** Read in data from he_transprob.dat *****************/ 01087 if( trace.lgTrace ) 01088 fprintf( ioQQQ," HelikeTransProbSetup opening he_transprob.dat:"); 01089 01090 ioDATA = open_data( "he_transprob.dat", "r" ); 01091 01092 /* check that magic number is ok */ 01093 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01094 { 01095 fprintf( ioQQQ, " HelikeTransProbSetup could not read first line of he_transprob.dat.\n"); 01096 cdEXIT(EXIT_FAILURE); 01097 } 01098 i = 1; 01099 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 01100 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 01101 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB ) 01102 { 01103 fprintf( ioQQQ, 01104 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" ); 01105 fprintf( ioQQQ, 01106 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" , 01107 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 ); 01108 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 01109 cdEXIT(EXIT_FAILURE); 01110 } 01111 01112 /* Initialize TransProbs[nelem][ipLo][ipHi] before filling it in. */ 01113 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 01114 { 01115 for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ ) 01116 { 01117 for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ ) 01118 { 01119 TransProbs[nelem][ipHi][ipLo] = -1.; 01120 } 01121 } 01122 } 01123 01124 for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ ) 01125 { 01126 char *chTemp; 01127 01128 /* get next line image */ 01129 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01130 BadRead(); 01131 01132 while( chLine[0]=='#' ) 01133 { 01134 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01135 BadRead(); 01136 } 01137 01138 i3 = 1; 01139 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL); 01140 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL); 01141 /* check that these numbers are correct */ 01142 if( i1<0 || i2<=i1 ) 01143 { 01144 fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n"); 01145 cdEXIT(EXIT_FAILURE); 01146 } 01147 01148 chTemp = chLine; 01149 01150 /* skip over 2 tabs to start of data */ 01151 for( i=0; i<1; ++i ) 01152 { 01153 if( (chTemp = strchr( chTemp, '\t' )) == NULL ) 01154 { 01155 fprintf( ioQQQ, " HelikeTransProbSetup could not init he_transprob\n" ); 01156 cdEXIT(EXIT_FAILURE); 01157 } 01158 ++chTemp; 01159 } 01160 01161 /* now read in the data */ 01162 for( nelem = ipHELIUM; nelem < LIMELM; nelem++ ) 01163 { 01164 if( (chTemp = strchr( chTemp, '\t' )) == NULL ) 01165 { 01166 fprintf( ioQQQ, " HelikeTransProbSetup could not scan he_transprob\n" ); 01167 cdEXIT(EXIT_FAILURE); 01168 } 01169 ++chTemp; 01170 01171 sscanf( chTemp , "%le" , &TransProbs[nelem][i2][i1] ); 01172 01174 if( lgEOL ) 01175 { 01176 fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n"); 01177 cdEXIT(EXIT_FAILURE); 01178 } 01179 } 01180 } 01181 01182 /* check that ending magic number is ok */ 01183 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01184 { 01185 fprintf( ioQQQ, " HelikeTransProbSetup could not read last line of he_transprob.dat.\n"); 01186 cdEXIT(EXIT_FAILURE); 01187 } 01188 i = 1; 01189 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 01190 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 01191 if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB ) 01192 { 01193 fprintf( ioQQQ, 01194 " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" ); 01195 fprintf( ioQQQ, 01196 " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" , 01197 TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 ); 01198 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 01199 cdEXIT(EXIT_FAILURE); 01200 } 01201 01202 /* close the data file */ 01203 fclose( ioDATA ); 01204 return; 01205 }