cloudy trunk
Loading...
Searching...
No Matches
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. */
25static double ***TransProbs;
26
27/*ritoa converts the square of the radial integral for a transition
28 * (calculated by scqdri) to the transition probability, Aul. */
29STATIC double ritoa( long li, long lf, long nelem, double k, double RI2 );
30
31/* handles all forbidden transition probabilities. */
32STATIC double ForbiddenAuls( long ipHi, long ipLo, long nelem );
33
34/*
35static long RI_ipHi, RI_ipLo, RI_ipLev;
36static long globalZ;
37*/
38
39/* used as parameters in qg32 integration */
40static double vJint , zJint;
41
42/* the integrand used in the AngerJ function below. */
43STATIC 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. */
58STATIC 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/******************************************************************************/
146double 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
221STATIC 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. */
501double 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
894void 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. */
1027STATIC 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
1059realnum 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. */
1132void 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}
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
STATIC
#define STATIC
Definition: cddefines.h:97
MIN2
#define MIN2
Definition: cddefines.h:761
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
qg32
double qg32(double, double, double(*)(double))
Definition: service.cpp:1053
BadRead
NORETURN void BadRead(void)
Definition: service.cpp:901
LIMELM
const int LIMELM
Definition: cddefines.h:258
POW2
#define POW2
Definition: cddefines.h:929
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
realnum
float realnum
Definition: cddefines.h:103
MAX2
#define MAX2
Definition: cddefines.h:782
strchr_s
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1439
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipNEON
const int ipNEON
Definition: cddefines.h:314
ipARGON
const int ipARGON
Definition: cddefines.h:322
fixit
void fixit(void)
Definition: service.cpp:991
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:567
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_iso_sp::st
qList st
Definition: iso.h:453
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
dense
t_dense dense
Definition: dense.cpp:24
TRANSPROBMAGIC
#define TRANSPROBMAGIC
Definition: helike.h:22
helike_quantum_defect
double helike_quantum_defect(long int nelem, long int ipLev)
Definition: helike_energy.cpp:243
chLine_LENGTH
#define chLine_LENGTH
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:894
Jint
STATIC double Jint(double theta)
Definition: helike_einsta.cpp:43
AngerJ
STATIC double AngerJ(double vv, double zz)
Definition: helike_einsta.cpp:58
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1132
vJint
static double vJint
Definition: helike_einsta.cpp:40
scqdri
double scqdri(double nstar, long int l, double npstar, long int lp, double iz)
Definition: helike_einsta.cpp:146
ForbiddenAuls
STATIC double ForbiddenAuls(long ipHi, long ipLo, long nelem)
Definition: helike_einsta.cpp:221
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
TransProbs
static double *** TransProbs
Definition: helike_einsta.cpp:25
zJint
static double zJint
Definition: helike_einsta.cpp:40
ritoa
STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
Definition: helike_einsta.cpp:1027
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1059
helike_einsta.h
MAX_TP_INDEX
#define MAX_TP_INDEX
Definition: helike_einsta.h:9
N_HE1_TRANS_PROB
#define N_HE1_TRANS_PROB
Definition: helike_einsta.h:7
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
hydro_bauman.h
HydroEinstA
double HydroEinstA(long int n1, long int n2)
Definition: hydroeinsta.cpp:12
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
IPRAD
#define IPRAD
Definition: iso.h:86
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:49
ipHe2s1S
const int ipHe2s1S
Definition: iso.h:45
ipHe3p3P
const int ipHe3p3P
Definition: iso.h:54
ipHe3s1S
const int ipHe3s1S
Definition: iso.h:53
N_
#define N_(A_)
Definition: iso.h:20
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
S_
#define S_(A_)
Definition: iso.h:22
ipHe3p1P
const int ipHe3p1P
Definition: iso.h:57
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
ipHe3s3S
const int ipHe3s3S
Definition: iso.h:52
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
L_
#define L_(A_)
Definition: iso.h:21
iso_put_error
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
eina
double eina(double gf, double enercm, double gup)
Definition: lines_service.cpp:84
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
H_BAR
UNUSED const double H_BAR
Definition: physconst.h:144
TRANS_PROB_CONST
UNUSED const double TRANS_PROB_CONST
Definition: physconst.h:237
PI
UNUSED const double PI
Definition: physconst.h:29
BOHR_RADIUS_CM
UNUSED const double BOHR_RADIUS_CM
Definition: physconst.h:222
SPEEDLIGHT
UNUSED const double SPEEDLIGHT
Definition: physconst.h:100
ELECTRON_MASS
UNUSED const double ELECTRON_MASS
Definition: physconst.h:91
RYD_INF
UNUSED const double RYD_INF
Definition: physconst.h:115
ATOMIC_MASS_UNIT
UNUSED const double ATOMIC_MASS_UNIT
Definition: physconst.h:88
g
static double * g
Definition: species2.cpp:28
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
trace
t_trace trace
Definition: trace.cpp:5
Generated by doxygen 1.9.5