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:
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");
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. */
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");
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" ,
1181 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
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");
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" );
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" );
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");
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");
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" ,
1271 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1273 }
1274
1275 /* close the data file */
1276 fclose( ioDATA );
1277 return;
1278}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
double qg32(double, double, double(*)(double))
Definition service.cpp:1053
NORETURN void BadRead(void)
Definition service.cpp:901
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipNEON
Definition cddefines.h:314
const int ipARGON
Definition cddefines.h:322
void fixit(void)
Definition service.cpp:991
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
t_dense dense
Definition dense.cpp:24
#define TRANSPROBMAGIC
Definition helike.h:22
double helike_quantum_defect(long int nelem, long int ipLev)
#define chLine_LENGTH
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
STATIC double Jint(double theta)
STATIC double AngerJ(double vv, double zz)
void HelikeTransProbSetup(void)
static double vJint
double scqdri(double nstar, long int l, double npstar, long int lp, double iz)
STATIC double ForbiddenAuls(long ipHi, long ipLo, long nelem)
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)
static double *** TransProbs
static double zJint
STATIC double ritoa(long li, long lf, long nelem, double k, double RI2)
realnum helike_transprob(long nelem, long ipHi, long ipLo)
#define MAX_TP_INDEX
#define N_HE1_TRANS_PROB
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
double HydroEinstA(long int n1, long int n2)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
#define IPRAD
Definition iso.h:86
const int ipHe2s3S
Definition iso.h:44
const int ipHe2p1P
Definition iso.h:49
const int ipHe2s1S
Definition iso.h:45
const int ipHe3p3P
Definition iso.h:54
const int ipHe3s1S
Definition iso.h:53
#define N_(A_)
Definition iso.h:20
const int ipHE_LIKE
Definition iso.h:63
const int ipHe2p3P1
Definition iso.h:47
#define S_(A_)
Definition iso.h:22
const int ipHe3p1P
Definition iso.h:57
const int ipHe2p3P0
Definition iso.h:46
const int ipHe3s3S
Definition iso.h:52
const int ipHe2p3P2
Definition iso.h:48
#define L_(A_)
Definition iso.h:21
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
double eina(double gf, double enercm, double gup)
#define S(I_, J_)
UNUSED const double H_BAR
Definition physconst.h:144
UNUSED const double TRANS_PROB_CONST
Definition physconst.h:237
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOHR_RADIUS_CM
Definition physconst.h:222
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double RYD_INF
Definition physconst.h:115
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
static double * g
Definition species2.cpp:28
t_trace trace
Definition trace.cpp:5