cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
helike_einsta.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*he_1trans compute Aul for given line */
4 /*ritoa - converts the square of the radial integral for a transition
5  * (calculated by scqdri) to the transition probability, Aul. */
6 /*ForbiddenAuls calculates transition probabilities for forbidden transitions. */
7 /*scqdri - stands for Semi-Classical Quantum Defect Radial Integral */
8 /*Jint - used by scqdri */
9 /*AngerJ - used by scqdri */
10 /*DoFSMixing - applies a fine structure mixing approximation to A's. To be replaced by
11  * method that treats the entire rate matrix. */
12 
13 #include "cddefines.h"
14 #include "physconst.h"
15 #include "taulines.h"
16 #include "dense.h"
17 #include "trace.h"
18 #include "hydro_bauman.h"
19 #include "iso.h"
20 #include "helike.h"
21 #include "helike_einsta.h"
22 #include "hydroeinsta.h"
23 
24 /* the array of transitions probabilities read from data file. */
25 static double ***TransProbs;
26 
27 /*ritoa converts the square of the radial integral for a transition
28  * (calculated by scqdri) to the transition probability, Aul. */
29 STATIC double ritoa( long li, long lf, long nelem, double k, double RI2 );
30 
31 /* handles all forbidden transition probabilities. */
32 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem );
33 
34 /*
35 static long RI_ipHi, RI_ipLo, RI_ipLev;
36 static long globalZ;
37 */
38 
39 /* used as parameters in qg32 integration */
40 static double vJint , zJint;
41 
42 /* the integrand used in the AngerJ function below. */
43 STATIC double Jint( double theta )
44 {
45  /* [ cos[vx - z sin[x]] ] */
46  double
47  d0 = ( 1.0 / PI ),
48  d1 = vJint * theta,
49  d2 = zJint * sin(theta),
50  d3 = (d1 - d2),
51  d4 = cos(d3),
52  d5 = (d0 * d4);
53 
54  return( d5 );
55 }
56 
57 /* AngerJ function. */
58 STATIC double AngerJ( double vv, double zz )
59 {
60  long int rep = 0, ddiv, divsor;
61 
62  double y = 0.0;
63 
64  DEBUG_ENTRY( "AngerJ()" );
65 
66  /* Estimate number of peaks in integrand. */
67  /* Divide region of integration by number */
68  /* peaks in region. */
69  if( (fabs(vv)) - (int)(fabs(vv)) > 0.5 )
70  ddiv = (int)(fabs(vv)) + 1;
71  else
72  ddiv = (int)(fabs(vv));
73 
74  divsor = ((ddiv == 0) ? 1 : ddiv);
75  vJint = vv;
76  zJint = zz;
77 
78  for( rep = 0; rep < divsor; rep++ )
79  {
80  double
81  rl = (((double) rep)/((double) divsor)),
82  ru = (((double) (rep+1))/((double) divsor)),
83  x_low = (PI * rl),
84  x_up = (PI * ru);
85 
86  y += qg32( x_low, x_up, Jint );
87  }
88 
89  return( y );
90 }
91 
92 /******************************************************************************/
93 /******************************************************************************/
94 /* */
95 /* Semi-Classical Quantum Defect Radial Integral */
96 /* */
97 /* See for example */
98 /* Atomic, Molecular & Optical Physics Handbook */
99 /* Gordon W. F. Drake; Editor */
100 /* AIP Press */
101 /* Woddbury, New York. */
102 /* 1996 */
103 /* */
104 /* NOTE:: we do not include the Bohr Radius a_o in the */
105 /* definition of of R(n,L;n'L') as per Drake. */
106 /* */
107 /* */
108 /* 1 (n_c)^2 | { D_l max(L,L') } */
109 /* R(n,L;n'L') = --- ------- | { 1 - ------------- } J( D_n-1; -x ) - */
110 /* Z 2 D_n | { n_c } */
111 /* */
112 /* */
113 /* { D_L max(L,L') } */
114 /* - { 1 + ------------- } J( D_n+1; -x ) */
115 /* { n_c } */
116 /* */
117 /* */
118 /* 2 | */
119 /* + --- sin(Pi D_n)(1-e) | */
120 /* Pi | */
121 /* */
122 /* where */
123 /* n_c = (2n*'n*)/(n*'+n*) */
124 /* */
125 /* Here is the quantity Drake gives... */
126 /* n_c = ( 2.0 * nstar * npstar ) / ( nstar + npstar ); */
127 /* */
128 /* while V.A. Davidkin uses */
129 /* n_c = sqrt( nstar * npstar ); */
130 /* */
131 /* D_n = n*' - n* */
132 /* */
133 /* D_L = L*' - L* */
134 /* */
135 /* x = e D_n */
136 /* */
137 /* Lmx = max(L',L) */
138 /* */
139 /* e = sqrt( 1 - {Lmx/n_c}^2 ) */
140 /* */
141 /* */
142 /* Here n* = n - qd where qd is the quantum defect */
143 /* */
144 /******************************************************************************/
145 /******************************************************************************/
146 double scqdri(/* upper and lower quantum numbers...n's are effective */
147  double nstar, long int l,
148  double npstar, long int lp,
149  double iz
150  )
151 {
152  double n_c = ((2.0 * nstar * npstar ) / ( nstar + npstar ));
153 
154  double D_n = (nstar - npstar);
155  double D_l = (double) ( l - lp );
156  double lg = (double) ( (lp > l) ? lp : l);
157 
158  double h = (lg/n_c);
159  double g = h*h;
160  double f = ( 1.0 - g );
161  double e = (( f >= 0.0) ? sqrt( f ) : 0.0 );
162 
163  double x = (e * D_n);
164  double z = (-1.0 * x);
165  double v1 = (D_n + 1.0);
166  double v2 = (D_n - 1.0);
167 
168  double d1,d2,d7,d8,d9,d34,d56,d6_1;
169 
170  DEBUG_ENTRY( "scqdri()" );
171 
172  if( iz == 0.0 )
173  iz += 1.0;
174 
175  if( D_n == 0.0 )
176  {
177  return( -1.0 );
178  }
179 
180  if( D_n < 0.0 )
181  {
182  return( -1.0 );
183  }
184 
185  if( f < 0.0 )
186  {
187  /* This can happen for certain quantum defects */
188  /* in the lower n=1:l=0 state. In which case you */
189  /* probably should be using some other alogrithm */
190  /* or theory to calculate the dipole moment. */
191  return( -1.0 );
192  }
193 
194  d1 = ( 1.0 / iz );
195 
196  d2 = (n_c * n_c)/(2.0 * D_n);
197 
198  d34 = (1.0 - ((D_l * lg)/n_c)) * AngerJ( v1, z );
199 
200  d56 = (1.0 + ((D_l * lg)/n_c)) * AngerJ( v2, z );
201 
202  d6_1 = PI * D_n;
203 
204  d7 = (2./PI) * sin( d6_1 ) * (1.0 - e);
205 
206  d8 = d1 * d2 * ( (d34) - (d56) + d7 );
207 
208  d9 = d8 * d8;
209 
210  ASSERT( D_n > 0.0 );
211  ASSERT( l >= 0 );
212  ASSERT( lp >= 0 );
213  ASSERT( (l == lp + 1) || ( l == lp - 1) );
214  ASSERT( n_c != 0.0 );
215  ASSERT( f >= 0.0 );
216  ASSERT( d9 > 0.0 );
217 
218  return( d9 );
219 }
220 
221 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem )
222 {
223  double A;
224  /* >>refer Helike 2pho Derevianko, A., & Johnson, W.R. 1997, Phys. Rev. A 56, 1288
225  * numbers are not explicitly given in this paper for Z=21-24,26,27,and 29.
226  * So numbers given here are interpolated. */
227  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,
228  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,
229  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};
230  /* Important clarification, according to Derevianko & Johnson (see ref above), 2^3S can decay
231  * to ground in one of two ways: through a two-photon process, or through a single-photon M1 decay,
232  * but the M1 rates are about 10^4 greater that the two-photon decays throughout the entire
233  * sequence. Thus these numbers, are much weaker than the effective decay rate, but should probably
234  * be treated in as a two-photon decay at some point */
235  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,
236  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,
237  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};
238 
239  DEBUG_ENTRY( "ForbiddenAuls()" );
240 
241  int ipISO = ipHE_LIKE;
242 
243  if( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
244  {
245  if( nelem == ipHELIUM )
246  {
247  /* All of these but the second and third one (values 51.02 and 1E-20) are from
248  * >>refer HeI As Lach, G, & Pachucki, K, 2001, Phys. Rev. A 64, 042510
249  ,* 1E-20 is made up
250  * 51.3 is from the Derevianko & Johnson paper cited above. */
251  double ForbiddenHe[5] = { 1.272E-4, 51.02, 1E-20, 177.58, 0.327 };
252 
253  fixit(); // adding the 2-photon decay 2^3S - 1^1S may be important in early universe
254  A = ForbiddenHe[ipHi - 1];
255  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
256  }
257  else
258  {
259  switch ( (int)ipHi )
260  {
261  case 1: /* Parameters for 2^3S to ground transition. */
262  /* >>refer Helike As Lin, C.D., Johnson, W.R., and Dalgarno, A. 1977,
263  * >>refercon Phys. Rev. A 15, 1, 010015 */
264  A = (3.9061E-7) * pow( (double)nelem+1., 10.419 ) + As2nuFrom3S[nelem-2];
265  break;
266  case 2: /* Parameters for 2^1S to ground transition. */
267  A = As2nuFrom1S[nelem-2];
268  break;
269  case 3: /* Parameters for 2^3P0 to ground transition. */
270  A = iso.SmallA;
271  break;
272  case 4: /* Parameters for 2^3P1 to ground transition. */
273  /* >>chng 06 jul 25, only use the fit to Johnson et al. values up to and
274  * including argon, where there calculation stops. For higher Z use below */
275  if( nelem <= ipARGON )
276  {
277  A = ( 11.431 * pow((double)nelem, 9.091) );
278  }
279  else
280  {
281  /* a fit to the Lin et al. 1977. values, which go past zinc. */
282  A = ( 383.42 * pow((double)nelem, 7.8901) );
283  }
284  break;
285  case 5: /* Parameters for 2^3P2 to ground transition. */
286  /* fit to Lin et al. 1977 values. This fit is good
287  * to 7% for the range from carbon to iron. The Lin et al. values
288  * differ from the Hata and Grant (1984) values (only given from
289  * calcium to copper) by less than 2%. */
290  A = ( 0.11012 * pow((double)nelem, 7.6954) );
291  break;
292  default:
293  TotalInsanity();
294  }
295  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
296  }
297  /* For now just just put 1% error for forbidden lines. */
298  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,.01f);
299  return A;
300  }
301 
302  /* The next two cases are fits to probabilities given in
303  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
304  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
305  /* originally astro.ph. 0201454 */
306  /* The involve Triplet P and Singlet S. Rates for Triplet S to Singlet P
307  * do not seem to be available. */
308 
309  /* Triplet P to Singlet S...Delta n not equal to zero! */
310  else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==3 &&
311  L_(ipLo)==0 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
312  {
313  A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
314  pow((double)N_(ipHi),2.877);
315  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
316  }
317 
318  /* Singlet S to Triplet P...Delta n not equal to zero! */
319  else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==1 &&
320  L_(ipLo)==1 && S_(ipLo)==3 && N_(ipLo) < N_(ipHi) )
321  {
322  A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
323 
324  if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
325  {
326  /* This is divided by 3 instead of 9, because value calculated is specifically for 2^3P1.
327  * Here we assume statistical population of the other two. */
328  A *= (2.*(ipLo-3)+1.0)/3.0;
329  }
330  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
331  }
332 
333  else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe3s3S ) )
334  {
335  double As_3TripS_to_2TripS[9] = {
336  7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
337  1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
338 
339  /* These M1 transitions given by
340  * >>refer He-like As Savukov, I.M., Labzowsky, and Johnson, W.R. 2005, PhysRevA, 72, 012504 */
341  if( nelem <= ipNEON )
342  {
343  A = As_3TripS_to_2TripS[nelem-1];
344  /* 20% error is based on difference between Savukov, Labzowsky, and Johnson (2005)
345  * and Lach and Pachucki (2001) for the helium transition. */
346  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f);
347  }
348  else
349  {
350  /* This is an extrapolation to the data given above. The fit reproduces
351  * the above values to 10% or better. */
352  A= 7.22E-9*pow((double)nelem, 9.33);
353  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f);
354  }
355  }
356 
357  else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
358  {
359  /* This transition,1.549 , given by Lach and Pachucki, 2001 for the atom */
360  if( nelem == ipHELIUM )
361  {
362  A = 1.549;
363  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
364  }
365  else
366  {
367  /* This is a fit to data given in
368  * >>refer He-like As Savukov, I.M., Johnson, W.R., & Safronova, U.I.
369  * >>refercon astro-ph 0205163 */
370  A= 0.1834*pow((double)nelem, 6.5735);
371  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
372  }
373  }
374 
375  else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 )
376  {
377  /* >>refer He As Bulatov, N.N. 1976, Soviet Astronomy, 20, 315 */
378  fixit();
379  /* This is the 29.6 GHz line that can be seen in radio galaxies. */
384  A = 5.78E-12;
385  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
386 
387  }
388 
389  else if( nelem==ipHELIUM && ipHi==5 && ipLo==4 )
390  {
391  fixit();
392  /* This is the 3 GHz line that can be seen in radio galaxies. */
397  A = 3.61E-15;
398  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
399  }
400 
401  else
402  {
403  /* Current transition is not supported. */
404  A = iso.SmallA;
405  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
406  }
407 
408  /* For now just put 1% error for forbidden lines. */
409  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,.01f);
410 
411  ASSERT( A > 0.);
412 
413  return A;
414 }
415 
416 /* Calculates Einstein A for a given transition. */
417 double he_1trans(
418  /* charge on c scale, Energy is wavenumbers, Einstein A */
419  long nelem , double Enerwn ,
420  /* quantum numbers of upper level: */
421  double Eff_nupper, long lHi, long sHi, long jHi,
422  /* and of lower level: */
423  double Eff_nlower, long lLo, long sLo, long jLo )
424  /* Note j is only necessary for 2 triplet P...for all other n,l,s,
425  * j is completely ignored. */
426 {
427  int ipISO = ipHE_LIKE;
428  double RI2, Aul;
429  long nHi, nLo, ipHi, ipLo;
430 
431  DEBUG_ENTRY( "he_1trans()" );
432 
433  ASSERT(nelem > ipHYDROGEN);
434 
435  /* Since 0.4 is bigger than any defect, adding that to the effective principle quantum number,
436  * and truncating to an integer will produce the principal quantum number. */
437  nHi = (int)(Eff_nupper + 0.4);
438  nLo = (int)(Eff_nlower + 0.4);
439 
440  /* Make sure this worked correctly. */
441  ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
442  ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );
443 
444  ipHi = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nHi][lHi][sHi];
445  if( (nHi==2) && (lHi==1) && (sHi==3) )
446  {
447  ASSERT( (jHi>=0) && (jHi<=2) );
448  ipHi -= (2 - jHi);
449  }
450 
451  ipLo = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nLo][lLo][sLo];
452  if( (nLo==2) && (lLo==1) && (sLo==3) )
453  {
454  ASSERT( (jLo>=0) && (jLo<=2) );
455  ipLo -= (2 - jLo);
456  }
457 
458  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n == nHi );
459  if( nHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
460  {
461  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].l == lHi );
462  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].S == sHi );
463  }
464  ASSERT( StatesElem[ipHE_LIKE][nelem][ipLo].n == nLo );
465  if( nLo <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
466  {
467  ASSERT( StatesElem[ipHE_LIKE][nelem][ipLo].l == lLo );
468  ASSERT( StatesElem[ipHE_LIKE][nelem][ipLo].S == sLo );
469  }
470 
471  /* First do allowed transitions */
472  if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
473  {
474  Aul = -2.;
475 
476  /* For clarity, let's do this in two separate chunks...one for helium, one for everything else. */
477  if( nelem == ipHELIUM )
478  {
479  /* Retrieve transition probabilities for Helium. */
480  /* >>refer He As Drake, G.W.F., Atomic, Molecular, and Optical Physics Handbook */
481  if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
482  {
483  /*Must not be accessed by collapsed levels! */
484  ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ) );
485  ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ) );
486  ASSERT( ipHi > 2 );
487 
488  Aul = TransProbs[nelem][ipHi][ipLo];
489 
490  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0005f);
491  }
492 
493  if( Aul < 0. )
494  {
495  /* Here are the Lyman transitions. */
496  if( ipLo == ipHe1s1S )
497  {
498  ASSERT( (lHi == 1) && (sHi == 1) );
499 
500  /* these fits calculated from Drake A's (1996) */
501  if( nLo == 1 )
502  Aul = (1.59208e10) / pow(Eff_nupper,3.0);
503  ASSERT( Aul > 0.);
504  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.005f);
505  }
506 
507  /* last resort for transitions involving significant defects,
508  * except that highest lLo are excluded */
509  else if( lHi>=2 && lLo>=2 && nHi>nLo )
510  {
511  Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
512  ASSERT( Aul > 0.);
513 
514  if( lHi + lLo >= 7 )
515  {
516  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f);
517  }
518  else
519  {
520  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f);
521  }
522  }
523  /* These fits come from extrapolations of Drake's oscillator strengths
524  * out to the series limit. We also use this method to obtain threshold
525  * photoionization cross-sections for the lower level of each transition here. */
526  /* See these two references :
527  * >>refer He As Hummer, D. G. \& Storey, P. J. 1998, MNRAS 297, 1073
528  * >>refer Seaton's Theorem Seaton, M. J. 1958, MNRAS 118, 504 */
529  else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
530  {
531  int paramSet=0;
532  double emisOscStr, x, a, b, c;
533  double extrapol_Params[2][4][4][3] = {
534  /* these are for singlets */
535  {
536  { /* these are P to S */
537  { 0.8267396 , 1.4837624 , -0.4615955 },
538  { 1.2738405 , 1.5841806 , -0.3022984 },
539  { 1.6128996 , 1.6842538 , -0.2393057 },
540  { 1.8855491 , 1.7709125 , -0.2115213 }},
541  { /* these are S to P */
542  { -1.4293664 , 2.3294080 , -0.0890470 },
543  { -0.3608082 , 2.3337636 , -0.0712380 },
544  { 0.3027974 , 2.3326252 , -0.0579008 },
545  { 0.7841193 , 2.3320138 , -0.0497094 }},
546  { /* these are D to P */
547  { 1.1341403 , 3.1702435 , -0.2085843 },
548  { 1.7915926 , 2.4942946 , -0.2266493 },
549  { 2.1979400 , 2.2785377 , -0.1518743 },
550  { 2.5018229 , 2.1925720 , -0.1081966 }},
551  { /* these are P to D */
552  { 0.0000000 , 0.0000000 , 0.0000000 },
553  { -2.6737396 , 2.9379143 , -0.3805367 },
554  { -1.4380124 , 2.7756396 , -0.2754625 },
555  { -0.6630196 , 2.6887253 , -0.2216493 }},
556  },
557  /* these are for triplets */
558  {
559  { /* these are P to S */
560  { 0.3075287 , 0.9087130 , -1.0387207 },
561  { 0.687069 , 1.1485864 , -0.6627317 },
562  { 0.9776064 , 1.3382024 , -0.5331906 },
563  { 1.2107725 , 1.4943721 , -0.4779232 }},
564  { /* these are S to P */
565  { -1.3659605 , 2.3262253 , -0.0306439 },
566  { -0.2899490 , 2.3279391 , -0.0298695 },
567  { 0.3678878 , 2.3266603 , -0.0240021 },
568  { 0.8427457 , 2.3249540 , -0.0194091 }},
569  { /* these are D to P */
570  { 1.3108281 , 2.8446367 , -0.1649923 },
571  { 1.8437692 , 2.2399326 , -0.2583398 },
572  { 2.1820792 , 2.0693762 , -0.1864091 },
573  { 2.4414052 , 2.0168255 , -0.1426083 }},
574  { /* these are P to D */
575  { 0.0000000 , 0.0000000 , 0.0000000 },
576  { -1.9219877 , 2.7689624 , -0.2536072 },
577  { -0.7818065 , 2.6595150 , -0.1895313 },
578  { -0.0665624 , 2.5955623 , -0.1522616 }},
579  }
580  };
581 
582  if( lLo==0 )
583  {
584  paramSet = 0;
585  }
586  else if( lLo==1 && lHi==0 )
587  {
588  paramSet = 1;
589  }
590  else if( lLo==1 && lHi==2 )
591  {
592  paramSet = 2;
593  }
594  else if( lLo==2 )
595  {
596  paramSet = 3;
597  ASSERT( lHi==1 );
598  }
599 
600  ASSERT( (int)((sHi-1)/2) == 0 || (int)((sHi-1)/2) == 1 );
601  a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
602  b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
603  c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
604  ASSERT( Enerwn > 0. );
605  x = log( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]*RYD_INF/Enerwn );
606 
607  emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
608  (2.*lLo+1)/(2.*lHi+1);
609 
610  Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;
611 
612  if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
613  {
614  Aul *= (2.*(ipLo-3)+1.0)/9.0;
615  }
616 
617  ASSERT( Aul > 0. );
618  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f);
619  }
620  else
621  {
622  /* Calculate the radial integral from the quantum defects. */
623  RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
624  ASSERT( RI2 > 0. );
625  /* Convert radial integral to Aul. */
626  Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
627  /* radial integral routine does not recognize fine structure.
628  * Here we split 2^3P. */
629  if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
630  {
631  Aul *= (2.*(ipLo-3)+1.0)/9.0;
632  }
633 
634  ASSERT( Aul > 0. );
635  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.03f);
636  }
637  }
638  }
639 
640  /* Heavier species */
641  else
642  {
643  /* Retrieve transition probabilities for Helike ions. */
644  /* >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
645  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J, originally astro.ph. 0201454 */
646  if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
647  {
648  /*Must not be accessed by collapsed levels! */
649  ASSERT( ipHi < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) );
650  ASSERT( ipLo < ( iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) );
651  ASSERT( ipHi > 2 );
652 
653  Aul = TransProbs[nelem][ipHi][ipLo];
654  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
655  }
656 
657  if( Aul < 0. )
658  {
659  /* Do same-n transitions. */
660  if( nLo==nHi && lHi<=2 && lLo<=2 )
661  {
662  /* These are 2p3Pj to 2s3S fits to (low Z) Porquet & Dubau (2000) &
663  * (high Z) NIST Atomic Spectra Database. */
664  if( ipLo == ipHe2s3S )
665  {
666  if(ipHi == ipHe2p3P0)
667  Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
668  else if(ipHi == ipHe2p3P1)
669  Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
670  else if(ipHi == ipHe2p3P2)
671  Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
672  else
673  {
674  fprintf(ioQQQ," PROBLEM Impossible value for ipHi in he_1trans.\n");
675  TotalInsanity();
676  }
677  }
678 
679  /* These are 2p1P to 2s1S fits to data from TOPbase. */
680  else if( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) )
681  {
682  Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
683  }
684 
685  else
686  {
687  /* This case should only be entered if n > 2. Those cases were done above. */
688  ASSERT( nLo > 2 );
689 
690  /* Triplet P to triplet S, delta n = 0 */
691  if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
692  {
693  Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
694  }
695  /* Singlet P to singlet D, delta n = 0 */
696  else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
697  {
698  Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
699  }
700  /* Singlet P to singlet S, delta n = 0 */
701  else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
702  {
703  Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077);
704  }
705  else
706  {
707  ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
708  Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
709  if( (lHi >2) || (lLo > 2) )
710  Aul *= (lHi/2.);
711  if(lLo > 2)
712  Aul *= (1./9.);
713  }
714  }
715  ASSERT( Aul > 0.);
716  }
717 
718  /* assume transitions involving F and higher orbitals are hydrogenic. */
719  else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
720  {
721  Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
722  ASSERT( Aul > 0.);
723  }
724 
725  /* These transitions are of great importance, but the below radial integral
726  * routine fails to achieve desirable accuracy, so these are fits as produced
727  * from He A's for nupper through 9. They are transitions to ground and
728  * 2, 3, and 4 triplet S. */
729  else if( ( ipLo == 0 ) || ( ipLo == ipHe2s1S ) || ( ipLo == ipHe2s3S )
730  || ( ipLo == ipHe3s3S ) || ( ipLo == ipHe4s3S ) )
731  {
732  /* Here are the Lyman transitions. */
733  if( ipLo == 0 )
734  {
735  ASSERT( (lHi == 1) && (sHi == 1) );
736 
737  /* In theory, this Z dependence should be Z^4, but values from TOPbase
738  * suggest 3.9 is a more accurate exponent. Values from
739  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
740  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
741  /* originally astro.ph. 0201454 */
742  Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
743  }
744 
745  /* Here are the Balmer transitions. */
746  else if( ipLo == ipHe2s1S )
747  {
748  ASSERT( (lHi == 1) && (sHi == 1) );
749 
750  /* More fits, as above. */
751  Aul = 5.0e8 * pow((double)nelem,4.) / pow((double)nHi, 2.889);
752  }
753 
754  /* Here are transitions down to triplet S */
755  else
756  {
757  ASSERT( (lHi == 1) && (sHi == 3) );
758 
759  /* These fits also as above. */
760  if( nLo == 2 )
761  Aul = 1.5 * 3.405E8 * pow((double)nelem,4.) / pow((double)nHi, 2.883);
762  else if( nLo == 3 )
763  Aul = 2.5 * 4.613E7 * pow((double)nelem,4.) / pow((double)nHi, 2.672);
764  else
765  Aul = 3.0 * 1.436E7 * pow((double)nelem,4.) / pow((double)nHi, 2.617);
766  }
767 
768  ASSERT( Aul > 0.);
769  }
770 
771  /* Every other allowed transition is calculated as follows. */
772  else
773  {
774  /* Calculate the radial integral from the quantum defects. */
775  RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem));
776  /* Convert radial integral to Aul. */
777  Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2);
778  /* radial integral routine does not recognize fine structure.
779  * Here we split 2^3P. */
780  if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) && (Aul > iso.SmallA) )
781  {
782  Aul *= (2.*(ipLo-3)+1.0)/9.0;
783  }
784 
785  }
786  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f);
787  }
788 
789  /* for now just give ions some a 5% error across the board */
790  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.05f);
791  }
792  }
793 
794  /* Now do forbidden transitions from 2-1 ... */
795  /* and those given by
796  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
797  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
798  /* originally astro.ph. 0201454
799  * for heavy elements. These are triplet P to singlet S,
800  * going either up or down...Triplet S to Singlet P are not included, as they are far weaker. */
801  else
802  {
803  ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
804  Aul = ForbiddenAuls(ipHi, ipLo, nelem);
805  ASSERT( Aul > 0. );
806  }
807 
808  Aul = MAX2( Aul, iso.SmallA );
809  ASSERT( Aul >= iso.SmallA );
810 
811  /* negative energy for a transition with substantial transition probability
812  * would be major logical error - but is ok for same-n l transitions */
813  if( Enerwn < 0. && Aul > iso.SmallA )
814  {
815  fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
816  }
817 
818  return Aul;
819 }
820 
821 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
822 {
823  /* Every bit of this routine is based upon the singlet-triplet mixing formalism given in
824  * >>refer He FSM Drake, G. W. F. 1996, in Atomic, Molecular, \& Optical Physics Handbook,
825  * >>refercon ed. G. W. F. Drake (New York: AIP Press).
826  * That formalism mixes the levels themselves, but since this code is not J-resolved, we simulate
827  * that by mixing only the transition probabilities. We find results comparable to those calculated
828  * in the fully J-resolved model spearheaded by Rob Bauman, and described in
829  * >>refer He FSM Bauman, R. P., Porter, R. L., Ferland, G. J., \& MacAdam, K. B. 2005, ApJ, accepted */
830  long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
831  double Ass, Att, Ast, Ats;
832  double SinHi, SinLo, CosHi, CosLo;
833  double HiMixingAngle, LoMixingAngle , error;
834  double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
835 
836  DEBUG_ENTRY( "DoFSMixing()" );
837 
838  nHi = StatesElem[ipHE_LIKE][nelem][ipHiSing].n;
839  lHi = StatesElem[ipHE_LIKE][nelem][ipHiSing].l;
840  sHi = StatesElem[ipHE_LIKE][nelem][ipHiSing].S;
841  nLo = StatesElem[ipHE_LIKE][nelem][ipLoSing].n;
842  lLo = StatesElem[ipHE_LIKE][nelem][ipLoSing].l;
843  sLo = StatesElem[ipHE_LIKE][nelem][ipLoSing].S;
844 
845  if( ( sHi == 3 || sLo == 3 ) ||
846  ( abs(lHi - lLo) != 1 ) ||
847  ( nLo < 2 ) ||
848  ( lHi <= 1 || lLo <= 1 ) ||
849  ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
850  ( nHi > nLo && lHi != 1 && lLo != 1 ) )
851  {
852  return;
853  }
854 
855  ASSERT( lHi > 0 );
856  /*ASSERT( (lHi > 1) && (lLo > 1) );*/
857 
858  ipHiTrip = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nHi][lHi][3];
859  ipLoTrip = iso.QuantumNumbers2Index[ipHE_LIKE][nelem][nLo][lLo][3];
860 
861  if( lHi == 2 )
862  {
863  HiMixingAngle = 0.01;
864  }
865  else if( lHi == 3 )
866  {
867  HiMixingAngle = 0.5;
868  }
869  else
870  {
871  HiMixingAngle = PI/4.;
872  }
873 
874  if( lLo == 2 )
875  {
876  LoMixingAngle = 0.01;
877  }
878  else if( lLo == 3 )
879  {
880  LoMixingAngle = 0.5;
881  }
882  else
883  {
884  LoMixingAngle = PI/4.;
885  }
886 
887  /* These would not work correctly if l<=1 were included in this treatment! */
888  ASSERT( ipHiTrip > ipLoTrip );
889  ASSERT( ipHiTrip > ipLoSing );
890  ASSERT( ipHiSing > ipLoTrip );
891  ASSERT( ipHiSing > ipLoSing );
892 
893  SinHi = sin( HiMixingAngle );
894  SinLo = sin( LoMixingAngle );
895  CosHi = cos( HiMixingAngle );
896  CosLo = cos( LoMixingAngle );
897 
898  Kss = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].EnergyWN;
899  Ktt = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].EnergyWN;
900  Kst = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].EnergyWN;
901  Kts = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].EnergyWN;
902 
903  fss = Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul/TRANS_PROB_CONST/POW2(Kss)/
904  Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Lo->g*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Hi->g;
905 
906  ftt = Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul/TRANS_PROB_CONST/POW2(Ktt)/
907  Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Lo->g*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Hi->g;
908 
909  temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
910  fssNew = Kss*POW2( temp );
911  temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
912  fttNew = Ktt*POW2( temp );
913  temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
914  fstNew = Kst*POW2( temp );
915  temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
916  ftsNew = Kts*POW2( temp );
917 
918  Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Lo->g/
919  Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Hi->g;
920 
921  Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Lo->g/
922  Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Hi->g;
923 
924  Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Lo->g/
925  Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Hi->g;
926 
927  Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Lo->g/
928  Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Hi->g;
929 
930  error = fabs( ( Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul+
931  Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul )/
932  (Ass+Ast+Ats+Att) - 1.f );
933 
934  if( error > 0.001 )
935  {
936  fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
937  ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
938  Ass/Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul,
939  Att/Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul,
940  Ast/Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Emis->Aul,
941  Ats/Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Emis->Aul );
942  }
943 
944  Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoSing].Emis->Aul = (realnum)Ass;
945  Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoTrip].Emis->Aul = (realnum)Att;
946  Transitions[ipHE_LIKE][nelem][ipHiSing][ipLoTrip].Emis->Aul = (realnum)Ast;
947  Transitions[ipHE_LIKE][nelem][ipHiTrip][ipLoSing].Emis->Aul = (realnum)Ats;
948 
949  return;
950 }
951 
952 /*ritoa converts the square of the radial integral for a transition
953  * (calculated by scqdri) to the transition probability, Aul. */
954 STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
955 {
956  /* Variables are as follows: */
957  /* lg = larger of li and lf */
958  /* fmean = mean oscillator strength */
959  /* for a given level. */
960  /* mu = reduced mass of optical electron. */
961  /* EinsteinA = Einstein emission coef. */
962  /* w = angular frequency of transition. */
963  /* RI2_cm = square of rad. int. in cm^2. */
964  long lg;
965  double fmean,mu,EinsteinA,w,RI2_cm;
966 
967  DEBUG_ENTRY( "ritoa()" );
968 
970 
971  w = 2. * PI * k * SPEEDLIGHT;
972 
973  RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM;
974 
975  lg = (lf > li) ? lf : li;
976 
977  fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));
978 
979  EinsteinA = TRANS_PROB_CONST*k*k*fmean;
980 
981  /* ASSERT( EinsteinA > SMALLFLOAT ); */
982 
983  return EinsteinA;
984 }
985 
986 realnum helike_transprob( long nelem, long ipHi, long ipLo )
987 {
988  double Aul, Aul1;
989  double Enerwn, n_eff_hi, n_eff_lo;
990  long ipISO = ipHE_LIKE;
991  /* charge to 4th power, needed for scaling laws for As*/
992  double z4 = POW2((double)nelem);
993  z4 *= z4;
994 
995  DEBUG_ENTRY( "helike_transprob()" );
996 
997  Enerwn = Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN;
998  n_eff_hi = N_(ipHi) - helike_quantum_defect(nelem,ipHi);
999  n_eff_lo = N_(ipLo) - helike_quantum_defect(nelem,ipLo);
1000 
1001  if( ipHi >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] )
1002  {
1003  if( ipLo >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] )
1004  {
1005  /* Neither upper nor lower is resolved...Aul's are easy. */
1006  Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
1007  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f);
1008 
1009  ASSERT( Aul > 0.);
1010  }
1011  else
1012  {
1013  /* Lower level resolved, upper not. First calculate Aul
1014  * from upper level with ang mom one higher. */
1015  Aul = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)+1,
1016  S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
1017 
1018  iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][0] = (realnum)Aul;
1019 
1020  Aul *= (2.*L_(ipLo)+3.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
1021 
1022  if( L_(ipLo) != 0 )
1023  {
1024  /* For all l>0, add in transitions from upper level with ang mom one lower. */
1025  Aul1 = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)-1,
1026  S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
1027 
1028  iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = (realnum)Aul1;
1029 
1030  /* now add in this part after multiplying by stat weight for lHi = lLo-1. */
1031  Aul += Aul1*(2.*L_(ipLo)-1.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
1032  }
1033  else
1034  iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = 0.f;
1035 
1036  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f);
1037  ASSERT( Aul > 0.);
1038  }
1039  }
1040  else
1041  {
1042  /* Both levels are resolved...do the normal bit. */
1043  if( Enerwn < 0. )
1044  {
1045  Aul = he_1trans( nelem, -1.*Enerwn, n_eff_lo,
1046  L_(ipLo), S_(ipLo), ipLo-3, n_eff_hi, L_(ipHi), S_(ipHi), ipHi-3);
1047  }
1048  else
1049  {
1050  Aul = he_1trans( nelem, Enerwn, n_eff_hi,
1051  L_(ipHi), S_(ipHi), ipHi-3, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
1052  }
1053  }
1054 
1055  return (realnum)Aul;
1056 }
1057 
1058 /* This routine is pure infrastructure and bookkeeping, no physics. */
1060 {
1061 
1062  const int chLine_LENGTH = 1000;
1063  char chLine[chLine_LENGTH];
1064 
1065  FILE *ioDATA;
1066  bool lgEOL;
1067 
1068  long nelem, ipLo, ipHi, i, i1, i2, i3;
1069 
1070  DEBUG_ENTRY( "HelikeTransProbSetup()" );
1071 
1072  TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
1073 
1074  for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
1075  {
1076 
1077  TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) );
1078 
1079  for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo )
1080  {
1081  TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX );
1082  }
1083  }
1084 
1085  /********************************************************************/
1086  /*************** Read in data from he_transprob.dat *****************/
1087  if( trace.lgTrace )
1088  fprintf( ioQQQ," HelikeTransProbSetup opening he_transprob.dat:");
1089 
1090  ioDATA = open_data( "he_transprob.dat", "r" );
1091 
1092  /* check that magic number is ok */
1093  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1094  {
1095  fprintf( ioQQQ, " HelikeTransProbSetup could not read first line of he_transprob.dat.\n");
1096  cdEXIT(EXIT_FAILURE);
1097  }
1098  i = 1;
1099  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1100  i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1101  if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
1102  {
1103  fprintf( ioQQQ,
1104  " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1105  fprintf( ioQQQ,
1106  " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1107  TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
1108  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1109  cdEXIT(EXIT_FAILURE);
1110  }
1111 
1112  /* Initialize TransProbs[nelem][ipLo][ipHi] before filling it in. */
1113  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
1114  {
1115  for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
1116  {
1117  for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
1118  {
1119  TransProbs[nelem][ipHi][ipLo] = -1.;
1120  }
1121  }
1122  }
1123 
1124  for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
1125  {
1126  char *chTemp;
1127 
1128  /* get next line image */
1129  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1130  BadRead();
1131 
1132  while( chLine[0]=='#' )
1133  {
1134  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1135  BadRead();
1136  }
1137 
1138  i3 = 1;
1139  i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
1140  i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
1141  /* check that these numbers are correct */
1142  if( i1<0 || i2<=i1 )
1143  {
1144  fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1145  cdEXIT(EXIT_FAILURE);
1146  }
1147 
1148  chTemp = chLine;
1149 
1150  /* skip over 2 tabs to start of data */
1151  for( i=0; i<1; ++i )
1152  {
1153  if( (chTemp = strchr( chTemp, '\t' )) == NULL )
1154  {
1155  fprintf( ioQQQ, " HelikeTransProbSetup could not init he_transprob\n" );
1156  cdEXIT(EXIT_FAILURE);
1157  }
1158  ++chTemp;
1159  }
1160 
1161  /* now read in the data */
1162  for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
1163  {
1164  if( (chTemp = strchr( chTemp, '\t' )) == NULL )
1165  {
1166  fprintf( ioQQQ, " HelikeTransProbSetup could not scan he_transprob\n" );
1167  cdEXIT(EXIT_FAILURE);
1168  }
1169  ++chTemp;
1170 
1171  sscanf( chTemp , "%le" , &TransProbs[nelem][i2][i1] );
1172 
1174  if( lgEOL )
1175  {
1176  fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1177  cdEXIT(EXIT_FAILURE);
1178  }
1179  }
1180  }
1181 
1182  /* check that ending magic number is ok */
1183  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1184  {
1185  fprintf( ioQQQ, " HelikeTransProbSetup could not read last line of he_transprob.dat.\n");
1186  cdEXIT(EXIT_FAILURE);
1187  }
1188  i = 1;
1189  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1190  i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
1191  if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
1192  {
1193  fprintf( ioQQQ,
1194  " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1195  fprintf( ioQQQ,
1196  " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1197  TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
1198  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1199  cdEXIT(EXIT_FAILURE);
1200  }
1201 
1202  /* close the data file */
1203  fclose( ioDATA );
1204  return;
1205 }

Generated for cloudy by doxygen 1.8.3.1