cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "dense.h"
15 #include "trace.h"
16 #include "hydro_bauman.h"
17 #include "iso.h"
18 #include "helike.h"
19 #include "helike_einsta.h"
20 #include "freebound.h"
21 #include "lines_service.h"
22 #include "integrate.h"
23 
24 /* the array of transitions probabilities read from data file. */
25 static double ***TransProbs;
26 
27 STATIC double helike_transprob_collapsed_to_collapsed( long nelem, long nHi, long nLo, double Enerwn );
28 STATIC double helike_transprob_collapsed_to_resolved( long nelem, long nHi, long nLo, long lLo, long sLo, long jLo, double Enerwn );
29 
30 /*ritoa converts the square of the radial integral for a transition
31  * (calculated by scqdri) to the transition probability, Aul. */
32 STATIC double ritoa( long li, long lf, long nelem, double k, double RI2 );
33 
34 /* handles all forbidden transition probabilities. */
35 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem );
36 
37 /*
38 static long RI_ipHi, RI_ipLo, RI_ipLev;
39 static long globalZ;
40 */
41 
42 /* used as parameters in qg32 integration */
43 static double vJint , zJint;
44 
45 /* the integrand used in the AngerJ function below. */
46 STATIC double Jint( double theta )
47 {
48  /* [ cos[vx - z sin[x]] ] */
49  double
50  d1 = vJint * theta,
51  d2 = zJint * sin(theta),
52  d3 = (d1 - d2),
53  d4 = cos(d3);
54 
55  return( d4 );
56 }
57 
58 /* AngerJ function. */
59 // See also http://dlmf.nist.gov/11.10 for series expansion and recurrences
60 STATIC double AngerJ_mac( double vv, double zz )
61 {
62  // Maclaurin series for J_nu(z)
63  // http://dlmf.nist.gov/11.10#iii
64  // Series is absolutely convergent, but can take
65  // numerous terms for large z
66  // Problems when vv=2 -- should fall back to Bessel functions?
67  // Worth trying the acceleration techinques of http://dlmf.nist.gov/3.9
68  double g1=tgamma(1+0.5*vv), g2=tgamma(1-0.5*vv),
69  g3=tgamma(1.5+0.5*vv), g4=tgamma(1.5-0.5*vv);
70  double y1 = 0., y2=0., zp = 1.;
71  for (long k=0; k<100; ++k)
72  {
73  double y10 = y1;
74  y1 += zp/(g1*g2);
75  if (y1 == y10)
76  break;
77  zp *= 0.5*zz;
78  // Gamma(t+1) = t*Gamma(t)
79  g1 *= k+0.5*vv+1;
80  g2 *= k-0.5*vv+1;
81  y2 += zp/(g3*g4);
82  zp *= -0.5*zz;
83  g3 *= k+0.5*vv+1.5;
84  g4 *= k-0.5*vv+1.5;
85  }
86  double cosfac=cos(0.5*PI*vv), sinfac=sin(0.5*PI*vv);
87  return y1*cosfac+y2*sinfac;
88 }
89 
90 STATIC double bessjnu( double vv, double zz )
91 {
92  // DLMF 10.17.3
93  // Asymptotic series for Bessel function for non-integer index
94  double rz=1./zz, akrz=1., bff=0., bgg=0., bfp=0.;
95  for (long k=0; k<100; ++k)
96  {
97  if (k > 0 && fabs(akrz) > fabs(bfp))
98  break;
99  bfp = akrz;
100  bff += akrz;
101  akrz *= (4*vv*vv-(4*k+1)*(4*k+1))/(8*(2*k+1)*zz);
102  bgg += akrz;
103  akrz *= -(4*vv*vv-(4*k+3)*(4*k+3))/(8*(2*k+2)*zz);
104  }
105  double omega = zz-(2*vv+1)*PI/4;
106  return sqrt(2/PI*rz)*(cos(omega)*bff-sin(omega)*bgg);
107 }
108 STATIC double AngerJ_asymp( double vv, double zz )
109 {
110  double ff = 0., gg=0., fk=1., gk=1., rz=1./zz, rzsq = rz*rz,
111  fp=0.;
112  for (long k=0; k<100; ++k)
113  {
114  if (k > 0)
115  {
116  fk *= (vv*vv-(2*k-1)*(2*k-1))*rzsq;
117  gk *= (vv*vv-4*k*k)*rzsq;
118  if (fk > fp)
119  break;
120  }
121  fp = fk;
122  ff += fk;
123  gg += gk;
124  }
125  return bessjnu(vv,zz)+sin(PI*vv)/(PI*zz)*(ff-(vv*rz)*gg);
126 }
127 STATIC double AngerJ( double vv, double zz )
128 {
129  long int rep = 0, ddiv, divsor;
130 
131  double y = 0.0;
132 
133  DEBUG_ENTRY( "AngerJ()" );
134 
135  /* Estimate number of peaks in integrand. */
136  /* Divide region of integration by number */
137  /* peaks in region. */
138  if( (fabs(vv)) - (int)(fabs(vv)) > 0.5 )
139  ddiv = (int)(fabs(vv)) + 1;
140  else
141  ddiv = (int)(fabs(vv));
142 
143  divsor = ((ddiv == 0) ? 1 : ddiv);
144  vJint = vv;
145  zJint = zz;
146 
147  for( rep = 0; rep < divsor; rep++ )
148  {
149  double
150  rl = (((double) rep)/((double) divsor)),
151  ru = (((double) (rep+1))/((double) divsor)),
152  x_low = (PI * rl),
153  x_up = (PI * ru);
154 
155  y += qg32( x_low, x_up, Jint );
156  }
157 
158  return( y / PI );
159 }
160 
161 /******************************************************************************/
162 /******************************************************************************/
163 /* */
164 /* Semi-Classical Quantum Defect Radial Integral */
165 /* */
166 /* See for example */
167 /* Atomic, Molecular & Optical Physics Handbook */
168 /* Gordon W. F. Drake; Editor */
169 /* AIP Press */
170 /* Woddbury, New York. */
171 /* 1996 */
172 /* */
173 /* NOTE:: we do not include the Bohr Radius a_o in the */
174 /* definition of of R(n,L;n'L') as per Drake. */
175 /* */
176 /* */
177 /* 1 (n_c)^2 | { D_l max(L,L') } */
178 /* R(n,L;n'L') = --- ------- | { 1 - ------------- } J( D_n-1; -x ) - */
179 /* Z 2 D_n | { n_c } */
180 /* */
181 /* */
182 /* { D_L max(L,L') } */
183 /* - { 1 + ------------- } J( D_n+1; -x ) */
184 /* { n_c } */
185 /* */
186 /* */
187 /* 2 | */
188 /* + --- sin(Pi D_n)(1-e) | */
189 /* Pi | */
190 /* */
191 /* where */
192 /* n_c = (2n*'n*)/(n*'+n*) */
193 /* */
194 /* Here is the quantity Drake gives... */
195 /* n_c = ( 2.0 * nstar * npstar ) / ( nstar + npstar ); */
196 /* */
197 /* while V.A. Davidkin uses */
198 /* n_c = sqrt( nstar * npstar ); */
199 /* */
200 /* D_n = n*' - n* */
201 /* */
202 /* D_L = L*' - L* */
203 /* */
204 /* x = e D_n */
205 /* */
206 /* Lmx = max(L',L) */
207 /* */
208 /* e = sqrt( 1 - {Lmx/n_c}^2 ) */
209 /* */
210 /* */
211 /* Here n* = n - qd where qd is the quantum defect */
212 /* */
213 /******************************************************************************/
214 /******************************************************************************/
215 STATIC double scqdri(/* upper and lower quantum numbers...n's are effective */
216  double nstar, long int l,
217  double npstar, long int lp,
218  double iz
219  )
220 {
221  double n_c = ((2.0 * nstar * npstar ) / ( nstar + npstar ));
222 
223  double D_n = (nstar - npstar);
224  double D_l = (double) ( l - lp );
225  double lg = (double) ( (lp > l) ? lp : l);
226 
227  double h = (lg/n_c);
228  double g = h*h;
229  double f = ( 1.0 - g );
230  double e = (( f >= 0.0) ? sqrt( f ) : 0.0 );
231 
232  double x = (e * D_n);
233  double z = (-1.0 * x);
234  double v1 = (D_n + 1.0);
235  double v2 = (D_n - 1.0);
236 
237  double d1,d2,d7,d8,d9,d34,d56,d6_1;
238 
239  DEBUG_ENTRY( "scqdri()" );
240 
241  if( iz == 0.0 )
242  iz += 1.0;
243 
244  if( D_n == 0.0 )
245  {
246  return( -1.0 );
247  }
248 
249  if( D_n < 0.0 )
250  {
251  return( -1.0 );
252  }
253 
254  if( f < 0.0 )
255  {
256  /* This can happen for certain quantum defects */
257  /* in the lower n=1:l=0 state. In which case you */
258  /* probably should be using some other alogrithm */
259  /* or theory to calculate the dipole moment. */
260  return( -1.0 );
261  }
262 
263  d1 = ( 1.0 / iz );
264 
265  d2 = (n_c * n_c)/(2.0 * D_n);
266 
267  if (0)
268  {
269  // Confirm accuracy of result of AngerJ function
270 
271  // relative accuracy is ~1e-10 in the worst case, where this error
272  // may well be in summing series representation
273  double vv = v2, zz = z, y = AngerJ(vv,zz);
274  double yt = AngerJ_mac(vv,zz);
275  double vp = (zz > 0.) ? vv : -vv;
276  double ya = fabs(zz) < fabs(vv) ? 0. : AngerJ_asymp(vp,fabs(zz));
277  printf("ANGER %15.8g %15.8g: %15.8g %15.8g (%15.8g) asymp %15.8g\n",
278  zz,vv,yt,y,fabs(yt-y)/(1e-37+fabs(yt)),ya);
279  }
280 
281  d34 = (1.0 - ((D_l * lg)/n_c)) * AngerJ( v1, z );
282 
283  d56 = (1.0 + ((D_l * lg)/n_c)) * AngerJ( v2, z );
284 
285  d6_1 = PI * D_n;
286 
287  d7 = (2./PI) * sin( d6_1 ) * (1.0 - e);
288 
289  d8 = d1 * d2 * ( (d34) - (d56) + d7 );
290 
291  d9 = d8 * d8;
292 
293  ASSERT( D_n > 0.0 );
294  ASSERT( l >= 0 );
295  ASSERT( lp >= 0 );
296  ASSERT( (l == lp + 1) || ( l == lp - 1) );
297  ASSERT( n_c != 0.0 );
298  ASSERT( f >= 0.0 );
299  ASSERT( d9 > 0.0 );
300 
301  return( d9 );
302 }
303 
304 STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem )
305 {
306  double A;
307 
308  DEBUG_ENTRY( "ForbiddenAuls()" );
309 
310  int ipISO = ipHE_LIKE;
311 
312  if( (ipLo == ipHe1s1S) && (N_(ipHi) == 2) )
313  {
314  if( nelem == ipHELIUM )
315  {
316  // All of these but the second and third one (with values 1E-20) are from
317  // >>refer HeI As Lach, G, & Pachucki, K, 2001, Phys. Rev. A 64, 042510
318  // 1E-20 is made up
319  double ForbiddenHe[5] = { 1.272E-4, 1E-20, 1E-20, 177.58, 0.327 };
320 
321  A = ForbiddenHe[ipHi - 1];
322  iso_put_error(ipHE_LIKE,nelem,ipHe2s3S ,ipHe1s1S,IPRAD, 0.01f, 0.20f);
323  iso_put_error(ipHE_LIKE,nelem,ipHe2s1S ,ipHe1s1S,IPRAD, 0.01f, 0.05f);
324  iso_put_error(ipHE_LIKE,nelem,ipHe2p3P0,ipHe1s1S,IPRAD, 0.00f, 0.00f);
325  iso_put_error(ipHE_LIKE,nelem,ipHe2p3P1,ipHe1s1S,IPRAD, 0.01f, 0.05f);
326  iso_put_error(ipHE_LIKE,nelem,ipHe2p3P2,ipHe1s1S,IPRAD, 0.01f, 0.01f);
327  }
328  else
329  {
330  switch ( (int)ipHi )
331  {
332  case 1: /* Parameters for 2^3S to ground transition. */
333  /* >>refer Helike As Lin, C.D., Johnson, W.R., and Dalgarno, A. 1977,
334  * >>refercon Phys. Rev. A 15, 1, 010015 */
335  // 2nu is treated elsewhere
336  A = (3.9061E-7) * pow( (double)nelem+1., 10.419 );
337  break;
338  case 2: /* Parameters for 2^1S to ground transition. */
339  A = iso_ctrl.SmallA; // 2nu is treated elsewhere
340  break;
341  case 3: /* Parameters for 2^3P0 to ground transition. */
342  A = iso_ctrl.SmallA;
343  break;
344  case 4: /* Parameters for 2^3P1 to ground transition. */
345  /* >>chng 06 jul 25, only use the fit to Johnson et al. values up to and
346  * including argon, where there calculation stops. For higher Z use below */
347  if( nelem <= ipARGON )
348  {
349  A = ( 11.431 * pow((double)nelem, 9.091) );
350  }
351  else
352  {
353  /* a fit to the Lin et al. 1977. values, which go past zinc. */
354  A = ( 383.42 * pow((double)nelem, 7.8901) );
355  }
356  break;
357  case 5: /* Parameters for 2^3P2 to ground transition. */
358  /* fit to Lin et al. 1977 values. This fit is good
359  * to 7% for the range from carbon to iron. The Lin et al. values
360  * differ from the Hata and Grant (1984) values (only given from
361  * calcium to copper) by less than 2%. */
362  A = ( 0.11012 * pow((double)nelem, 7.6954) );
363  break;
364  default:
365  TotalInsanity();
366  }
367  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
368  }
369 
370  return A;
371  }
372 
373  /* The next two cases are fits to probabilities given in
374  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
375  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
376  /* originally astro.ph. 0201454 */
377  /* The involve Triplet P and Singlet S. Rates for Triplet S to Singlet P
378  * do not seem to be available. */
379 
380  /* Triplet P to Singlet S...Delta n not equal to zero! */
381  else if( nelem>ipHELIUM && L_(ipHi)==1 && S_(ipHi)==3 &&
382  L_(ipLo)==0 && S_(ipLo)==1 && N_(ipLo) < N_(ipHi) )
383  {
384  A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) /
385  pow((double)N_(ipHi),2.877);
386  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
387  }
388 
389  /* Singlet S to Triplet P...Delta n not equal to zero! */
390  else if( nelem > ipHELIUM && L_(ipHi)==0 && S_(ipHi)==1 &&
391  L_(ipLo)==1 && S_(ipLo)==3 && N_(ipLo) < N_(ipHi) )
392  {
393  A = 2.416 * exp(-0.256*N_(ipLo)) * pow((double)nelem,9.159) / pow((double)N_(ipHi),3.111);
394 
395  if( ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P2) ) )
396  {
397  /* This is divided by 3 instead of 9, because value calculated is specifically for 2^3P1.
398  * Here we assume statistical population of the other two. */
399  A *= (2.*(ipLo-3)+1.0)/3.0;
400  }
401  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
402  }
403 
404  else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe3s3S ) )
405  {
406  double As_3TripS_to_2TripS[9] = {
407  7.86E-9, 4.59E-6, 1.90E-4, 2.76E-3, 2.25E-2,
408  1.27E-1, 5.56E-1, 2.01E0, 6.26E0 };
409 
410  /* These M1 transitions given by
411  * >>refer He-like As Savukov, I.M., Labzowsky, and Johnson, W.R. 2005, PhysRevA, 72, 012504 */
412  if( nelem <= ipNEON )
413  {
414  A = As_3TripS_to_2TripS[nelem-1];
415  /* 20% error is based on difference between Savukov, Labzowsky, and Johnson (2005)
416  * and Lach and Pachucki (2001) for the helium transition. */
417  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f,0.2f);
418  }
419  else
420  {
421  /* This is an extrapolation to the data given above. The fit reproduces
422  * the above values to 10% or better. */
423  A= 7.22E-9*pow((double)nelem, 9.33);
424  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f,0.3f);
425  }
426  }
427 
428  else if( ( ipLo == ipHe2s3S ) && ( ipHi == ipHe2p1P ) )
429  {
430  /* This transition,1.549 , given by Lach and Pachucki, 2001 for the atom */
431  if( nelem == ipHELIUM )
432  {
433  A = 1.549;
434  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
435  }
436  else
437  {
438  /* This is a fit to data given in
439  * >>refer He-like As Savukov, I.M., Johnson, W.R., & Safronova, U.I.
440  * >>refercon astro-ph 0205163 */
441  A= 0.1834*pow((double)nelem, 6.5735);
442  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
443  }
444  }
445 
446  else if( nelem==ipHELIUM && ipHi==4 && ipLo==3 )
447  {
448  /* >>refer He As Bulatov, N.N. 1976, Soviet Astronomy, 20, 315 */
449  fixit("find a transition probability for this 2^3P0 - 2^3P1 transition.");
450  /* This is the 29.6 GHz line that can be seen in radio galaxies. */
455  A = 5.78E-12;
456  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
457 
458  }
459 
460  else if( nelem==ipHELIUM && ipHi==5 && ipLo==4 )
461  {
462  fixit("find a transition probability for this 2^3P1 - 2^3P2 transition.");
463  /* This is the 3 GHz line that can be seen in radio galaxies. */
468  A = 3.61E-15;
469  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
470  }
471 
472  else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe1s1S )
473  {
474  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
475  A = 44.326 * (1./3.) + 0.1146547 * (5./9.);
476  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
477  }
478 
479  else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe2s1S )
480  {
481  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
482  A = 1.459495 * (1./3.) + 3.6558e-5 * (5./9.);
483  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
484  }
485 
486  else if( nelem==ipHELIUM && ipHi==ipHe3p3P && ipLo==ipHe3s1S )
487  {
488  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
489  A = 2.2297e-3 * (1./3.);
490  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
491  }
492 
493  else if( nelem==ipHELIUM && ipHi==ipHe3p1P && ipLo==ipHe2s3S )
494  {
495  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
496  A = 0.1815436;
497  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
498  }
499 
500  else if( nelem==ipHELIUM && ipHi==ipHe3p1P && ipLo==ipHe3s3S )
501  {
502  // >>refer He As Morton, D. \& Drake, G.W.F. 2011, PhysRevA, 83, 042503
503  A = 0.14895852;
504  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
505  }
506 
507  else if( ( ipLo==0 && L_(ipHi)==2 && S_(ipHi)==1 ) ||
508  ( ipLo==ipHe2s3S && L_(ipHi)==2 && S_(ipHi)==3 ) ||
509  ( ipLo==ipHe2s1S && L_(ipHi)==2 && S_(ipHi)==1 ) ||
510  ( N_(ipLo)==2 && L_(ipLo)==1 && S_(ipLo)==3 && N_(ipHi)>=3 && L_(ipHi)==1 && S_(ipHi)==3 ) ||
511  ( ipLo==ipHe2p1P && L_(ipHi)==1 && S_(ipHi)==1 ) )
512  {
513  // >>refer Helike QOS Cann, N. M. \& Thakkar, A. J. 2002, J.Phys.B, 35, 421
514 
515  static const double f_params[5][4][3] = {
516  {
517  /* 1,0,3,2,1 */ {9.360591E-007, -3.1937E-006, 3.5186E-006},
518  /* 1,0,4,2,1 */ {4.631435E-007, -1.4973E-006, 1.4848E-006},
519  /* 1,0,5,2,1 */ {2.493912E-007, -7.8658E-007, 7.3994E-007},
520  /* 1,0,6,2,1 */ {1.476742E-007, -4.5953E-007, 4.1932E-007}},
521  {
522  /* 2,0,3,2,3 */ {1.646733E-006, -2.0028E-006, -1.3552E-006},
523  /* 2,0,4,2,3 */ {9.120593E-008, 3.1301E-007, -3.2190E-007},
524  /* 2,0,5,2,3 */ {1.360965E-008, 1.1467E-007, 8.6977E-008},
525  /* 2,0,6,2,3 */ {3.199421E-009, 4.5485E-008, 1.1016E-007}},
526  {
527  /* 2,0,3,2,1 */ {1.646733E-006, -2.9720E-006, 1.5367E-006},
528  /* 2,0,4,2,1 */ {9.120593E-008, -3.9037E-008, 3.9156E-008},
529  /* 2,0,5,2,1 */ {1.360965E-008, 1.4671E-008, 1.5626E-008},
530  /* 2,0,6,2,1 */ {3.199421E-009, 8.9806E-009, 1.2436E-008}},
531  {
532  /* 2,1,3,1,3 */ {1.543812E-007, -2.8220E-007, 3.0261E-008},
533  /* 2,1,4,1,3 */ {3.648237E-008, -6.6824E-008, 4.5879E-009},
534  /* 2,1,5,1,3 */ {1.488556E-008, -2.7304E-008, 1.6628E-009},
535  /* 2,1,6,1,3 */ {7.678610E-009, -1.4112E-008, 6.8453E-010}},
536  {
537  /* 2,1,3,1,1 */ {1.543812E-007, -2.8546E-007, 4.6605E-008},
538  /* 2,1,4,1,1 */ {3.648237E-008, -6.8422E-008, 1.7038E-008},
539  /* 2,1,5,1,1 */ {1.488556E-008, -2.8170E-008, 8.5012E-009},
540  /* 2,1,6,1,1 */ {7.678610E-009, -1.4578E-008, 4.6686E-009}}
541  };
542 
543  ASSERT( ipLo <= 6 );
544  long index_lo;
545  if( ipLo==ipHe2p1P )
546  index_lo = 4;
547  else if( ipLo==ipHe2p3P1 || ipLo==ipHe2p3P2 )
548  index_lo = 3;
549  else
550  index_lo = ipLo;
551 
552  ASSERT( N_(ipHi) >= 3 );
553  long index_hi = MIN2( N_(ipHi), 6 ) - 3;
554  double f_lu = POW2(nelem+1) * (
555  f_params[index_lo][index_hi][0] +
556  f_params[index_lo][index_hi][1]/(nelem+1) +
557  f_params[index_lo][index_hi][2]/POW2(nelem+1) );
558 
559  // extrapolate for higher upper n
560  if( N_(ipHi) > 6 )
561  f_lu *= pow3( 6. / N_(ipHi) );
562 
563  double gLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].g();
564  double gHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].g();
565  double eWN = iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN();
566 
567  A = eina( gLo * f_lu, eWN, gHi );
568  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
569  }
570 
571  else
572  {
573  /* Current transition is not supported. */
574  A = iso_ctrl.SmallA;
575  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f);
576  }
577 
578  ASSERT( A > 0.);
579 
580  return A;
581 }
582 
583 /* Calculates Einstein A for a given transition. */
584 double he_1trans(
585  /* charge on c scale, Energy is wavenumbers, Einstein A */
586  long nelem , double Enerwn ,
587  /* quantum numbers of upper level: */
588  double Eff_nupper, long lHi, long sHi, long jHi,
589  /* and of lower level: */
590  double Eff_nlower, long lLo, long sLo, long jLo,
591  realnum *error1, realnum *error2 )
592  /* Note j is only necessary for 2 triplet P...for all other n,l,s,
593  * j is completely ignored. */
594 {
595  double RI2, Aul;
596  long nHi, nLo, ipHi, ipLo;
597 
598  DEBUG_ENTRY( "he_1trans()" );
599 
600  ASSERT(nelem > ipHYDROGEN);
601 
602  /* Since 0.4 is bigger than any defect, adding that to the effective principle quantum number,
603  * and truncating to an integer will produce the principal quantum number. */
604  nHi = (int)(Eff_nupper + 0.4);
605  nLo = (int)(Eff_nlower + 0.4);
606 
607  /* Make sure this worked correctly. */
608  ASSERT( fabs(Eff_nupper-(double)nHi) < 0.4 );
609  ASSERT( fabs(Eff_nlower-(double)nLo) < 0.4 );
610 
611  ipHi = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][sHi];
612  if( (nHi==2) && (lHi==1) && (sHi==3) )
613  {
614  ASSERT( (jHi>=0) && (jHi<=2) );
615  ipHi -= (2 - jHi);
616  }
617 
618  ipLo = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][sLo];
619  if( (nLo==2) && (lLo==1) && (sLo==3) )
620  {
621  ASSERT( (jLo>=0) && (jLo<=2) );
622  ipLo -= (2 - jLo);
623  }
624 
625  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == nHi );
626  if( nHi <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
627  {
628  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].l() == lHi );
629  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipHi].S() == sHi );
630  }
631  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].n() == nLo );
632  if( nLo <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
633  {
634  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].l() == lLo );
635  ASSERT( iso_sp[ipHE_LIKE][nelem].st[ipLo].S() == sLo );
636  }
637 
638  /* First do allowed transitions */
639  if( (sHi == sLo) && (abs((int)(lHi - lLo)) == 1) )
640  {
641  Aul = -2.;
642 
643  /* For clarity, let's do this in two separate chunks...one for helium, one for everything else. */
644  if( nelem == ipHELIUM )
645  {
646  /* Retrieve transition probabilities for Helium. */
647  /* >>refer He As Drake, G.W.F., Atomic, Molecular, and Optical Physics Handbook */
648  if( ipHi <= MAX_TP_INDEX && nHi <= iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max )
649  {
650  /*Must not be accessed by collapsed levels! */
651  ASSERT( ipHi < ( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max ) );
652  ASSERT( ipLo < ( iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max ) );
653  ASSERT( ipHi > 2 );
654 
655  Aul = TransProbs[nelem][ipHi][ipLo];
656 
657  *error1 = 0.0001f;
658  *error2 = 0.002f;
659  }
660  else
661  {
662  ASSERT( nHi > 2 );
663  long ipHiTemp = iso_sp[ipHE_LIKE][nelem].IndexIfAllResolved[nHi][lHi][sHi];
664  long ipLoTemp = iso_sp[ipHE_LIKE][nelem].IndexIfAllResolved[nLo][lLo][sLo];
665  if( (nLo==2) && (lLo==1) && (sLo==3) )
666  {
667  ASSERT( (jLo>=0) && (jLo<=2) );
668  ipLoTemp -= (2 - jLo);
669  }
670  ASSERT( ipLoTemp < ipHiTemp );
671 
672  if( ipHiTemp <= MAX_TP_INDEX )
673  {
674  Aul = TransProbs[nelem][ipHiTemp][ipLoTemp];
675  *error1 = 0.0001f;
676  *error2 = 0.002f;
677  }
678  }
679 
680  if( Aul < 0. )
681  {
682  /* Here are the Lyman transitions. */
683  if( nLo==1 )
684  {
685  ASSERT( (lHi == 1) && (sHi == 1) );
686 
687  /* these fits calculated from Drake A's (1996) */
688  Aul = (1.59208e10) / pow3(Eff_nupper);
689  ASSERT( Aul > 0.);
690  *error1 = 0.0002f;
691  *error2 = 0.002f;
692  }
693 
694  /* last resort for transitions involving significant defects,
695  * except that highest lLo are excluded */
696  else if( lHi>=2 && lLo>=2 && nHi>nLo )
697  {
698  Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
699  ASSERT( Aul > 0.);
700  *error1 = 0.006f;
701  *error2 = 0.04f;
702  }
703  /* These fits come from extrapolations of Drake's oscillator strengths
704  * out to the series limit. We also use this method to obtain threshold
705  * photoionization cross-sections for the lower level of each transition here. */
706  /* See these two references :
707  * >>refer He As Hummer, D. G. \& Storey, P. J. 1998, MNRAS 297, 1073
708  * >>refer Seaton's Theorem Seaton, M. J. 1958, MNRAS 118, 504 */
709  else if( nHi>10 && nLo<=5 && lHi<=2 && lLo<=2 )
710  {
711  int paramSet=0;
712  double emisOscStr, x, a, b, c;
713  static const double extrapol_Params[2][4][4][3] = {
714  /* these are for singlets */
715  {
716  { /* these are P to S */
717  { 0.8267396 , 1.4837624 , -0.4615955 },
718  { 1.2738405 , 1.5841806 , -0.3022984 },
719  { 1.6128996 , 1.6842538 , -0.2393057 },
720  { 1.8855491 , 1.7709125 , -0.2115213 }},
721  { /* these are S to P */
722  { -1.4293664 , 2.3294080 , -0.0890470 },
723  { -0.3608082 , 2.3337636 , -0.0712380 },
724  { 0.3027974 , 2.3326252 , -0.0579008 },
725  { 0.7841193 , 2.3320138 , -0.0497094 }},
726  { /* these are D to P */
727  { 1.1341403 , 3.1702435 , -0.2085843 },
728  { 1.7915926 , 2.4942946 , -0.2266493 },
729  { 2.1979400 , 2.2785377 , -0.1518743 },
730  { 2.5018229 , 2.1925720 , -0.1081966 }},
731  { /* these are P to D */
732  { 0.0000000 , 0.0000000 , 0.0000000 },
733  { -2.6737396 , 2.9379143 , -0.3805367 },
734  { -1.4380124 , 2.7756396 , -0.2754625 },
735  { -0.6630196 , 2.6887253 , -0.2216493 }},
736  },
737  /* these are for triplets */
738  {
739  { /* these are P to S */
740  { 0.3075287 , 0.9087130 , -1.0387207 },
741  { 0.687069 , 1.1485864 , -0.6627317 },
742  { 0.9776064 , 1.3382024 , -0.5331906 },
743  { 1.2107725 , 1.4943721 , -0.4779232 }},
744  { /* these are S to P */
745  { -1.3659605 , 2.3262253 , -0.0306439 },
746  { -0.2899490 , 2.3279391 , -0.0298695 },
747  { 0.3678878 , 2.3266603 , -0.0240021 },
748  { 0.8427457 , 2.3249540 , -0.0194091 }},
749  { /* these are D to P */
750  { 1.3108281 , 2.8446367 , -0.1649923 },
751  { 1.8437692 , 2.2399326 , -0.2583398 },
752  { 2.1820792 , 2.0693762 , -0.1864091 },
753  { 2.4414052 , 2.0168255 , -0.1426083 }},
754  { /* these are P to D */
755  { 0.0000000 , 0.0000000 , 0.0000000 },
756  { -1.9219877 , 2.7689624 , -0.2536072 },
757  { -0.7818065 , 2.6595150 , -0.1895313 },
758  { -0.0665624 , 2.5955623 , -0.1522616 }},
759  }
760  };
761 
762  if( lLo==0 )
763  {
764  paramSet = 0;
765  }
766  else if( lLo==1 && lHi==0 )
767  {
768  paramSet = 1;
769  }
770  else if( lLo==1 && lHi==2 )
771  {
772  paramSet = 2;
773  }
774  else if( lLo==2 )
775  {
776  paramSet = 3;
777  ASSERT( lHi==1 );
778  }
779 
780  ASSERT( (int)((sHi-1)/2) == 0 || (int)((sHi-1)/2) == 1 );
781  a = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][0];
782  b = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][1];
783  c = extrapol_Params[(int)((sHi-1)/2)][paramSet][nLo-2][2];
784  ASSERT( Enerwn > 0. );
785  x = log( iso_sp[ipHE_LIKE][nelem].fb[ipLo].xIsoLevNIonRyd*RYD_INF/Enerwn );
786 
787  emisOscStr = exp(a+b*x+c*x*x)/pow3(Eff_nupper)*
788  (2.*lLo+1)/(2.*lHi+1);
789 
790  Aul = TRANS_PROB_CONST*Enerwn*Enerwn*emisOscStr;
791 
792  if( nLo==2 && lLo==1 && sLo==3 )
793  {
794  ASSERT( jLo == ipLo-3 );
795  Aul *= (2.*jLo+1.0)/9.0;
796  }
797 
798  ASSERT( Aul > 0. );
799  *error1 = 0.0002f;
800  *error2 = 0.002f;
801  }
802  else
803  {
804  /* Calculate the radial integral from the quantum defects. */
805  RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(ipHELIUM));
806  ASSERT( RI2 > 0. );
807  /* Convert radial integral to Aul. */
808  Aul = ritoa(lHi,lLo,ipHELIUM,Enerwn,RI2);
809  /* radial integral routine does not recognize fine structure.
810  * Here we split 2^3P. */
811  if( nLo==2 && lLo==1 && sLo==3 )
812  {
813  ASSERT( jLo == ipLo-3 );
814  Aul *= (2.*jLo+1.0)/9.0;
815  }
816 
817  ASSERT( Aul > 0. );
818  *error1 = 0.01f;
819  *error2 = 0.07f;
820  }
821  }
822  }
823 
824  /* Heavier species */
825  else
826  {
827  /* Retrieve transition probabilities for Helike ions. */
828  /* >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
829  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J, originally astro.ph. 0201454 */
830  if( ipHi <= MAX_TP_INDEX && nHi <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
831  {
832  /*Must not be accessed by collapsed levels! */
833  ASSERT( ipHi < ( iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) );
834  ASSERT( ipLo < ( iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max ) );
835  ASSERT( ipHi > 2 );
836 
837  Aul = TransProbs[nelem][ipHi][ipLo];
838  *error1 = 0.1f;
839  *error2 = 0.1f;
840  }
841  else
842  {
843  ASSERT( nHi > 2 );
844  long ipHiTemp = iso_sp[ipHE_LIKE][nelem].IndexIfAllResolved[nHi][lHi][sHi];
845  long ipLoTemp = iso_sp[ipHE_LIKE][nelem].IndexIfAllResolved[nLo][lLo][sLo];
846  if( (nLo==2) && (lLo==1) && (sLo==3) )
847  {
848  ASSERT( (jLo>=0) && (jLo<=2) );
849  ipLoTemp -= (2 - jLo);
850  }
851  ASSERT( ipLoTemp < ipHiTemp );
852 
853  if( ipHiTemp <= MAX_TP_INDEX )
854  {
855  Aul = TransProbs[nelem][ipHiTemp][ipLoTemp];
856  *error1 = 0.1f;
857  *error2 = 0.1f;
858  }
859  }
860 
861  if( Aul < 0. )
862  {
863  /* Do same-n transitions. */
864  if( nLo==nHi && lHi<=2 && lLo<=2 )
865  {
866  /* These are 2p3Pj to 2s3S fits to (low Z) Porquet & Dubau (2000) &
867  * (high Z) NIST Atomic Spectra Database. */
868  if( nLo==2 && lLo==0 && sLo==3 )
869  {
870  ASSERT( nHi==2 && lHi==1 && sHi==3 );
871  ASSERT( jHi==ipHi-3 );
872 
873  if( jHi==0 )
874  Aul = 3.31E7 + 1.13E6 * pow((double)nelem+1.,1.76);
875  else if( jHi==1 )
876  Aul = 2.73E7 + 1.31E6 * pow((double)nelem+1.,1.76);
877  else if( jHi==2 )
878  Aul = 3.68E7 + 1.04E7 * exp(((double)nelem+1.)/5.29);
879  else
880  {
881  fprintf(ioQQQ," PROBLEM Impossible value for jHi in he_1trans.\n");
882  TotalInsanity();
883  }
884  }
885 
886  /* These are 2p1P to 2s1S fits to data from TOPbase. */
887  else if( ( nLo==2 && lLo==0 && sLo==1 ) && ( nHi==2 && lHi==1 && sHi==1 ) )
888  {
889  Aul = 5.53E6 * exp( 0.171*(nelem+1.) );
890  }
891 
892  else
893  {
894  /* This case should only be entered if n > 2. Those cases were done above. */
895  ASSERT( nLo > 2 );
896 
897  /* Triplet P to triplet S, delta n = 0 */
898  if( (lHi == 1) && (sHi == 3) && (lLo == 0) && (sLo == 3))
899  {
900  Aul = 0.4 * 3.85E8 * pow((double)nelem,1.6)/pow((double)nHi,5.328);
901  }
902  /* Singlet P to singlet D, delta n = 0 */
903  else if( (lHi == 1) && (sHi == 1) && (lLo == 2) && (sLo == 1))
904  {
905  Aul = 1.95E4 * pow((double)nelem,1.6) / pow((double)nHi, 4.269);
906  }
907  /* Singlet P to singlet S, delta n = 0 */
908  else if( (lHi == 1) && (sHi == 1) && (lLo == 0) )
909  {
910  Aul = 6.646E7 * powpq((double)nelem,3,2) / pow((double)nHi, 5.077);
911  }
912  else
913  {
914  ASSERT( (lHi == 2) && (sHi == 3) && (lLo == 1) );
915  Aul = 3.9E6 * pow((double)nelem,1.6) / pow((double)nHi, 4.9);
916  if( (lHi >2) || (lLo > 2) )
917  Aul *= (lHi/2.);
918  if(lLo > 2)
919  Aul *= (1./9.);
920  }
921  }
922  ASSERT( Aul > 0.);
923  }
924 
925  /* assume transitions involving F and higher orbitals are hydrogenic. */
926  else if( (nHi > nLo) && ((lHi > 2) || (lLo > 2)) )
927  {
928  Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem);
929  ASSERT( Aul > 0.);
930  }
931 
932  /* These transitions are of great importance, but the below radial integral
933  * routine fails to achieve desirable accuracy, so these are fits as produced
934  * from He A's for nupper through 9. They are transitions to ground and
935  * 2, 3, and 4 triplet S. */
936  else if( ( nLo==1 ) || ( nLo==2 && lLo==0 ) ||
937  ( nLo==3 && lLo==0 && sLo==3 ) ||
938  ( nLo==4 && lLo==0 && sLo==3 ) )
939  {
940  /* Here are the Lyman transitions. */
941  if( nLo==1 )
942  {
943  ASSERT( (lHi == 1) && (sHi == 1) );
944 
945  /* In theory, this Z dependence should be Z^4, but values from TOPbase
946  * suggest 3.9 is a more accurate exponent. Values from
947  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
948  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
949  /* originally astro.ph. 0201454 */
950  Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
951  }
952 
953  /* Here are the Balmer transitions. */
954  else if( nLo==2 && lLo==0 && sLo==1 )
955  {
956  ASSERT( (lHi == 1) && (sHi == 1) );
957 
958  /* More fits, as above. */
959  Aul = 5.0e8 * pow4((double)nelem) / pow((double)nHi, 2.889);
960  }
961 
962  /* Here are transitions down to triplet S */
963  else
964  {
965  ASSERT( (lLo==0) && (sLo==3) );
966  ASSERT( (lHi==1) && (sHi==3) );
967 
968  /* These fits also as above. */
969  if( nLo == 2 )
970  Aul = 1.5 * 3.405E8 * pow4((double)nelem) / pow((double)nHi, 2.883);
971  else if( nLo == 3 )
972  Aul = 2.5 * 4.613E7 * pow4((double)nelem) / pow((double)nHi, 2.672);
973  else
974  Aul = 3.0 * 1.436E7 * pow4((double)nelem) / pow((double)nHi, 2.617);
975  }
976 
977  ASSERT( Aul > 0.);
978  }
979 
980  /* Every other allowed transition is calculated as follows. */
981  else
982  {
983  /* Calculate the radial integral from the quantum defects. */
984  RI2 = scqdri(Eff_nupper,lHi,Eff_nlower,lLo,(double)(nelem));
985  /* Convert radial integral to Aul. */
986  Aul = ritoa(lHi,lLo,nelem,Enerwn,RI2);
987  /* radial integral routine does not recognize fine structure.
988  * Here we split 2^3P. */
989  if( nLo==2 && lLo==1 && sLo==3 && Aul > iso_ctrl.SmallA )
990  {
991  ASSERT( jLo == ipLo-3 );
992  Aul *= (2.*jLo+1.0)/9.0;
993  }
994 
995  }
996  *error1 = 0.1f;
997  *error2 = 0.1f;
998  }
999  }
1000  }
1001 
1002  /* Now do forbidden transitions from 2-1 ... */
1003  /* and those given by
1004  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
1005  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
1006  /* originally astro.ph. 0201454
1007  * for heavy elements. These are triplet P to singlet S,
1008  * going either up or down...Triplet S to Singlet P are not included, as they are far weaker. */
1009  else
1010  {
1011  ASSERT( (sHi != sLo) || (abs((int)(lHi - lLo)) != 1) );
1012  Aul = ForbiddenAuls(ipHi, ipLo, nelem);
1013  ASSERT( Aul > 0. );
1014  }
1015 
1016  Aul = MAX2( Aul, iso_ctrl.SmallA );
1017  ASSERT( Aul >= iso_ctrl.SmallA );
1018 
1019  /* negative energy for a transition with substantial transition probability
1020  * would be major logical error - but is ok for same-n l transitions */
1021  if( Enerwn < 0. && Aul > iso_ctrl.SmallA )
1022  {
1023  fprintf( ioQQQ," he_1trans hit negative energy, nelem=%li, val was %f \n", nelem ,Enerwn );
1024  }
1025 
1026  return Aul;
1027 }
1028 
1029 STATIC double helike_transprob_collapsed_to_collapsed( long nelem, long nHi, long nLo, double Enerwn )
1030 {
1031  DEBUG_ENTRY( "helike_transprob_collapsed_to_collapsed()" );
1032 
1033  double Aul = 0.;
1034  ASSERT( Enerwn > 0. );
1035  ASSERT( nHi > nLo );
1036 
1037  for( long lLo=0; lLo < nLo; ++lLo )
1038  {
1039  double Aul_sing = helike_transprob_collapsed_to_resolved( nelem, nHi, nLo, lLo, 1, -1, Enerwn );
1040  double Aul_trip = helike_transprob_collapsed_to_resolved( nelem, nHi, nLo, lLo, 3, -1, Enerwn );
1041  Aul += Aul_sing + Aul_trip;
1042  }
1043 
1044  return Aul;
1045 }
1046 
1047 STATIC double helike_transprob_collapsed_to_resolved( long nelem, long nHi, long nLo, long lLo, long sLo, long jLo, double Enerwn )
1048 {
1049  DEBUG_ENTRY( "helike_transprob_collapsed_to_resolved()" );
1050 
1051  double n_eff_hi = nHi - helike_quantum_defect(nelem,nHi,-1,-1,-1);
1052  double n_eff_lo = nLo - helike_quantum_defect(nelem,nLo,lLo,sLo,jLo);
1053  realnum error1, error2;
1054 
1055  long ipISO = ipHE_LIKE;
1056  t_iso_sp *sp = &iso_sp[ipISO][nelem];
1057  long nResMax = sp->n_HighestResolved_max;
1058  long ipLoRes = sp->IndexIfAllResolved[nLo][lLo][sLo];
1059  // make sure 2^3Pj terms are fully specified
1060  if( nLo==2 && lLo==1 && sLo==3 )
1061  {
1062  ASSERT( jLo>=0 && jLo<=2 );
1063  ipLoRes -= (2 - jLo);
1064  }
1065  ASSERT( nLo > nResMax || ipLoRes <= sp->numLevels_max - sp->nCollapsed_max );
1066 
1067  /* Lower level resolved, upper not. First calculate Aul
1068  * from upper level with ang mom one higher. */
1069  double Aul = he_1trans( nelem, Enerwn,
1070  n_eff_hi, lLo+1, sLo, -1, n_eff_lo, lLo, sLo, jLo,
1071  &error1, &error2);
1072 
1073  sp->CachedAs[ nHi-nResMax-1 ][ ipLoRes ][0] = (realnum)Aul;
1074 
1075  Aul *= (2.*(lLo+1.)+1.) * sLo / (4.*(double)nHi*(double)nHi);
1076 
1077  if( lLo != 0 )
1078  {
1079  /* For all l>0, add in transitions from upper level with ang mom one lower. */
1080  double Aul1 = he_1trans( nelem, Enerwn,
1081  n_eff_hi, lLo-1, sLo, -1, n_eff_lo, lLo, sLo, jLo,
1082  &error1, &error2);
1083 
1084  sp->CachedAs[ nHi-nResMax-1 ][ ipLoRes ][1] = (realnum)Aul1;
1085 
1086  /* now add in this part after multiplying by stat weight for lHi = lLo-1. */
1087  Aul += Aul1*(2.*(lLo-1.)+1.) * sLo / (4.*(double)nHi*(double)nHi);
1088  }
1089  else
1090  sp->CachedAs[ nHi-nResMax-1 ][ ipLoRes ][1] = 0.f;
1091 
1092  ASSERT( Aul > 0.);
1093 
1094  return Aul;
1095 }
1096 
1097 void DoFSMixing( long nelem, long ipLoSing, long ipHiSing )
1098 {
1099  /* Every bit of this routine is based upon the singlet-triplet mixing formalism given in
1100  * >>refer He FSM Drake, G. W. F. 1996, in Atomic, Molecular, \& Optical Physics Handbook,
1101  * >>refercon ed. G. W. F. Drake (New York: AIP Press).
1102  * That formalism mixes the levels themselves, but since this code is not J-resolved, we simulate
1103  * that by mixing only the transition probabilities. We find results comparable to those calculated
1104  * in the fully J-resolved model spearheaded by Rob Bauman, and described in
1105  * >>refer He FSM Bauman, R. P., Porter, R. L., Ferland, G. J., \& MacAdam, K. B. 2005, ApJ, accepted */
1106  long int nHi, lHi, sHi, nLo, lLo, sLo, ipHiTrip, ipLoTrip;
1107  double Ass, Att, Ast, Ats;
1108  double SinHi, SinLo, CosHi, CosLo;
1109  double HiMixingAngle, LoMixingAngle , error;
1110  double Kss, Ktt, Kts, Kst, fss, ftt, fssNew, fttNew, ftsNew, fstNew, temp;
1111 
1112  DEBUG_ENTRY( "DoFSMixing()" );
1113 
1114  nHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].n();
1115  lHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].l();
1116  sHi = iso_sp[ipHE_LIKE][nelem].st[ipHiSing].S();
1117  nLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].n();
1118  lLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].l();
1119  sLo = iso_sp[ipHE_LIKE][nelem].st[ipLoSing].S();
1120 
1121  if( ( sHi == 3 || sLo == 3 ) ||
1122  ( abs(lHi - lLo) != 1 ) ||
1123  ( nLo < 2 ) ||
1124  ( lHi <= 1 || lLo <= 1 ) ||
1125  ( nHi == nLo && lHi == 1 && lLo == 2 ) ||
1126  ( nHi > nLo && lHi != 1 && lLo != 1 ) )
1127  {
1128  return;
1129  }
1130 
1131  ASSERT( lHi > 0 );
1132  /*ASSERT( (lHi > 1) && (lLo > 1) );*/
1133 
1134  ipHiTrip = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][3];
1135  ipLoTrip = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][3];
1136 
1137  if( lHi == 2 )
1138  {
1139  HiMixingAngle = 0.01;
1140  }
1141  else if( lHi == 3 )
1142  {
1143  HiMixingAngle = 0.5;
1144  }
1145  else
1146  {
1147  HiMixingAngle = PI/4.;
1148  }
1149 
1150  if( lLo == 2 )
1151  {
1152  LoMixingAngle = 0.01;
1153  }
1154  else if( lLo == 3 )
1155  {
1156  LoMixingAngle = 0.5;
1157  }
1158  else
1159  {
1160  LoMixingAngle = PI/4.;
1161  }
1162 
1163  /* These would not work correctly if l<=1 were included in this treatment! */
1164  ASSERT( ipHiTrip > ipLoTrip );
1165  ASSERT( ipHiTrip > ipLoSing );
1166  ASSERT( ipHiSing > ipLoTrip );
1167  ASSERT( ipHiSing > ipLoSing );
1168 
1169  SinHi = sin( HiMixingAngle );
1170  SinLo = sin( LoMixingAngle );
1171  CosHi = cos( HiMixingAngle );
1172  CosLo = cos( LoMixingAngle );
1173 
1174  Kss = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).EnergyWN();
1175  Ktt = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).EnergyWN();
1176  Kst = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).EnergyWN();
1177  Kts = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).EnergyWN();
1178 
1179  fss = iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()/TRANS_PROB_CONST/POW2(Kss)/
1180  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Lo()).g()*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Hi()).g();
1181 
1182  ftt = iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul()/TRANS_PROB_CONST/POW2(Ktt)/
1183  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Lo()).g()*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Hi()).g();
1184 
1185  temp = sqrt(fss/Kss)*CosHi*CosLo + sqrt(ftt/Ktt)*SinHi*SinLo;
1186  fssNew = Kss*POW2( temp );
1187  temp = sqrt(fss/Kss)*SinHi*SinLo + sqrt(ftt/Ktt)*CosHi*CosLo;
1188  fttNew = Ktt*POW2( temp );
1189  temp = sqrt(fss/Kss)*CosHi*SinLo - sqrt(ftt/Ktt)*SinHi*CosLo;
1190  fstNew = Kst*POW2( temp );
1191  temp = sqrt(fss/Kss)*SinHi*CosLo - sqrt(ftt/Ktt)*CosHi*SinLo;
1192  ftsNew = Kts*POW2( temp );
1193 
1194  Ass = TRANS_PROB_CONST*POW2(Kss)*fssNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Lo()).g()/
1195  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Hi()).g();
1196 
1197  Att = TRANS_PROB_CONST*POW2(Ktt)*fttNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Lo()).g()/
1198  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Hi()).g();
1199 
1200  Ast = TRANS_PROB_CONST*POW2(Kst)*fstNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Lo()).g()/
1201  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Hi()).g();
1202 
1203  Ats = TRANS_PROB_CONST*POW2(Kts)*ftsNew*(*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Lo()).g()/
1204  (*iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Hi()).g();
1205 
1206  error = fabs( ( iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul()+
1207  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul() )/
1208  (Ass+Ast+Ats+Att) - 1.f );
1209 
1210  if( error > 0.001 )
1211  {
1212  fprintf( ioQQQ, "FSM error: %e\t Element: %ld\t States: LS %li HS %li LT %li HT %li\t New values: Ass %e Att %e Ast %e Ats %e\n", error,
1213  nelem, ipLoSing, ipHiSing, ipLoTrip, ipHiTrip, Ass, Att, Ast, Ats );
1214  }
1215 
1216  if( !iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).hasEmis() )
1217  {
1218  iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).AddLine2Stack();
1219  iso_setRedisFun(ipHE_LIKE, nelem, ipLoSing, ipHiSing);
1220  iso_setOpacity (ipHE_LIKE, nelem, ipLoSing, ipHiSing);
1221  }
1222  if( !iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).hasEmis() )
1223  {
1224  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).AddLine2Stack();
1225  iso_setRedisFun(ipHE_LIKE, nelem, ipLoTrip, ipHiTrip);
1226  iso_setOpacity (ipHE_LIKE, nelem, ipLoTrip, ipHiTrip);
1227  }
1228  if( !iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).hasEmis() )
1229  {
1230  iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).AddLine2Stack();
1231  iso_setRedisFun(ipHE_LIKE, nelem, ipLoTrip, ipHiSing);
1232  iso_setOpacity (ipHE_LIKE, nelem, ipLoTrip, ipHiSing);
1233  }
1234  if( !iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).hasEmis() )
1235  {
1236  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).AddLine2Stack();
1237  iso_setRedisFun(ipHE_LIKE, nelem, ipLoSing, ipHiTrip);
1238  iso_setOpacity (ipHE_LIKE, nelem, ipLoSing, ipHiTrip);
1239  }
1240 
1241  iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoSing).Emis().Aul() = (realnum)Ass;
1242  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoTrip).Emis().Aul() = (realnum)Att;
1243  iso_sp[ipHE_LIKE][nelem].trans(ipHiSing,ipLoTrip).Emis().Aul() = (realnum)Ast;
1244  iso_sp[ipHE_LIKE][nelem].trans(ipHiTrip,ipLoSing).Emis().Aul() = (realnum)Ats;
1245 
1246  return;
1247 }
1248 
1249 /*ritoa converts the square of the radial integral for a transition
1250  * (calculated by scqdri) to the transition probability, Aul. */
1251 STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
1252 {
1253  /* Variables are as follows: */
1254  /* lg = larger of li and lf */
1255  /* fmean = mean oscillator strength */
1256  /* for a given level. */
1257  /* mu = reduced mass of optical electron. */
1258  /* EinsteinA = Einstein emission coef. */
1259  /* w = angular frequency of transition. */
1260  /* RI2_cm = square of rad. int. in cm^2. */
1261  long lg;
1262  double fmean,mu,EinsteinA,w,RI2_cm;
1263 
1264  DEBUG_ENTRY( "ritoa()" );
1265 
1266  mu = ELECTRON_MASS/(1+ELECTRON_MASS/(dense.AtomicWeight[nelem]*ATOMIC_MASS_UNIT));
1267 
1268  w = 2. * PI * k * SPEEDLIGHT;
1269 
1270  RI2_cm = RI2 * BOHR_RADIUS_CM * BOHR_RADIUS_CM;
1271 
1272  lg = max(lf, li);
1273 
1274  fmean = 2.0*mu*w*lg*RI2_cm/((3.0*H_BAR) * (2.0*li+1.0));
1275 
1276  EinsteinA = TRANS_PROB_CONST*k*k*fmean;
1277 
1278  /* ASSERT( EinsteinA > SMALLFLOAT ); */
1279 
1280  return EinsteinA;
1281 }
1282 
1283 realnum helike_transprob( long nelem, long ipHi, long ipLo )
1284 {
1285  double Aul;
1286  double Enerwn, n_eff_hi, n_eff_lo;
1287  long ipISO = ipHE_LIKE;
1288 
1289  DEBUG_ENTRY( "helike_transprob()" );
1290 
1291  Enerwn = iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN();
1292  n_eff_hi = N_(ipHi) - helike_quantum_defect(nelem,N_(ipHi),L_(ipHi),S_(ipHi),J_(ipHi));
1293  n_eff_lo = N_(ipLo) - helike_quantum_defect(nelem,N_(ipLo),L_(ipLo),S_(ipLo),J_(ipLo));
1294 
1295  if( ipHi >= iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max )
1296  {
1297  if( ipLo >= iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max )
1298  {
1299  /* Neither upper nor lower is resolved...Aul()'s are easy. */
1300  Aul = helike_transprob_collapsed_to_collapsed( nelem, N_(ipHi), N_(ipLo), Enerwn );
1301  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f,0.001f);
1302 
1303  ASSERT( Aul > 0.);
1304  }
1305  else
1306  {
1307  Aul = helike_transprob_collapsed_to_resolved( nelem, N_(ipHi), N_(ipLo), L_(ipLo), S_(ipLo), J_(ipLo), Enerwn );
1308  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.01f);
1309  ASSERT( Aul > 0.);
1310  }
1311  }
1312  else
1313  {
1314  realnum error1, error2;
1315  /* Both levels are resolved...do the normal bit. */
1316  if( Enerwn < 0. )
1317  {
1318  Aul = he_1trans( nelem, -1.*Enerwn, n_eff_lo,
1319  L_(ipLo), S_(ipLo), ipLo-3, n_eff_hi, L_(ipHi), S_(ipHi), ipHi-3, &error1, &error2);
1320  }
1321  else
1322  {
1323  Aul = he_1trans( nelem, Enerwn, n_eff_hi,
1324  L_(ipHi), S_(ipHi), ipHi-3, n_eff_lo, L_(ipLo), S_(ipLo), ipLo-3, &error1, &error2);
1325  }
1326  iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.01f);
1327  }
1328 
1329  return (realnum)Aul;
1330 }
1331 
1332 /* This routine is pure infrastructure and bookkeeping, no physics. */
1333 void HelikeTransProbSetup( void )
1334 {
1335 
1336  const int chLine_LENGTH = 1000;
1337  char chLine[chLine_LENGTH];
1338 
1339  FILE *ioDATA;
1340  bool lgEOL;
1341 
1342  long nelem, ipLo, ipHi, i, i1, i2, i3;
1343 
1344  DEBUG_ENTRY( "HelikeTransProbSetup()" );
1345 
1346  TransProbs = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
1347 
1348  for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
1349  {
1350 
1351  TransProbs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MAX_TP_INDEX+1) );
1352 
1353  for( ipLo=ipHe1s1S; ipLo <= MAX_TP_INDEX;++ipLo )
1354  {
1355  TransProbs[nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)MAX_TP_INDEX );
1356  }
1357  }
1358 
1359  /********************************************************************/
1360  /*************** Read in data from he_transprob.dat *****************/
1361  if( trace.lgTrace )
1362  fprintf( ioQQQ," HelikeTransProbSetup opening he_transprob.dat:");
1363 
1364  ioDATA = open_data( "he_transprob.dat", "r" );
1365 
1366  /* check that magic number is ok */
1367  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1368  {
1369  fprintf( ioQQQ, " HelikeTransProbSetup could not read first line of he_transprob.dat.\n");
1370  cdEXIT(EXIT_FAILURE);
1371  }
1372  i = 1;
1373  i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1374  i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1375  if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
1376  {
1377  fprintf( ioQQQ,
1378  " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1379  fprintf( ioQQQ,
1380  " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1381  TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
1382  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1383  cdEXIT(EXIT_FAILURE);
1384  }
1385 
1386  /* Initialize TransProbs[nelem][ipLo][ipHi] before filling it in. */
1387  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
1388  {
1389  for( ipHi=0; ipHi <= MAX_TP_INDEX; ipHi++ )
1390  {
1391  for( ipLo=0; ipLo < MAX_TP_INDEX; ipLo++ )
1392  {
1393  TransProbs[nelem][ipHi][ipLo] = -1.;
1394  }
1395  }
1396  }
1397 
1398  for( ipLo=1; ipLo <= N_HE1_TRANS_PROB; ipLo++ )
1399  {
1400  char *chTemp;
1401 
1402  /* get next line image */
1403  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1404  BadRead();
1405 
1406  while( chLine[0]=='#' )
1407  {
1408  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1409  BadRead();
1410  }
1411 
1412  i3 = 1;
1413  i1 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
1414  i2 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
1415  /* check that these numbers are correct */
1416  if( i1<0 || i2<=i1 )
1417  {
1418  fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1419  cdEXIT(EXIT_FAILURE);
1420  }
1421 
1422  chTemp = chLine;
1423 
1424  /* skip over 2 tabs to start of data */
1425  for( i=0; i<1; ++i )
1426  {
1427  if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
1428  {
1429  fprintf( ioQQQ, " HelikeTransProbSetup could not init he_transprob\n" );
1430  cdEXIT(EXIT_FAILURE);
1431  }
1432  ++chTemp;
1433  }
1434 
1435  /* now read in the data */
1436  for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
1437  {
1438  if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
1439  {
1440  fprintf( ioQQQ, " HelikeTransProbSetup could not scan he_transprob\n" );
1441  cdEXIT(EXIT_FAILURE);
1442  }
1443  ++chTemp;
1444 
1445  sscanf( chTemp , "%le" , &TransProbs[nelem][i2][i1] );
1446 
1448  if( lgEOL )
1449  {
1450  fprintf( ioQQQ, " HelikeTransProbSetup detected insanity in he_transprob.dat.\n");
1451  cdEXIT(EXIT_FAILURE);
1452  }
1453  }
1454  }
1455 
1456  /* check that ending magic number is ok */
1457  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1458  {
1459  fprintf( ioQQQ, " HelikeTransProbSetup could not read last line of he_transprob.dat.\n");
1460  cdEXIT(EXIT_FAILURE);
1461  }
1462  i = 1;
1463  i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1464  i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1465  if( i1 !=TRANSPROBMAGIC || i2 != N_HE1_TRANS_PROB )
1466  {
1467  fprintf( ioQQQ,
1468  " HelikeTransProbSetup: the version of he_transprob.dat is not the current version.\n" );
1469  fprintf( ioQQQ,
1470  " HelikeTransProbSetup: I expected to find the number %i %i and got %li %li instead.\n" ,
1471  TRANSPROBMAGIC, N_HE1_TRANS_PROB, i1, i2 );
1472  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1473  cdEXIT(EXIT_FAILURE);
1474  }
1475 
1476  /* close the data file */
1477  fclose( ioDATA );
1478  return;
1479 }
pow4
T pow4(T a)
Definition: cddefines.h:995
ForbiddenAuls
STATIC double ForbiddenAuls(long ipHi, long ipLo, long nelem)
Definition: helike_einsta.cpp:304
Jint
STATIC double Jint(double theta)
Definition: helike_einsta.cpp:46
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
t_iso_sp::st
qList st
Definition: iso.h:482
bessjnu
STATIC double bessjnu(double vv, double zz)
Definition: helike_einsta.cpp:90
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:65
J_
#define J_(A_)
Definition: iso.h:25
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
MAX_TP_INDEX
#define MAX_TP_INDEX
Definition: helike_einsta.h:9
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:596
ipARGON
const int ipARGON
Definition: cddefines.h:366
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:49
eina
double eina(double gf, double enercm, double gup)
Definition: lines_service.cpp:75
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:48
helike_transprob_collapsed_to_resolved
STATIC double helike_transprob_collapsed_to_resolved(long nelem, long nHi, long nLo, long lLo, long sLo, long jLo, double Enerwn)
Definition: helike_einsta.cpp:1047
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:46
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:518
helike_einsta.h
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:51
ipHe3p3P
const int ipHe3p3P
Definition: iso.h:56
pow3
T pow3(T a)
Definition: cddefines.h:988
vJint
static double vJint
Definition: helike_einsta.cpp:43
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:43
iso_put_error
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
AngerJ
STATIC double AngerJ(double vv, double zz)
Definition: helike_einsta.cpp:127
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
IPRAD
#define IPRAD
Definition: iso.h:88
t_iso_sp::IndexIfAllResolved
multi_arr< long, 3 > IndexIfAllResolved
Definition: iso.h:492
dense
t_dense dense
Definition: global.cpp:15
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:391
trace
t_trace trace
Definition: trace.cpp:5
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:554
lines_service.h
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, realnum *error1, realnum *error2)
Definition: helike_einsta.cpp:584
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:536
ipHe2s1S
const int ipHe2s1S
Definition: iso.h:47
L_
#define L_(A_)
Definition: iso.h:23
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:477
helike_quantum_defect
double helike_quantum_defect(long nelem, long n, long lqn, long S, long j)
Definition: helike_energy.cpp:245
POW2
#define POW2
Definition: cddefines.h:979
STATIC
#define STATIC
Definition: cddefines.h:118
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:447
N_
#define N_(A_)
Definition: iso.h:22
ipHe3p1P
const int ipHe3p1P
Definition: iso.h:59
realnum
float realnum
Definition: cddefines.h:124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
scqdri
STATIC double scqdri(double nstar, long int l, double npstar, long int lp, double iz)
Definition: helike_einsta.cpp:215
max
long max(int a, long b)
Definition: cddefines.h:817
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:435
ipHe3s3S
const int ipHe3s3S
Definition: iso.h:54
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
ipNEON
const int ipNEON
Definition: cddefines.h:358
S_
#define S_(A_)
Definition: iso.h:24
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1283
strchr_s
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1310
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
zJint
static double zJint
Definition: helike_einsta.cpp:43
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
t_iso_sp
Definition: iso.h:470
helike_transprob_collapsed_to_collapsed
STATIC double helike_transprob_collapsed_to_collapsed(long nelem, long nHi, long nLo, double Enerwn)
Definition: helike_einsta.cpp:1029
hydro_bauman.h
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:431
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:1097
qg32
double qg32(double, double, double(*)(double))
Definition: service.cpp:1175
LIMELM
const int LIMELM
Definition: cddefines.h:308
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:50
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:630
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:350
iso_setRedisFun
void iso_setRedisFun(long ipISO, long nelem, long ipLo, long ipHi)
Definition: iso_create.cpp:48
AngerJ_mac
STATIC double AngerJ_mac(double vv, double zz)
Definition: helike_einsta.cpp:60
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
chLine_LENGTH
#define chLine_LENGTH
ritoa
STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
Definition: helike_einsta.cpp:1251
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1826
N_HE1_TRANS_PROB
#define N_HE1_TRANS_PROB
Definition: helike_einsta.h:7
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
fixit
#define fixit(a)
Definition: cddefines.h:417
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1333
iso_setOpacity
void iso_setOpacity(long ipISO, long nelem, long ipLo, long ipHi)
Definition: iso_create.cpp:82
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:349
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:690
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:2356
AngerJ_asymp
STATIC double AngerJ_asymp(double vv, double zz)
Definition: helike_einsta.cpp:108
BadRead
NORETURN void BadRead(void)
Definition: service.cpp:986
ipHe3s1S
const int ipHe3s1S
Definition: iso.h:55
TRANSPROBMAGIC
#define TRANSPROBMAGIC
Definition: helike.h:22
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:631
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
TransProbs
static double *** TransProbs
Definition: helike_einsta.cpp:25
Generated by   doxygen 1.8.5