cloudy  trunk
helike_einsta.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 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 
225  DEBUG_ENTRY( "ForbiddenAuls()" );
226 
227  int ipISO = ipHE_LIKE;
228 
229  if( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
230  {
231  if( nelem == ipHELIUM )
232  {
233  // All of these but the second and third one (with values 1E-20) are from
234  // >>refer HeI As Lach, G, & Pachucki, K, 2001, Phys. Rev. A 64, 042510
235  // 1E-20 is made up
236  double ForbiddenHe[5] = { 1.272E-4, 1E-20, 1E-20, 177.58, 0.327 };
237 
238  A = ForbiddenHe[ipHi - 1];
239  iso_put_error(ipHE_LIKE,nelem,ipHe2s3S ,ipHe1s1S,IPRAD, 0.01f, 0.20f);
240  iso_put_error(ipHE_LIKE,nelem,ipHe2s1S ,ipHe1s1S,IPRAD, 0.01f, 0.05f);
241  iso_put_error(ipHE_LIKE,nelem,ipHe2p3P0,ipHe1s1S,IPRAD, 0.00f, 0.00f);
242  iso_put_error(ipHE_LIKE,nelem,ipHe2p3P1,ipHe1s1S,IPRAD, 0.01f, 0.05f);
243  iso_put_error(ipHE_LIKE,nelem,ipHe2p3P2,ipHe1s1S,IPRAD, 0.01f, 0.01f);
244  }
245  else
246  {
247  switch ( (int)ipHi )
248  {
249  case 1: /* Parameters for 2^3S to ground transition. */
250  /* >>refer Helike As Lin, C.D., Johnson, W.R., and Dalgarno, A. 1977,
251  * >>refercon Phys. Rev. A 15, 1, 010015 */
252  // 2nu is treated elsewhere
253  A = (3.9061E-7) * pow( (double)nelem+1., 10.419 );
254  break;
255  case 2: /* Parameters for 2^1S to ground transition. */
256  A = iso_ctrl.SmallA; // 2nu is treated elsewhere
257  break;
258  case 3: /* Parameters for 2^3P0 to ground transition. */
259  A = iso_ctrl.SmallA;
260  break;
261  case 4: /* Parameters for 2^3P1 to ground transition. */
262  /* >>chng 06 jul 25, only use the fit to Johnson et al. values up to and
263  * including argon, where there calculation stops. For higher Z use below */
264  if( nelem <= ipARGON )
265  {
266  A = ( 11.431 * pow((double)nelem, 9.091) );
267  }
268  else
269  {
270  /* a fit to the Lin et al. 1977. values, which go past zinc. */
271  A = ( 383.42 * pow((double)nelem, 7.8901) );
272  }
273  break;
274  case 5: /* Parameters for 2^3P2 to ground transition. */
275  /* fit to Lin et al. 1977 values. This fit is good
276  * to 7% for the range from carbon to iron. The Lin et al. values
277  * differ from the Hata and Grant (1984) values (only given from
278  * calcium to copper) by less than 2%. */
279  A = ( 0.11012 * pow((double)nelem, 7.6954) );
280  break;
281  default:
282  TotalInsanity();
283  }
284  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
285  }
286 
287  return A;
288  }
289 
290  /* The next two cases are fits to probabilities given in
291  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
292  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
293  /* originally astro.ph. 0201454 */
294  /* The involve Triplet P and Singlet S. Rates for Triplet S to Singlet P
295  * do not seem to be available. */
296 
297  /* Triplet P to Singlet S...Delta n not equal to zero! */
298  else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==3 &&
299  L_(ipLo)==0 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
300  {
301  A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
302  pow((double)N_(ipHi),2.877);
303  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
304  }
305 
306  /* Singlet S to Triplet P...Delta n not equal to zero! */
307  else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==1 &&
308  L_(ipLo)==1 && S_(ipLo)==3 && N_(ipLo) < N_(ipHi) )
309  {
310  A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
311 
312  if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
313  {
314  /* This is divided by 3 instead of 9, because value calculated is specifically for 2^3P1.
315  * Here we assume statistical population of the other two. */
316  A *= (2.*(ipLo-3)+1.0)/3.0;
317  }
318  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
319  }
320 
321  else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe3s3S ) )
322  {
323  double As_3TripS_to_2TripS[9] = {
324  7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
325  1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
326 
327  /* These M1 transitions given by
328  * >>refer He-like As Savukov, I.M., Labzowsky, and Johnson, W.R. 2005, PhysRevA, 72, 012504 */
329  if( nelem <= ipNEON )
330  {
331  A = As_3TripS_to_2TripS[nelem-1];
332  /* 20% error is based on difference between Savukov, Labzowsky, and Johnson (2005)
333  * and Lach and Pachucki (2001) for the helium transition. */
334  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f,0.2f);
335  }
336  else
337  {
338  /* This is an extrapolation to the data given above. The fit reproduces
339  * the above values to 10% or better. */
340  A= 7.22E-9*pow((double)nelem, 9.33);
341  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f,0.3f);
342  }
343  }
344 
345  else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
346  {
347  /* This transition,1.549 , given by Lach and Pachucki, 2001 for the atom */
348  if( nelem == ipHELIUM )
349  {
350  A = 1.549;
351  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
352  }
353  else
354  {
355  /* This is a fit to data given in
356  * >>refer He-like As Savukov, I.M., Johnson, W.R., & Safronova, U.I.
357  * >>refercon astro-ph 0205163 */
358  A= 0.1834*pow((double)nelem, 6.5735);
359  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
360  }
361  }
362 
363  else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 )
364  {
365  /* >>refer He As Bulatov, N.N. 1976, Soviet Astronomy, 20, 315 */
366  fixit();
367  /* This is the 29.6 GHz line that can be seen in radio galaxies. */
372  A = 5.78E-12;
373  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
374 
375  }
376 
377  else if( nelem==ipHELIUM && ipHi==5 && ipLo==4 )
378  {
379  fixit();
380  /* This is the 3 GHz line that can be seen in radio galaxies. */
385  A = 3.61E-15;
386  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
387  }
388 
389  else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe1s1S )
390  {
391  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
392  A = 44.326 * (1./3.) + 0.1146547 * (5./9.);
393  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
394  }
395 
396  else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe2s1S )
397  {
398  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
399  A = 1.459495 * (1./3.) + 3.6558e-5 * (5./9.);
400  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
401  }
402 
403  else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe3s1S )
404  {
405  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
406  A = 2.2297e-3 * (1./3.);
407  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
408  }
409 
410  else if( nelem==ipHELIUM && ipHi==ipHe3p1P && ipLo==ipHe2s3S )
411  {
412  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
413  A = 0.1815436;
414  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
415  }
416 
417  else if( nelem==ipHELIUM && ipHi==ipHe3p1P && ipLo==ipHe3s3S )
418  {
419  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
420  A = 0.14895852;
421  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
422  }
423 
424  else if( ( ipLo==0 && L_(ipHi)==2 && S_(ipHi)==1 ) ||
425  ( ipLo==ipHe2s3S && L_(ipHi)==2 && S_(ipHi)==3 ) ||
426  ( ipLo==ipHe2s1S && L_(ipHi)==2 && S_(ipHi)==1 ) ||
427  ( N_(ipLo)==2 && L_(ipLo)==1 && S_(ipLo)==3 && N_(ipHi)>=3 && L_(ipHi)==1 && S_(ipHi)==3 ) ||
428  ( ipLo==ipHe2p1P && L_(ipHi)==1 && S_(ipHi)==1 ) )
429  {
430  // >>refer Helike QOS Cann, N. M. \& Thakkar, A. J. 2002, J.Phys.B, 35, 421
431 
432  double f_params[5][4][3] = {
433  {
434  /* 1,0,3,2,1 */ {9.360591E-007, -3.1937E-006, 3.5186E-006},
435  /* 1,0,4,2,1 */ {4.631435E-007, -1.4973E-006, 1.4848E-006},
436  /* 1,0,5,2,1 */ {2.493912E-007, -7.8658E-007, 7.3994E-007},
437  /* 1,0,6,2,1 */ {1.476742E-007, -4.5953E-007, 4.1932E-007}},
438  {
439  /* 2,0,3,2,3 */ {1.646733E-006, -2.0028E-006, -1.3552E-006},
440  /* 2,0,4,2,3 */ {9.120593E-008, 3.1301E-007, -3.2190E-007},
441  /* 2,0,5,2,3 */ {1.360965E-008, 1.1467E-007, 8.6977E-008},
442  /* 2,0,6,2,3 */ {3.199421E-009, 4.5485E-008, 1.1016E-007}},
443  {
444  /* 2,0,3,2,1 */ {1.646733E-006, -2.9720E-006, 1.5367E-006},
445  /* 2,0,4,2,1 */ {9.120593E-008, -3.9037E-008, 3.9156E-008},
446  /* 2,0,5,2,1 */ {1.360965E-008, 1.4671E-008, 1.5626E-008},
447  /* 2,0,6,2,1 */ {3.199421E-009, 8.9806E-009, 1.2436E-008}},
448  {
449  /* 2,1,3,1,3 */ {1.543812E-007, -2.8220E-007, 3.0261E-008},
450  /* 2,1,4,1,3 */ {3.648237E-008, -6.6824E-008, 4.5879E-009},
451  /* 2,1,5,1,3 */ {1.488556E-008, -2.7304E-008, 1.6628E-009},
452  /* 2,1,6,1,3 */ {7.678610E-009, -1.4112E-008, 6.8453E-010}},
453  {
454  /* 2,1,3,1,1 */ {1.543812E-007, -2.8546E-007, 4.6605E-008},
455  /* 2,1,4,1,1 */ {3.648237E-008, -6.8422E-008, 1.7038E-008},
456  /* 2,1,5,1,1 */ {1.488556E-008, -2.8170E-008, 8.5012E-009},
457  /* 2,1,6,1,1 */ {7.678610E-009, -1.4578E-008, 4.6686E-009}}
458  };
459 
460  ASSERT( ipLo <= 6 );
461  long index_lo;
462  if( ipLo==ipHe2p1P )
463  index_lo = 4;
464  else if( ipLo==ipHe2p3P1 || ipLo==ipHe2p3P2 )
465  index_lo = 3;
466  else
467  index_lo = ipLo;
468 
469  ASSERT( N_(ipHi) >= 3 );
470  long index_hi = MIN2( N_(ipHi), 6 ) - 3;
471  double f_lu = POW2(nelem+1) * (
472  f_params[index_lo][index_hi][0] +
473  f_params[index_lo][index_hi][1]/(nelem+1) +
474  f_params[index_lo][index_hi][2]/POW2(nelem+1) );
475 
476  // extrapolate for higher upper n
477  if( N_(ipHi) > 6 )
478  f_lu *= pow( 6. / N_(ipHi), 3. );
479 
480  double gLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].g();
481  double gHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].g();
482  double eWN = iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN();
483 
484  A = eina( gLo * f_lu, eWN, gHi );
485  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
486  }
487 
488  else
489  {
490  /* Current transition is not supported. */
491  A = iso_ctrl.SmallA;
492  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
493  }
494 
495  ASSERT( A > 0.);
496 
497  return A;
498 }
499 
500 /* Calculates Einstein A for a given transition. */
501 double he_1trans(
502  /* charge on c scale, Energy is wavenumbers, Einstein A */
503  long nelem , double Enerwn ,
504  /* quantum numbers of upper level: */
505  double Eff_nupper, long lHi, long sHi, long jHi,
506  /* and of lower level: */
507  double Eff_nlower, long lLo, long sLo, long jLo )
508  /* Note j is only necessary for 2 triplet P...for all other n,l,s,
509  * j is completely ignored. */
510 {
511  int ipISO = ipHE_LIKE;
512  double RI2, Aul;
513  long nHi, nLo, ipHi, ipLo;
514 
515  DEBUG_ENTRY( "he_1trans()" );
516 
517  ASSERT(nelem > ipHYDROGEN);
518 
519  /* Since 0.4 is bigger than any defect, adding that to the effective principle quantum number,
520  * and truncating to an integer will produce the principal quantum number. */
521  nHi = (int)(Eff_nupper + 0.4);
522  nLo = (int)(Eff_nlower + 0.4);
523 
524  /* Make sure this worked correctly. */
525  ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
526  ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );
527 
528  ipHi = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][sHi];
529  if( (nHi==2) && (lHi==1) && (sHi==3) )
530  {
531  ASSERT( (jHi>=0) && (jHi<=2) );
532  ipHi -= (2 - jHi);
533  }
534 
535  ipLo = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][sLo];
536  if( (nLo==2) && (lLo==1) && (sLo==3) )
537  {
538  ASSERT( (jLo>=0) && (jLo<=2) );
539  ipLo -= (2 - jLo);
540  }
541 
542  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == nHi );
543  if( nHi <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
544  {
545  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].l() == lHi );
546  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == sHi );
547  }
548  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].n() == nLo );
549  if( nLo <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
550  {
551  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() == lLo );
552  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].S() == sLo );
553  }
554 
555  /* First do allowed transitions */
556  if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
557  {
558  Aul = -2.;
559 
560  /* For clarity, let's do this in two separate chunks...one for helium, one for everything else. */
561  if( nelem == ipHELIUM )
562  {
563  /* Retrieve transition probabilities for Helium. */
564  /* >>refer He As Drake, G.W.F., Atomic, Molecular, and Optical Physics Handbook */
565  if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max )
566  {
567  /*Must not be accessed by collapsed levels! */
568  ASSERT( ipHi < ( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max ) );
569  ASSERT( ipLo < ( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max ) );
570  ASSERT( ipHi > 2 );
571 
572  Aul = TransProbs[nelem][ipHi][ipLo];
573 
574  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0001f,0.002f);
575  }
576 
577  if( Aul < 0. )
578  {
579  /* Here are the Lyman transitions. */
580  if( ipLo == ipHe1s1S )
581  {
582  ASSERT( (lHi == 1) && (sHi == 1) );
583 
584  /* these fits calculated from Drake A's (1996) */
585  if( nLo == 1 )
586  Aul = (1.59208e10) / pow(Eff_nupper,3.0);
587  ASSERT( Aul > 0.);
588  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f);
589  }
590 
591  /* last resort for transitions involving significant defects,
592  * except that highest lLo are excluded */
593  else if( lHi>=2 && lLo>=2 && nHi>nLo )
594  {
595  Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
596  ASSERT( Aul > 0.);
597  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.006f,0.04f);
598  }
599  /* These fits come from extrapolations of Drake's oscillator strengths
600  * out to the series limit. We also use this method to obtain threshold
601  * photoionization cross-sections for the lower level of each transition here. */
602  /* See these two references :
603  * >>refer He As Hummer, D. G. \& Storey, P. J. 1998, MNRAS 297, 1073
604  * >>refer Seaton's Theorem Seaton, M. J. 1958, MNRAS 118, 504 */
605  else if( N_(ipHi)>10 && N_(ipLo)<=5 && lHi<=2 && lLo<=2 )
606  {
607  int paramSet=0;
608  double emisOscStr, x, a, b, c;
609  double extrapol_Params[2][4][4][3] = {
610  /* these are for singlets */
611  {
612  { /* these are P to S */
613  { 0.8267396 , 1.4837624 , -0.4615955 },
614  { 1.2738405 , 1.5841806 , -0.3022984 },
615  { 1.6128996 , 1.6842538 , -0.2393057 },
616  { 1.8855491 , 1.7709125 , -0.2115213 }},
617  { /* these are S to P */
618  { -1.4293664 , 2.3294080 , -0.0890470 },
619  { -0.3608082 , 2.3337636 , -0.0712380 },
620  { 0.3027974 , 2.3326252 , -0.0579008 },
621  { 0.7841193 , 2.3320138 , -0.0497094 }},
622  { /* these are D to P */
623  { 1.1341403 , 3.1702435 , -0.2085843 },
624  { 1.7915926 , 2.4942946 , -0.2266493 },
625  { 2.1979400 , 2.2785377 , -0.1518743 },
626  { 2.5018229 , 2.1925720 , -0.1081966 }},
627  { /* these are P to D */
628  { 0.0000000 , 0.0000000 , 0.0000000 },
629  { -2.6737396 , 2.9379143 , -0.3805367 },
630  { -1.4380124 , 2.7756396 , -0.2754625 },
631  { -0.6630196 , 2.6887253 , -0.2216493 }},
632  },
633  /* these are for triplets */
634  {
635  { /* these are P to S */
636  { 0.3075287 , 0.9087130 , -1.0387207 },
637  { 0.687069 , 1.1485864 , -0.6627317 },
638  { 0.9776064 , 1.3382024 , -0.5331906 },
639  { 1.2107725 , 1.4943721 , -0.4779232 }},
640  { /* these are S to P */
641  { -1.3659605 , 2.3262253 , -0.0306439 },
642  { -0.2899490 , 2.3279391 , -0.0298695 },
643  { 0.3678878 , 2.3266603 , -0.0240021 },
644  { 0.8427457 , 2.3249540 , -0.0194091 }},
645  { /* these are D to P */
646  { 1.3108281 , 2.8446367 , -0.1649923 },
647  { 1.8437692 , 2.2399326 , -0.2583398 },
648  { 2.1820792 , 2.0693762 , -0.1864091 },
649  { 2.4414052 , 2.0168255 , -0.1426083 }},
650  { /* these are P to D */
651  { 0.0000000 , 0.0000000 , 0.0000000 },
652  { -1.9219877 , 2.7689624 , -0.2536072 },
653  { -0.7818065 , 2.6595150 , -0.1895313 },
654  { -0.0665624 , 2.5955623 , -0.1522616 }},
655  }
656  };
657 
658  if( lLo==0 )
659  {
660  paramSet = 0;
661  }
662  else if( lLo==1 && lHi==0 )
663  {
664  paramSet = 1;
665  }
666  else if( lLo==1 && lHi==2 )
667  {
668  paramSet = 2;
669  }
670  else if( lLo==2 )
671  {
672  paramSet = 3;
673  ASSERT( lHi==1 );
674  }
675 
676  ASSERT( (int)((sHi-1)/2) == 0 || (int)((sHi-1)/2) == 1 );
677  a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
678  b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
679  c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
680  ASSERT( Enerwn > 0. );
681  x = log( iso_sp[ipHE_LIKE][nelem].fb[ipLo].xIsoLevNIonRyd*RYD_INF/Enerwn );
682 
683  emisOscStr = exp(a+b*x+c*x*x)/pow(Eff_nupper,3.)*
684  (2.*lLo+1)/(2.*lHi+1);
685 
686  Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;
687 
688  if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
689  {
690  Aul *= (2.*(ipLo-3)+1.0)/9.0;
691  }
692 
693  ASSERT( Aul > 0. );
694  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f);
695  }
696  else
697  {
698  /* Calculate the radial integral from the quantum defects. */
699  RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
700  ASSERT( RI2 > 0. );
701  /* Convert radial integral to Aul. */
702  Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
703  /* radial integral routine does not recognize fine structure.
704  * Here we split 2^3P. */
705  if( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
706  {
707  Aul *= (2.*(ipLo-3)+1.0)/9.0;
708  }
709 
710  ASSERT( Aul > 0. );
711  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.07f);
712  }
713  }
714  }
715 
716  /* Heavier species */
717  else
718  {
719  /* Retrieve transition probabilities for Helike ions. */
720  /* >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
721  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J, originally astro.ph. 0201454 */
722  if( ipHi <= MAX_TP_INDEX && N_(ipHi) <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
723  {
724  /*Must not be accessed by collapsed levels! */
725  ASSERT( ipHi < ( iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) );
726  ASSERT( ipLo < ( iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) );
727  ASSERT( ipHi > 2 );
728 
729  Aul = TransProbs[nelem][ipHi][ipLo];
730  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
731  }
732 
733  if( Aul < 0. )
734  {
735  /* Do same-n transitions. */
736  if( nLo==nHi && lHi<=2 && lLo<=2 )
737  {
738  /* These are 2p3Pj to 2s3S fits to (low Z) Porquet & Dubau (2000) &
739  * (high Z) NIST Atomic Spectra Database. */
740  if( ipLo == ipHe2s3S )
741  {
742  if(ipHi == ipHe2p3P0)
743  Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
744  else if(ipHi == ipHe2p3P1)
745  Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
746  else if(ipHi == ipHe2p3P2)
747  Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
748  else
749  {
750  fprintf(ioQQQ," PROBLEM Impossible value for ipHi in he_1trans.\n");
751  TotalInsanity();
752  }
753  }
754 
755  /* These are 2p1P to 2s1S fits to data from TOPbase. */
756  else if( ( ipLo == ipHe2s1S ) && ( ipHi == ipHe2p1P) )
757  {
758  Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
759  }
760 
761  else
762  {
763  /* This case should only be entered if n > 2. Those cases were done above. */
764  ASSERT( nLo > 2 );
765 
766  /* Triplet P to triplet S, delta n = 0 */
767  if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
768  {
769  Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
770  }
771  /* Singlet P to singlet D, delta n = 0 */
772  else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
773  {
774  Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
775  }
776  /* Singlet P to singlet S, delta n = 0 */
777  else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
778  {
779  Aul = 6.646E7 * pow((double)nelem,1.5) / pow((double)nHi, 5.077);
780  }
781  else
782  {
783  ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
784  Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
785  if( (lHi >2) || (lLo > 2) )
786  Aul *= (lHi/2.);
787  if(lLo > 2)
788  Aul *= (1./9.);
789  }
790  }
791  ASSERT( Aul > 0.);
792  }
793 
794  /* assume transitions involving F and higher orbitals are hydrogenic. */
795  else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
796  {
797  Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
798  ASSERT( Aul > 0.);
799  }
800 
801  /* These transitions are of great importance, but the below radial integral
802  * routine fails to achieve desirable accuracy, so these are fits as produced
803  * from He A's for nupper through 9. They are transitions to ground and
804  * 2, 3, and 4 triplet S. */
805  else if( ( ipLo == 0 ) || ( ipLo == ipHe2s1S ) || ( ipLo == ipHe2s3S )
806  || ( ipLo == ipHe3s3S ) || ( nLo==4 && lLo==0 && sLo==3 ) )
807  {
808  /* Here are the Lyman transitions. */
809  if( ipLo == 0 )
810  {
811  ASSERT( (lHi == 1) && (sHi == 1) );
812 
813  /* In theory, this Z dependence should be Z^4, but values from TOPbase
814  * suggest 3.9 is a more accurate exponent. Values from
815  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
816  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
817  /* originally astro.ph. 0201454 */
818  Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
819  }
820 
821  /* Here are the Balmer transitions. */
822  else if( ipLo == ipHe2s1S )
823  {
824  ASSERT( (lHi == 1) && (sHi == 1) );
825 
826  /* More fits, as above. */
827  Aul = 5.0e8 * pow((double)nelem,4.) / pow((double)nHi, 2.889);
828  }
829 
830  /* Here are transitions down to triplet S */
831  else
832  {
833  ASSERT( (lHi == 1) && (sHi == 3) );
834 
835  /* These fits also as above. */
836  if( nLo == 2 )
837  Aul = 1.5 * 3.405E8 * pow((double)nelem,4.) / pow((double)nHi, 2.883);
838  else if( nLo == 3 )
839  Aul = 2.5 * 4.613E7 * pow((double)nelem,4.) / pow((double)nHi, 2.672);
840  else
841  Aul = 3.0 * 1.436E7 * pow((double)nelem,4.) / pow((double)nHi, 2.617);
842  }
843 
844  ASSERT( Aul > 0.);
845  }
846 
847  /* Every other allowed transition is calculated as follows. */
848  else
849  {
850  /* Calculate the radial integral from the quantum defects. */
851  RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem));
852  /* Convert radial integral to Aul. */
853  Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2);
854  /* radial integral routine does not recognize fine structure.
855  * Here we split 2^3P. */
856  if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) ) && (Aul > iso_ctrl.SmallA) )
857  {
858  Aul *= (2.*(ipLo-3)+1.0)/9.0;
859  }
860 
861  }
862  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
863  }
864  }
865  }
866 
867  /* Now do forbidden transitions from 2-1 ... */
868  /* and those given by
869  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
870  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
871  /* originally astro.ph. 0201454
872  * for heavy elements. These are triplet P to singlet S,
873  * going either up or down...Triplet S to Singlet P are not included, as they are far weaker. */
874  else
875  {
876  ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
877  Aul = ForbiddenAuls(ipHi, ipLo, nelem);
878  ASSERT( Aul > 0. );
879  }
880 
881  Aul = MAX2( Aul, iso_ctrl.SmallA );
882  ASSERT( Aul >= iso_ctrl.SmallA );
883 
884  /* negative energy for a transition with substantial transition probability
885  * would be major logical error - but is ok for same-n l transitions */
886  if( Enerwn < 0. && Aul > iso_ctrl.SmallA )
887  {
888  fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
889  }
890 
891  return Aul;
892 }
893 
894 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
895 {
896  /* Every bit of this routine is based upon the singlet-triplet mixing formalism given in
897  * >>refer He FSM Drake, G. W. F. 1996, in Atomic, Molecular, \& Optical Physics Handbook,
898  * >>refercon ed. G. W. F. Drake (New York: AIP Press).
899  * That formalism mixes the levels themselves, but since this code is not J-resolved, we simulate
900  * that by mixing only the transition probabilities. We find results comparable to those calculated
901  * in the fully J-resolved model spearheaded by Rob Bauman, and described in
902  * >>refer He FSM Bauman, R. P., Porter, R. L., Ferland, G. J., \& MacAdam, K. B. 2005, ApJ, accepted */
903  long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
904  double Ass, Att, Ast, Ats;
905  double SinHi, SinLo, CosHi, CosLo;
906  double HiMixingAngle, LoMixingAngle , error;
907  double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
908 
909  DEBUG_ENTRY( "DoFSMixing()" );
910 
911  nHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].n();
912  lHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].l();
913  sHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].S();
914  nLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].n();
915  lLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].l();
916  sLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].S();
917 
918  if( ( sHi == 3 || sLo == 3 ) ||
919  ( abs(lHi - lLo) != 1 ) ||
920  ( nLo < 2 ) ||
921  ( lHi <= 1 || lLo <= 1 ) ||
922  ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
923  ( nHi > nLo && lHi != 1 && lLo != 1 ) )
924  {
925  return;
926  }
927 
928  ASSERT( lHi > 0 );
929  /*ASSERT( (lHi > 1) && (lLo > 1) );*/
930 
931  ipHiTrip = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][3];
932  ipLoTrip = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][3];
933 
934  if( lHi == 2 )
935  {
936  HiMixingAngle = 0.01;
937  }
938  else if( lHi == 3 )
939  {
940  HiMixingAngle = 0.5;
941  }
942  else
943  {
944  HiMixingAngle = PI/4.;
945  }
946 
947  if( lLo == 2 )
948  {
949  LoMixingAngle = 0.01;
950  }
951  else if( lLo == 3 )
952  {
953  LoMixingAngle = 0.5;
954  }
955  else
956  {
957  LoMixingAngle = PI/4.;
958  }
959 
960  /* These would not work correctly if l<=1 were included in this treatment! */
961  ASSERT( ipHiTrip > ipLoTrip );
962  ASSERT( ipHiTrip > ipLoSing );
963  ASSERT( ipHiSing > ipLoTrip );
964  ASSERT( ipHiSing > ipLoSing );
965 
966  SinHi = sin( HiMixingAngle );
967  SinLo = sin( LoMixingAngle );
968  CosHi = cos( HiMixingAngle );
969  CosLo = cos( LoMixingAngle );
970 
971  Kss = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).EnergyWN();
972  Ktt = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).EnergyWN();
973  Kst = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).EnergyWN();
974  Kts = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).EnergyWN();
975 
976  fss = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()/TRANS_PROB_CONST/POW2(Kss)/
977  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Lo()).g()*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Hi()).g();
978 
979  ftt = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul()/TRANS_PROB_CONST/POW2(Ktt)/
980  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Lo()).g()*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Hi()).g();
981 
982  temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
983  fssNew = Kss*POW2( temp );
984  temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
985  fttNew = Ktt*POW2( temp );
986  temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
987  fstNew = Kst*POW2( temp );
988  temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
989  ftsNew = Kts*POW2( temp );
990 
991  Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Lo()).g()/
992  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Hi()).g();
993 
994  Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Lo()).g()/
995  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Hi()).g();
996 
997  Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Lo()).g()/
998  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Hi()).g();
999 
1000  Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Lo()).g()/
1001  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Hi()).g();
1002 
1003  error = fabs( ( iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()+
1004  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul() )/
1005  (Ass+Ast+Ats+Att) - 1.f );
1006 
1007  if( error > 0.001 )
1008  {
1009  fprintf( ioQQQ, "FSM error %e LS %li HS %li LT %li HT %li Ratios Ass %e Att %e Ast %e Ats %e\n", error,
1010  ipLoSing, ipHiSing, ipLoTrip, ipHiTrip,
1011  Ass/iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul(),
1012  Att/iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul(),
1013  Ast/iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Emis().Aul(),
1014  Ats/iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Emis().Aul() );
1015  }
1016 
1017  iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul() = (realnum)Ass;
1018  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul() = (realnum)Att;
1019  iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Emis().Aul() = (realnum)Ast;
1020  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Emis().Aul() = (realnum)Ats;
1021 
1022  return;
1023 }
1024 
1025 /*ritoa converts the square of the radial integral for a transition
1026  * (calculated by scqdri) to the transition probability, Aul. */
1027 STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
1028 {
1029  /* Variables are as follows: */
1030  /* lg = larger of li and lf */
1031  /* fmean = mean oscillator strength */
1032  /* for a given level. */
1033  /* mu = reduced mass of optical electron. */
1034  /* EinsteinA = Einstein emission coef. */
1035  /* w = angular frequency of transition. */
1036  /* RI2_cm = square of rad. int. in cm^2. */
1037  long lg;
1038  double fmean,mu,EinsteinA,w,RI2_cm;
1039 
1040  DEBUG_ENTRY( "ritoa()" );
1041 
1042  mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT));
1043 
1044  w = 2. * PI * k * SPEEDLIGHT;
1045 
1046  RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM;
1047 
1048  lg = (lf > li) ? lf : li;
1049 
1050  fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));
1051 
1052  EinsteinA = TRANS_PROB_CONST*k*k*fmean;
1053 
1054  /* ASSERT( EinsteinA > SMALLFLOAT ); */
1055 
1056  return EinsteinA;
1057 }
1058 
1059 realnum helike_transprob( long nelem, long ipHi, long ipLo )
1060 {
1061  double Aul, Aul1;
1062  double Enerwn, n_eff_hi, n_eff_lo;
1063  long ipISO = ipHE_LIKE;
1064  /* charge to 4th power, needed for scaling laws for As*/
1065  double z4 = POW2((double)nelem);
1066  z4 *= z4;
1067 
1068  DEBUG_ENTRY( "helike_transprob()" );
1069 
1070  Enerwn = iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN();
1071  n_eff_hi = N_(ipHi) - helike_quantum_defect(nelem,ipHi);
1072  n_eff_lo = N_(ipLo) - helike_quantum_defect(nelem,ipLo);
1073 
1074  if( ipHi >= iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max )
1075  {
1076  if( ipLo >= iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max )
1077  {
1078  /* Neither upper nor lower is resolved...Aul()'s are easy. */
1079  Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
1080  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f,0.001f);
1081 
1082  ASSERT( Aul > 0.);
1083  }
1084  else
1085  {
1086  /* Lower level resolved, upper not. First calculate Aul
1087  * from upper level with ang mom one higher. */
1088  Aul = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)+1,
1089  S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
1090 
1091  iso_sp[ipISO][nelem].CachedAs[ N_(ipHi)-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][0] = (realnum)Aul;
1092 
1093  Aul *= (2.*L_(ipLo)+3.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
1094 
1095  if( L_(ipLo) != 0 )
1096  {
1097  /* For all l>0, add in transitions from upper level with ang mom one lower. */
1098  Aul1 = he_1trans( nelem, Enerwn, n_eff_hi, L_(ipLo)-1,
1099  S_(ipLo), -1, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
1100 
1101  iso_sp[ipISO][nelem].CachedAs[ N_(ipHi)-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] = (realnum)Aul1;
1102 
1103  /* now add in this part after multiplying by stat weight for lHi = lLo-1. */
1104  Aul += Aul1*(2.*L_(ipLo)-1.) * S_(ipLo) / (4.*(double)N_(ipHi)*(double)N_(ipHi));
1105  }
1106  else
1107  iso_sp[ipISO][nelem].CachedAs[ N_(ipHi)-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] = 0.f;
1108 
1109  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.01f);
1110  ASSERT( Aul > 0.);
1111  }
1112  }
1113  else
1114  {
1115  /* Both levels are resolved...do the normal bit. */
1116  if( Enerwn < 0. )
1117  {
1118  Aul = he_1trans( nelem, -1.*Enerwn, n_eff_lo,
1119  L_(ipLo), S_(ipLo), ipLo-3, n_eff_hi, L_(ipHi), S_(ipHi), ipHi-3);
1120  }
1121  else
1122  {
1123  Aul = he_1trans( nelem, Enerwn, n_eff_hi,
1124  L_(ipHi), S_(ipHi), ipHi-3, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3);
1125  }
1126  }
1127 
1128  return (realnum)Aul;
1129 }
1130 
1131 /* This routine is pure infrastructure and bookkeeping, no physics. */
1132 void HelikeTransProbSetup( void )
1133 {
1134 
1135  const int chLine_LENGTH = 1000;
1136  char chLine[chLine_LENGTH];
1137 
1138  FILE *ioDATA;
1139  bool lgEOL;
1140 
1141  long nelem, ipLo, ipHi, i, i1, i2, i3;
1142 
1143  DEBUG_ENTRY( "HelikeTransProbSetup()" );
1144 
1145  TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
1146 
1147  for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
1148  {
1149 
1150  TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) );
1151 
1152  for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo )
1153  {
1154  TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX );
1155  }
1156  }
1157 
1158  /********************************************************************/
1159  /*************** Read in data from he_transprob.dat *****************/
1160  if( trace.lgTrace )
1161  fprintf( ioQQQ," HelikeTransProbSetup opening he_transprob.dat:");
1162 
1163  ioDATA = open_data( "he_transprob.dat", "r" );
1164 
1165  /* check that magic number is ok */
1166  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1167  {
1168  fprintf( ioQQQ, " HelikeTransProbSetup could not read first line of he_transprob.dat.\n");
1169  cdEXIT(EXIT_FAILURE);
1170  }
1171  i = 1;
1172  i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1173  i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1174  if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
1175  {
1176  fprintf( ioQQQ,
1177  " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1178  fprintf( ioQQQ,
1179  " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1180  TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
1181  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1182  cdEXIT(EXIT_FAILURE);
1183  }
1184 
1185  /* Initialize TransProbs[nelem][ipLo][ipHi] before filling it in. */
1186  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
1187  {
1188  for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
1189  {
1190  for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
1191  {
1192  TransProbs[nelem][ipHi][ipLo] = -1.;
1193  }
1194  }
1195  }
1196 
1197  for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
1198  {
1199  char *chTemp;
1200 
1201  /* get next line image */
1202  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1203  BadRead();
1204 
1205  while( chLine[0]=='#' )
1206  {
1207  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1208  BadRead();
1209  }
1210 
1211  i3 = 1;
1212  i1 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
1213  i2 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
1214  /* check that these numbers are correct */
1215  if( i1<0 || i2<=i1 )
1216  {
1217  fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1218  cdEXIT(EXIT_FAILURE);
1219  }
1220 
1221  chTemp = chLine;
1222 
1223  /* skip over 2 tabs to start of data */
1224  for( i=0; i<1; ++i )
1225  {
1226  if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
1227  {
1228  fprintf( ioQQQ, " HelikeTransProbSetup could not init he_transprob\n" );
1229  cdEXIT(EXIT_FAILURE);
1230  }
1231  ++chTemp;
1232  }
1233 
1234  /* now read in the data */
1235  for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
1236  {
1237  if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
1238  {
1239  fprintf( ioQQQ, " HelikeTransProbSetup could not scan he_transprob\n" );
1240  cdEXIT(EXIT_FAILURE);
1241  }
1242  ++chTemp;
1243 
1244  sscanf( chTemp , "%le" , &TransProbs[nelem][i2][i1] );
1245 
1247  if( lgEOL )
1248  {
1249  fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1250  cdEXIT(EXIT_FAILURE);
1251  }
1252  }
1253  }
1254 
1255  /* check that ending magic number is ok */
1256  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1257  {
1258  fprintf( ioQQQ, " HelikeTransProbSetup could not read last line of he_transprob.dat.\n");
1259  cdEXIT(EXIT_FAILURE);
1260  }
1261  i = 1;
1262  i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1263  i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1264  if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
1265  {
1266  fprintf( ioQQQ,
1267  " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1268  fprintf( ioQQQ,
1269  " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1270  TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
1271  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1272  cdEXIT(EXIT_FAILURE);
1273  }
1274 
1275  /* close the data file */
1276  fclose( ioDATA );
1277  return;
1278 }
chLine_LENGTH
#define chLine_LENGTH
helike.h
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
iso_put_error
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
H_Einstein_A
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
Definition: hydro_bauman.cpp:2324
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
dense
t_dense dense
Definition: dense.cpp:24
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
BadRead
NORETURN void BadRead(void)
Definition: service.cpp:901
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
TransProbs
static double *** TransProbs
Definition: helike_einsta.cpp:25
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
realnum
float realnum
Definition: cddefines.h:103
helike_quantum_defect
double helike_quantum_defect(long int nelem, long int ipLev)
Definition: helike_energy.cpp:243
ipHe3s3S
const int ipHe3s3S
Definition: iso.h:52
STATIC
#define STATIC
Definition: cddefines.h:97
vJint
static double vJint
Definition: helike_einsta.cpp:40
Jint
STATIC double Jint(double theta)
Definition: helike_einsta.cpp:43
hydro_bauman.h
ipHe3s1S
const int ipHe3s1S
Definition: iso.h:53
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
L_
#define L_(A_)
Definition: iso.h:21
trace.h
ipNEON
const int ipNEON
Definition: cddefines.h:314
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
scqdri
double scqdri(double nstar, long int l, double npstar, long int lp, double iz)
Definition: helike_einsta.cpp:146
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
AngerJ
STATIC double AngerJ(double vv, double zz)
Definition: helike_einsta.cpp:58
POW2
#define POW2
Definition: cddefines.h:929
MIN2
#define MIN2
Definition: cddefines.h:761
ATOMIC_MASS_UNIT
const UNUSED double ATOMIC_MASS_UNIT
Definition: physconst.h:88
PI
const UNUSED double PI
Definition: physconst.h:29
helike_einsta.h
MAX_TP_INDEX
#define MAX_TP_INDEX
Definition: helike_einsta.h:9
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
dense.h
eina
double eina(double gf, double enercm, double gup)
Definition: lines_service.cpp:84
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
ipHe3p1P
const int ipHe3p1P
Definition: iso.h:57
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
H_BAR
const UNUSED double H_BAR
Definition: physconst.h:144
IPRAD
#define IPRAD
Definition: iso.h:86
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
SPEEDLIGHT
const UNUSED double SPEEDLIGHT
Definition: physconst.h:100
ELECTRON_MASS
const UNUSED double ELECTRON_MASS
Definition: physconst.h:91
MAX2
#define MAX2
Definition: cddefines.h:782
LIMELM
const int LIMELM
Definition: cddefines.h:258
RYD_INF
const UNUSED double RYD_INF
Definition: physconst.h:115
N_HE1_TRANS_PROB
#define N_HE1_TRANS_PROB
Definition: helike_einsta.h:7
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
HydroEinstA
double HydroEinstA(long int n1, long int n2)
Definition: hydroeinsta.cpp:12
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:49
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:567
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:894
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
ipHe3p3P
const int ipHe3p3P
Definition: iso.h:54
physconst.h
zJint
static double zJint
Definition: helike_einsta.cpp:40
he_1trans
double he_1trans(long nelem, double Enerwn, double Eff_nupper, long lHi, long sHi, long jHi, double Eff_nlower, long lLo, long sLo, long jLo)
Definition: helike_einsta.cpp:501
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
fixit
void fixit(void)
Definition: service.cpp:991
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
TRANSPROBMAGIC
#define TRANSPROBMAGIC
Definition: helike.h:22
TRANS_PROB_CONST
const UNUSED double TRANS_PROB_CONST
Definition: physconst.h:237
taulines.h
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1059
ritoa
STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
Definition: helike_einsta.cpp:1027
ipHe2s1S
const int ipHe2s1S
Definition: iso.h:45
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1132
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
qg32
double qg32(double, double, double(*)(double))
Definition: service.cpp:1053
S_
#define S_(A_)
Definition: iso.h:22
ipARGON
const int ipARGON
Definition: cddefines.h:322
ForbiddenAuls
STATIC double ForbiddenAuls(long ipHi, long ipLo, long nelem)
Definition: helike_einsta.cpp:221
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
BOHR_RADIUS_CM
const UNUSED double BOHR_RADIUS_CM
Definition: physconst.h:222
hydroeinsta.h
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
N_
#define N_(A_)
Definition: iso.h:20
strchr_s
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1439
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17