cloudy trunk
Loading...
Searching...
No Matches
iso_create.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/*iso_create create data for hydrogen and helium, 1 per coreload, called by ContCreatePointers
4 * in turn called after commands parsed */
5/*iso_zero zero data for hydrogen and helium */
6#include "cddefines.h"
7#include "atmdat.h"
8#include "dense.h"
9#include "elementnames.h"
10#include "helike.h"
11#include "helike_einsta.h"
12#include "hydro_bauman.h"
13#include "hydrogenic.h"
14#include "hydroeinsta.h"
15#include "iso.h"
16#include "lines_service.h"
17#include "opacity.h"
18#include "phycon.h"
19#include "physconst.h"
20#include "secondaries.h"
21#include "taulines.h"
22#include "thirdparty.h"
23
24/*iso_zero zero data for hydrogen and helium */
25STATIC void iso_zero(void);
26
27/* allocate memory for iso sequence structures */
28STATIC void iso_allocate(void);
29
30/* define levels of iso sequences and assign quantum numbers to those levels */
31STATIC void iso_assign_quantum_numbers(void);
32
33STATIC void FillExtraLymanLine( const TransitionList::iterator& t, long ipISO, long nelem, long nHi );
34
35STATIC void iso_satellite( void );
36
37char chL[21]={'S','P','D','F','G','H','I','K','L','M','N','O','Q','R','T','U','V','W','X','Y','Z'};
38
39void iso_create(void)
40{
41 long int ipHi,
42 ipLo;
43
44 static int nCalled = 0;
45
46 double HIonPoten;
47
48 DEBUG_ENTRY( "iso_create()" );
49
50 /* > 1 if not first call, then just zero arrays out */
51 if( nCalled > 0 )
52 {
53 iso_zero();
54 return;
55 }
56
57 /* this is first call, increment the nCalled counterso never do this again */
58 ++nCalled;
59
60 /* these are the statistical weights of the ions */
61 iso_ctrl.stat_ion[ipH_LIKE] = 1.f;
62 iso_ctrl.stat_ion[ipHE_LIKE] = 2.f;
63
64 /* this routine allocates all the memory
65 * needed for the iso sequence structures */
66 iso_allocate();
67
68 /* loop over iso sequences and assign quantum numbers to all levels */
69 iso_assign_quantum_numbers();
70
71 /* this is a dummy line, junk it too. */
72 (*TauDummy).Junk();
73 (*TauDummy).AddHiState();
74 (*TauDummy).AddLoState();
75 (*TauDummy).AddLine2Stack();
76
77 /********************************************/
78 /********** Line and level energies ********/
79 /********************************************/
80 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
81 {
82 /* main hydrogenic arrays, fill with sane values */
83 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
84 {
85 /* must always do helium even if turned off */
86 if( nelem < 2 || dense.lgElmtOn[nelem] )
87 {
88 /* Dima's array has ionization potentials in eV, but not on same
89 * scale as cloudy itself*/
90 /* extra factor accounts for this */
91 HIonPoten = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
92 ASSERT(HIonPoten > 0.);
93
94 double EnergyRydGround = 0.;
95 /* go from ground to the highest level */
96 for( ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
97 {
98 double EnergyWN, EnergyRyd;
99
100 if( ipISO == ipH_LIKE )
101 {
102 EnergyRyd = HIonPoten/POW2((double)N_(ipHi));
103 }
104 else if( ipISO == ipHE_LIKE )
105 {
106 EnergyRyd = helike_energy( nelem, ipHi ) * WAVNRYD;
107 }
108 else
109 {
110 /* Other iso sequences don't exist yet. */
111 TotalInsanity();
112 }
113
114 /* >>chng 02 feb 09, change test to >= 0 since we now use 0 for 2s-2p */
115 ASSERT(EnergyRyd >= 0.);
116
117 iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd = EnergyRyd;
118 if (ipHi == 0)
119 EnergyRydGround = EnergyRyd;
120 iso_sp[ipISO][nelem].st[ipHi].energy().set(EnergyRydGround-EnergyRyd);
121
122 /* now loop from ground to level below ipHi */
123 for( ipLo=0; ipLo < ipHi; ipLo++ )
124 {
125 EnergyWN = RYD_INF * (iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd -
126 iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
127
128 /* This is the minimum line energy we will allow. */
129 /* \todo 2 wire this to lowest energy of code. */
130 if( EnergyWN==0 && ipISO==ipHE_LIKE )
131 EnergyWN = 0.0001;
132
133 if( EnergyWN < 0. )
134 EnergyWN = -1.0 * EnergyWN;
135
136 /* transition energy in various units: */
137 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() = (realnum)EnergyWN;
138
139 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() >= 0.);
140 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() >= 0.);
141 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyK() >= 0.);
142
144 if( N_(ipLo)==N_(ipHi) && ipISO==ipH_LIKE )
145 {
146 iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() = 0.;
147 }
148 else
149 {
150 /* make following an air wavelength */
151 iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() =
152 (realnum)(1.0e8/
153 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()/
154 RefIndex( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()));
155 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() > 0.);
156 }
157
158 }
159 }
160
161 /* fill the extra Lyman lines */
162 for( ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
163 {
164 FillExtraLymanLine( ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi], ipISO, nelem, ipHi );
165 }
166 }
167 }
168 }
169
170 /***************************************************************/
171 /***** Set up recombination tables for later interpolation *****/
172 /***************************************************************/
173 /* NB - the above is all we need if we are compiling recombination tables. */
174 iso_recomb_malloc();
175 iso_recomb_setup( ipH_LIKE );
176 iso_recomb_setup( ipHE_LIKE );
177 iso_recomb_auxiliary_free();
178
179 /* set up helium collision tables */
180 HeCollidSetup();
181
182 /***********************************************************************************/
183 /********** Transition Probabilities, Redistribution Functions, Opacitites ********/
184 /***********************************************************************************/
185 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
186 {
187 if( ipISO == ipH_LIKE )
188 {
189 /* do nothing here */
190 }
191 else if( ipISO == ipHE_LIKE )
192 {
193 /* This routine reads in transition probabilities from a file. */
194 HelikeTransProbSetup();
195 }
196 else
197 {
198 TotalInsanity();
199 }
200
201 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
202 {
203 /* must always do helium even if turned off */
204 if( nelem < 2 || dense.lgElmtOn[nelem] )
205 {
206 for( ipLo=ipH1s; ipLo < (iso_sp[ipISO][nelem].numLevels_max - 1); ipLo++ )
207 {
208 for( ipHi=ipLo + 1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
209 {
210 realnum Aul;
211
212 /* transition prob, EinstA uses current H atom indices */
213 if( ipISO == ipH_LIKE )
214 {
215 Aul = hydro_transprob( nelem, ipHi, ipLo );
216 }
217 else if( ipISO == ipHE_LIKE )
218 {
219 Aul = helike_transprob(nelem, ipHi, ipLo);
220 }
221 else
222 {
223 TotalInsanity();
224 }
225
226 if( Aul <= iso_ctrl.SmallA )
227 iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipEmis() = -1;
228 else
229 iso_sp[ipISO][nelem].trans(ipHi,ipLo).AddLine2Stack();
230
231 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = Aul;
232
233 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0.);
234
235 if( ipLo == 0 && ipHi == iso_ctrl.nLyaLevel[ipISO] )
236 {
237 long redis = iso_ctrl.ipLyaRedist[ipISO];
238 // H LyA has a special redistribution function
239 if( ipISO==ipH_LIKE && nelem==ipHYDROGEN )
240 redis = ipLY_A;
241 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = redis;
242 }
243 else if( ipLo == 0 )
244 {
245 /* these are rest of Lyman lines,
246 * complete redistribution, doppler core only, K2 core, default ipCRD */
247 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
248 }
249 else
250 {
251 /* all lines coming from excited states, default is complete
252 * redis with wings, ipCRDW*/
253 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = iso_ctrl.ipSubRedist[ipISO];
254 }
255
256 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA ||
257 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0.)
258 {
259 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() = 0.;
260 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() = 0.;
261 }
262 else
263 {
264 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() =
265 (realnum)(GetGF(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(),
266 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
267 iso_sp[ipISO][nelem].st[ipHi].g()));
268 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() > 0.);
269
270 /* derive the abs coef, call to function is gf, wl (A), g_low */
271 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() =
272 (realnum)(abscf(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf(),
273 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
274 iso_sp[ipISO][nelem].st[ipLo].g()));
275 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() > 0.);
276 }
277 }
278 }
279 }
280 }
281 }
282
283 /************************************************/
284 /********** Fine Structure Mixing - FSM ********/
285 /************************************************/
286 if( iso_ctrl.lgFSM[ipHE_LIKE] )
287 {
288 /* set some special optical depth values */
289 for( long nelem=ipHE_LIKE; nelem < LIMELM; nelem++ )
290 {
291 /* must always do helium even if turned off */
292 if( nelem < 2 || dense.lgElmtOn[nelem] )
293 {
294 for( ipHi=ipHe2s3S; ipHi<iso_sp[ipHE_LIKE][nelem].numLevels_max; ipHi++ )
295 {
296 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
297 {
298 DoFSMixing( nelem, ipLo, ipHi );
299 }
300 }
301 }
302 }
303 }
304
305 /* following comes out very slightly off, correct here */
306 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3s,ipH2s).WLAng() = 1640.f;
307 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3s,ipH2p).WLAng() = 1640.f;
308 if( iso_sp[ipH_LIKE][ipHELIUM].n_HighestResolved_max >=3 )
309 {
310 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3p,ipH2s).WLAng() = 1640.f;
311 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3p,ipH2p).WLAng() = 1640.f;
312 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3d,ipH2s).WLAng() = 1640.f;
313 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3d,ipH2p).WLAng() = 1640.f;
314 }
315
316 /****************************************************/
317 /********** lifetimes and damping constants ********/
318 /****************************************************/
319 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
320 {
321 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
322 {
323 /* define these for H and He always */
324 if( nelem < 2 || dense.lgElmtOn[nelem] )
325 {
326 /* these are not defined and must never be used */
327 iso_sp[ipISO][nelem].st[0].lifetime() = -FLT_MAX;
328
329 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
330 {
331 iso_sp[ipISO][nelem].st[ipHi].lifetime() = SMALLFLOAT;
332
333 for( ipLo=0; ipLo < ipHi; ipLo++ )
334 {
335 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
336 continue;
337
338 iso_sp[ipISO][nelem].st[ipHi].lifetime() += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
339 }
340
341 /* sum of A's was just stuffed, now invert for lifetime. */
342 iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./iso_sp[ipISO][nelem].st[ipHi].lifetime();
343
344 for( ipLo=0; ipLo < ipHi; ipLo++ )
345 {
346 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
347 continue;
348
349 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
350 continue;
351
352 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
353 (1.f/iso_sp[ipISO][nelem].st[ipHi].lifetime())/
354 PI4/iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
355
356 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
357 }
358 }
359 }
360 }
361 }
362
363 /* zero out some line information */
364 iso_zero();
365
366 /* loop over iso sequences */
367 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
368 {
369 for( long nelem = ipISO; nelem < LIMELM; nelem++ )
370 {
371 /* must always do helium even if turned off */
372 if( nelem == ipISO || dense.lgElmtOn[nelem] )
373 {
374 /* calculate cascade probabilities, branching ratios, and associated errors. */
375 iso_cascade( ipISO, nelem);
376 }
377 }
378 }
379
380 iso_satellite();
381
382 for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
383 iso_satellite_update( nelem );
384
385 /***************************************/
386 /********** Stark Broadening **********/
387 /***************************************/
388 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
389 {
390 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
391 {
392 if( nelem < 2 || dense.lgElmtOn[nelem] )
393 {
394 for( long ipHi= 1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
395 {
396 for( long ipLo=0; ipLo < ipHi; ++ipLo )
397 {
398 iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk = 0.;
399 iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up = 0.;
400 }
401 }
402 }
403 }
404 }
405
406 // arrays set up, do not allow number of levels to change in later sims
407 lgHydroMalloc = true;
408
409 return;
410}
411
412/* ============================================================================== */
413STATIC void iso_zero(void)
414{
415 DEBUG_ENTRY( "iso_zero()" );
416
417 hydro.HLineWidth = 0.;
418
419 /****************************************************/
420 /********** initialize some variables **********/
421 /****************************************************/
422 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
423 {
424 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
425 {
426 if( nelem < 2 || dense.lgElmtOn[nelem] )
427 {
428 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
429 {
430 iso_sp[ipISO][nelem].st[ipHi].Pop() = 0.;
431 iso_sp[ipISO][nelem].fb[ipHi].Reset();
432 }
433 if (ipISO <= nelem)
434 iso_sp[ipISO][nelem].st[0].Pop() =
435 dense.xIonDense[nelem][nelem-ipISO];
436 }
437
438 if( nelem < 2 )
439 {
440 iso_collapsed_bnl_set( ipISO, nelem );
441 //iso_collapsed_bnl_print( ipISO, nelem );
442 iso_collapsed_Aul_update( ipISO, nelem );
443 iso_collapsed_lifetimes_update( ipISO, nelem );
444 }
445 }
446 }
447
448 /* ground state of H and He is different since totally determine
449 * their own opacities */
450 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ConOpacRatio = 1e-5;
451 iso_sp[ipH_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
452 iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
453
454 return;
455}
456
457STATIC void iso_allocate(void)
458{
459
460 DEBUG_ENTRY( "iso_allocate()" );
461
462 /* the hydrogen and helium like iso-sequences */
463 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
464 {
465 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
466 {
467 /* only grab core for elements that are turned on */
468 if( nelem < 2 || dense.lgElmtOn[nelem] )
469 {
470 t_iso_sp *sp = &iso_sp[ipISO][nelem];
471 sp->numLevels_malloc = sp->numLevels_max;
472
473 ASSERT( sp->numLevels_max > 0 );
474 ASSERT( iso_ctrl.nLyman_malloc[ipISO] == iso_ctrl.nLyman[ipISO] );
475
476 sp->CachedAs.reserve( MAX2(1, sp->nCollapsed_max) );
477
478 sp->ipTrans.reserve( sp->numLevels_max );
479 sp->ex.reserve( sp->numLevels_max );
480 sp->CascadeProb.reserve( sp->numLevels_max );
481 sp->BranchRatio.reserve( sp->numLevels_max );
482 //sp->st.resize( sp->numLevels_max );
483 sp->fb.resize( sp->numLevels_max );
484
485 for( long i = 0; i < sp->nCollapsed_max; ++i )
486 {
487 sp->CachedAs.reserve( i, sp->numLevels_max - sp->nCollapsed_max );
488 for( long i1 = 0; i1 < sp->numLevels_max - sp->nCollapsed_max; ++i1 )
489 {
490 /* allocate two spaces delta L +/- 1 */
491 sp->CachedAs.reserve( i, i1, 2 );
492 }
493 }
494
495 sp->QuantumNumbers2Index.reserve( sp->n_HighestResolved_max + sp->nCollapsed_max + 1 );
496
497 for( long i = 1; i <= sp->n_HighestResolved_max + sp->nCollapsed_max; ++i )
498 {
499 /* Allocate proper number of angular momentum quantum number. */
500 sp->QuantumNumbers2Index.reserve( i, i );
501
502 for( long i1 = 0; i1 < i; ++i1 )
503 {
504 /* This may have to change for other iso sequences. */
505 ASSERT( NISO == 2 );
506 /* Allocate 4 spaces for multiplicity. H-like will be accessed with "2" for doublet,
507 * He-like will be accessed via "1" for singlet or "3" for triplet. "0" will not be used. */
508 sp->QuantumNumbers2Index.reserve( i, i1, 4 );
509 }
510 }
511
512 for( long n=1; n < sp->numLevels_max; ++n )
513 {
514 sp->ipTrans.reserve( n, n );
515 }
516
517 for( long n=0; n < sp->numLevels_max; ++n )
518 {
519 sp->ex.reserve( n, sp->numLevels_max );
520 sp->CascadeProb.reserve( n, sp->numLevels_max );
521 sp->BranchRatio.reserve( n, sp->numLevels_max );
522 }
523
524 sp->ipTrans.alloc();
525 sp->ex.alloc();
526 sp->CascadeProb.alloc();
527 sp->BranchRatio.alloc();
528
529 sp->CachedAs.alloc();
530 sp->QuantumNumbers2Index.alloc();
531 sp->bnl_effective.alloc( sp->QuantumNumbers2Index.clone() );
532 sp->QuantumNumbers2Index.invalidate();
533 sp->bnl_effective.invalidate();
534 }
535 }
536 }
537
538 ipSatelliteLines.reserve( NISO );
539 ipExtraLymanLines.reserve( NISO );
540
541 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
542 {
543 ipSatelliteLines.reserve( ipISO, LIMELM );
544 ipExtraLymanLines.reserve( ipISO, LIMELM );
545
546 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
547 {
548 /* only grab core for elements that are turned on */
549 if( nelem < 2 || dense.lgElmtOn[nelem] )
550 {
551 ASSERT( iso_sp[ipISO][nelem].numLevels_max > 0 );
552
553 ipSatelliteLines.reserve( ipISO, nelem, iso_sp[ipISO][nelem].numLevels_max );
554 ipExtraLymanLines.reserve( ipISO, nelem, iso_ctrl.nLyman_malloc[ipISO] );
555 }
556 }
557 }
558
559 ipSatelliteLines.alloc();
560 ipExtraLymanLines.alloc();
561
562 Transitions.resize(NISO);
563 SatelliteLines.resize(NISO);
564 ExtraLymanLines.resize(NISO);
565 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
566 {
567 Transitions[ipISO].reserve(LIMELM);
568 SatelliteLines[ipISO].reserve(LIMELM);
569 ExtraLymanLines[ipISO].reserve(LIMELM);
570 for( long nelem=0; nelem < ipISO; ++nelem )
571 {
572 Transitions[ipISO].push_back(
573 TransitionList("Insanity",&AnonStates));
574 SatelliteLines[ipISO].push_back(
575 TransitionList("Insanity",&AnonStates));
576 ExtraLymanLines[ipISO].push_back(
577 TransitionList("Insanity",&AnonStates));
578 }
579 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
580 {
581 if( nelem < 2 || dense.lgElmtOn[nelem] )
582 {
583 Transitions[ipISO].push_back(
584 TransitionList("Isosequence",&iso_sp[ipISO][nelem].st));
585 SatelliteLines[ipISO].push_back(
586 TransitionList("SatelliteLines",&iso_sp[ipISO][nelem].st));
587 ExtraLymanLines[ipISO].push_back(
588 TransitionList("ExtraLymanLines",&iso_sp[ipISO][nelem].st));
589 }
590 else
591 {
592 Transitions[ipISO].push_back(
593 TransitionList("Insanity",&AnonStates));
594 SatelliteLines[ipISO].push_back(
595 TransitionList("Insanity",&AnonStates));
596 ExtraLymanLines[ipISO].push_back(
597 TransitionList("Insanity",&AnonStates));
598 }
599 }
600 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
601 {
602 /* only grab core for elements that are turned on */
603 if( nelem < 2 || dense.lgElmtOn[nelem] )
604 {
605 if( iso_ctrl.lgDielRecom[ipISO] )
606 {
607 SatelliteLines[ipISO][nelem].resize( iso_sp[ipISO][nelem].numLevels_max );
608 AllTransitions.push_back(SatelliteLines[ipISO][nelem]);
609 unsigned int nLine = 0;
610 for( long ipLo=0; ipLo<iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
611 {
612 /* Upper level is continuum, use a generic state
613 * lower level is the same as the index. */
614 ipSatelliteLines[ipISO][nelem][ipLo] = nLine;
615 SatelliteLines[ipISO][nelem][nLine].Junk();
616 long ipHi = iso_sp[ipISO][nelem].numLevels_max;
617 SatelliteLines[ipISO][nelem][nLine].setHi(ipHi);
618 SatelliteLines[ipISO][nelem][nLine].setLo(ipLo);
619 SatelliteLines[ipISO][nelem][nLine].AddLine2Stack();
620 ++nLine;
621 }
622 ASSERT(SatelliteLines[ipISO][nelem].size() == nLine);
623 }
624
625 //iso_sp[ipISO][nelem].tr.resize( iso_sp[ipISO][nelem].ipTrans.size() );
626 //iso_sp[ipISO][nelem].tr.states() = &iso_sp[ipISO][nelem].st;
627 Transitions[ipISO][nelem].resize( iso_sp[ipISO][nelem].ipTrans.size() );
628 AllTransitions.push_back(Transitions[ipISO][nelem]);
629 unsigned int nTransition=0;
630 for( long ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
631 {
632 for( long ipLo=0; ipLo < ipHi; ipLo++ )
633 {
634 /* set ENTIRE array to impossible values, in case of bad pointer */
635 iso_sp[ipISO][nelem].ipTrans[ipHi][ipLo] = nTransition;
636 Transitions[ipISO][nelem][nTransition].Junk();
637 Transitions[ipISO][nelem][nTransition].setHi(ipHi);
638 Transitions[ipISO][nelem][nTransition].setLo(ipLo);
639 ++nTransition;
640 }
641 }
642 ASSERT(Transitions[ipISO][nelem].size() == nTransition);
643 iso_sp[ipISO][nelem].tr = &Transitions[ipISO][nelem];
644
645 /* junk the extra Lyman lines */
646 AllTransitions.push_back(ExtraLymanLines[ipISO][nelem]);
647 ExtraLymanLines[ipISO][nelem].resize(iso_ctrl.nLyman_malloc[ipISO]-2);
648 ExtraLymanLines[ipISO][nelem].states() = &iso_sp[ipISO][nelem].st;
649 unsigned int nExtraLyman = 0;
650 for( long ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
651 {
652 ipExtraLymanLines[ipISO][nelem][ipHi] = nExtraLyman;
653 ExtraLymanLines[ipISO][nelem][nExtraLyman].Junk();
654 long ipHi_offset = iso_sp[ipISO][nelem].numLevels_max + ipHi - 2;
655 if( iso_ctrl.lgDielRecom[ipISO] )
656 ipHi_offset += 1;
657 ExtraLymanLines[ipISO][nelem][nExtraLyman].setHi(ipHi_offset);
658 /* lower level is just ground state of the ion */
659 ExtraLymanLines[ipISO][nelem][nExtraLyman].setLo(0);
660 ExtraLymanLines[ipISO][nelem][nExtraLyman].AddLine2Stack();
661 ++nExtraLyman;
662 }
663 ASSERT(ExtraLymanLines[ipISO][nelem].size() == nExtraLyman);
664 }
665 }
666 }
667
668 // associate line and level stacks with species
669 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
670 {
671 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
672 {
673 if( dense.lgElmtOn[nelem] )
674 {
675 long ion = nelem - ipISO;
676 ASSERT( ion >= 0 && ion <= nelem );
677 char chLabel[6] = {'\0'}, chTemp[4] = {'\0'};
678 sprintf( chLabel, "%s", elementnames.chElementSym[nelem] );
679 if( chLabel[1]==' ' )
680 chLabel[1] = '\0';
681 else
682 chLabel[2] = '\0';
683 if( ion==1 )
684 sprintf( chTemp, "+" );
685 else if( ion>0 )
686 sprintf( chTemp, "+%li", ion );
687 strcat( chLabel, chTemp );
688
689 molecule *spmole = findspecies(chLabel);
690 ASSERT( spmole != null_mole );
691 mole.species[ spmole->index ].levels = &iso_sp[ipISO][nelem].st;
692 mole.species[ spmole->index ].lines = &Transitions[ipISO][nelem];
693 }
694 }
695 }
696
697 return;
698}
699
700STATIC void iso_assign_quantum_numbers(void)
701{
702 long int
703 ipLo,
704 level,
705 i,
706 in,
707 il,
708 is,
709 ij;
710
711 DEBUG_ENTRY( "iso_assign_quantum_numbers()" );
712
713 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
714 {
715 long ipISO = ipH_LIKE;
716 /* only check elements that are turned on */
717 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
718 {
719 i = 0;
720
721 /* 2 for doublet */
722 is = ipDOUBLET;
723
724 /* this loop is over quantum number n */
725 for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
726 {
727 for( il = 0L; il < in; ++il )
728 {
729 iso_sp[ipISO][nelem].st[i].n() = in;
730 iso_sp[ipISO][nelem].st[i].S() = is;
731 iso_sp[ipISO][nelem].st[i].l() = il;
732 iso_sp[ipISO][nelem].st[i].j() = -1;
733 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
734 ++i;
735 }
736 }
737 /* now do the collapsed levels */
738 in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
739 for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
740 {
741 iso_sp[ipISO][nelem].st[level].n() = in;
742 iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
743 iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
744 iso_sp[ipISO][nelem].st[level].j() = -1;
745 /* Point every l to same index for collapsed levels. */
746 for( il = 0; il < in; ++il )
747 {
748 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
749 }
750 ++in;
751 }
752 --in;
753
754 /* confirm that we did not overrun the array */
755 ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
756
757 /* confirm that n is positive and not greater than the max n. */
758 ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
759
760 /* Verify states and QuantumNumbers2Index agree in all cases */
761 for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
762 {
763 for( il = 0L; il < in; ++il )
764 {
765 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
766 if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
767 {
768 /* Must only check these for resolved levels...
769 * collapsed levels have pointers for l and s that will blow if used. */
770 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
771 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
772 }
773 }
774 }
775 }
776 }
777
778 /* then do he-like */
779 for( long nelem=ipHELIUM; nelem < LIMELM; nelem++ )
780 {
781 long ipISO = ipHE_LIKE;
782 /* only check elements that are turned on */
783 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
784 {
785 i = 0;
786
787 /* this loop is over quantum number n */
788 for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
789 {
790 for( il = 0L; il < in; ++il )
791 {
792 for( is = 3L; is >= 1L; is -= 2 )
793 {
794 /* All levels except singlet P follow the ordering scheme: */
795 /* lower l's have lower energy */
796 /* triplets have lower energy */
797 if( (il == 1L) && (is == 1L) )
798 continue;
799 /* n = 1 has no triplet, of course. */
800 if( (in == 1L) && (is == 3L) )
801 continue;
802
803 /* singlets */
804 if( is == 1 )
805 {
806 iso_sp[ipISO][nelem].st[i].n() = in;
807 iso_sp[ipISO][nelem].st[i].S() = is;
808 iso_sp[ipISO][nelem].st[i].l() = il;
809 /* this is not a typo, J=L for singlets. */
810 iso_sp[ipISO][nelem].st[i].j() = il;
811 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
812 ++i;
813 }
814 /* 2 triplet P is j-resolved */
815 else if( (in == 2) && (il == 1) && (is == 3) )
816 {
817 ij = 0;
818 do
819 {
820 iso_sp[ipISO][nelem].st[i].n() = in;
821 iso_sp[ipISO][nelem].st[i].S() = is;
822 iso_sp[ipISO][nelem].st[i].l() = il;
823 iso_sp[ipISO][nelem].st[i].j() = ij;
824 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
825 ++i;
826 ++ij;
827 /* repeat this for the separate j-levels within 2^3P. */
828 } while ( ij < 3 );
829 }
830 else
831 {
832 iso_sp[ipISO][nelem].st[i].n() = in;
833 iso_sp[ipISO][nelem].st[i].S() = is;
834 iso_sp[ipISO][nelem].st[i].l() = il;
835 iso_sp[ipISO][nelem].st[i].j() = -1L;
836 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
837 ++i;
838 }
839 }
840 }
841 /* Insert singlet P at the end of every sequence for a given n. */
842 if( in > 1L )
843 {
844 iso_sp[ipISO][nelem].st[i].n() = in;
845 iso_sp[ipISO][nelem].st[i].S() = 1L;
846 iso_sp[ipISO][nelem].st[i].l() = 1L;
847 iso_sp[ipISO][nelem].st[i].j() = 1L;
848 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][1][1] = i;
849 ++i;
850 }
851 }
852 /* now do the collapsed levels */
853 in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
854 for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
855 {
856 iso_sp[ipISO][nelem].st[level].n() = in;
857 iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
858 iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
859 iso_sp[ipISO][nelem].st[level].j() = -1;
860 /* Point every l and s to same index for collapsed levels. */
861 for( il = 0; il < in; ++il )
862 {
863 for( is = 1; is <= 3; is += 2 )
864 {
865 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
866 }
867 }
868 ++in;
869 }
870 --in;
871
872 /* confirm that we did not overrun the array */
873 ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
874
875 /* confirm that n is positive and not greater than the max n. */
876 ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
877
878 /* Verify states and QuantumNumbers2Index agree in all cases */
879 for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
880 {
881 for( il = 0L; il < in; ++il )
882 {
883 for( is = 3L; is >= 1; is -= 2 )
884 {
885 /* Ground state is not triplicate. */
886 if( (in == 1L) && (is == 3L) )
887 continue;
888
889 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
890 if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
891 {
892 /* Must only check these for resolved levels...
893 * collapsed levels have pointers for l and s that will blow if used. */
894 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
895 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
896 }
897 }
898 }
899 }
900 }
901 }
902
903 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
904 {
905 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
906 {
907 /* must always do helium even if turned off */
908 if( nelem < 2 || dense.lgElmtOn[nelem] )
909 {
910 for( ipLo=ipH1s; ipLo < iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
911 {
912 iso_sp[ipISO][nelem].st[ipLo].nelem() = (int)(nelem+1);
913 iso_sp[ipISO][nelem].st[ipLo].IonStg() = (int)(nelem+1-ipISO);
914
915 if( iso_sp[ipISO][nelem].st[ipLo].j() >= 0 )
916 {
917 iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*iso_sp[ipISO][nelem].st[ipLo].j()+1.f;
918 }
919 else if( iso_sp[ipISO][nelem].st[ipLo].l() >= 0 )
920 {
921 iso_sp[ipISO][nelem].st[ipLo].g() = (2.f*iso_sp[ipISO][nelem].st[ipLo].l()+1.f) *
922 iso_sp[ipISO][nelem].st[ipLo].S();
923 }
924 else
925 {
926 if( ipISO == ipH_LIKE )
927 iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
928 else if( ipISO == ipHE_LIKE )
929 iso_sp[ipISO][nelem].st[ipLo].g() = 4.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
930 else
931 {
932 /* replace this with correct thing if more sequences are added. */
933 TotalInsanity();
934 }
935 }
936 char chConfiguration[11] = " ";
937 long nCharactersWritten = 0;
938
939 ASSERT( iso_sp[ipISO][nelem].st[ipLo].n() < 1000 );
940
941 /* include j only if defined. */
942 if( iso_sp[ipISO][nelem].st[ipLo].n() > iso_sp[ipISO][nelem].n_HighestResolved_max )
943 {
944 nCharactersWritten = sprintf( chConfiguration, "n=%3li",
945 iso_sp[ipISO][nelem].st[ipLo].n() );
946 }
947 else if( iso_sp[ipISO][nelem].st[ipLo].j() > 0 )
948 {
949 nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c_%li",
950 iso_sp[ipISO][nelem].st[ipLo].n(),
951 iso_sp[ipISO][nelem].st[ipLo].S(),
952 chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l() ) ],
953 iso_sp[ipISO][nelem].st[ipLo].j() );
954 }
955 else
956 {
957 nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c",
958 iso_sp[ipISO][nelem].st[ipLo].n(),
959 iso_sp[ipISO][nelem].st[ipLo].S(),
960 chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l()) ] );
961 }
962
963 ASSERT( nCharactersWritten <= 10 );
964 chConfiguration[10] = '\0';
965
966 strncpy( iso_sp[ipISO][nelem].st[ipLo].chConfig(), chConfiguration, 10 );
967 }
968 }
969 }
970 }
971 return;
972}
973
974#if defined(__ICC) && defined(__i386)
975#pragma optimization_level 1
976#endif
977STATIC void FillExtraLymanLine( const TransitionList::iterator& t, long ipISO, long nelem, long nHi )
978{
979 double Enerwn, Aul;
980
981 DEBUG_ENTRY( "FillExtraLymanLine()" );
982
983 /* atomic number or charge and stage: */
984 (*(*t).Hi()).nelem() = (int)(nelem+1);
985 (*(*t).Hi()).IonStg() = (int)(nelem+1-ipISO);
986
987 (*(*t).Hi()).n() = nHi;
988
989 /* statistical weight is same as statistical weight of corresponding LyA. */
990 (*(*t).Hi()).g() = iso_sp[ipISO][nelem].st[iso_ctrl.nLyaLevel[ipISO]].g();
991
992 /* energies */
993 Enerwn = iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd * RYD_INF * ( 1. - 1./POW2((double)nHi) );
994
995 /* transition energy in various units:*/
996 (*t).EnergyWN() = (realnum)(Enerwn);
997 (*t).WLAng() = (realnum)(1.0e8/ Enerwn/ RefIndex(Enerwn));
998 (*(*t).Hi()).energy().set( Enerwn, "cm^-1" );
999
1000 if( ipISO == ipH_LIKE )
1001 {
1002 Aul = H_Einstein_A( nHi, 1, 1, 0, nelem+1 );
1003 }
1004 else
1005 {
1006 if( nelem == ipHELIUM )
1007 {
1008 /* A simple fit for the calculation of Helium lyman Aul's. */
1009 Aul = (1.508e10) / pow((double)nHi,2.975);
1010 }
1011 else
1012 {
1013 /* Fit to values given in
1014 * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
1015 * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
1016 /* originally astro.ph. 0201454 */
1017 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
1018 }
1019 }
1020
1021 (*t).Emis().Aul() = (realnum)Aul;
1022
1023 (*(*t).Hi()).lifetime() = iso_state_lifetime( ipISO, nelem, nHi, 1 );
1024
1025 (*t).Emis().dampXvel() = (realnum)( 1.f / (*(*t).Hi()).lifetime() / PI4 / (*t).EnergyWN() );
1026
1027 (*t).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
1028
1029 (*t).Emis().gf() = (realnum)(GetGF((*t).Emis().Aul(), (*t).EnergyWN(), (*(*t).Hi()).g()));
1030
1031 /* derive the abs coef, call to function is Emis().gf(), wl (A), g_low */
1032 (*t).Emis().opacity() = (realnum)(abscf((*t).Emis().gf(), (*t).EnergyWN(), (*(*t).Lo()).g()));
1033
1034 /* create array indices that will blow up */
1035 (*t).ipCont() = INT_MIN;
1036 (*t).Emis().ipFine() = INT_MIN;
1037
1038 {
1039 /* option to print particulars of some line when called
1040 * a prettier print statement is near where chSpin is defined below
1041 * search for "pretty print" */
1042 enum {DEBUG_LOC=false};
1043 if( DEBUG_LOC )
1044 {
1045 fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
1046 nelem+1,
1047 nHi,
1048 (*t).Emis().Aul() ,
1049 (*t).Emis().opacity()
1050 );
1051 }
1052 }
1053 return;
1054}
1055
1056/* calculate radiative lifetime of an individual iso state */
1057double iso_state_lifetime( long ipISO, long nelem, long n, long l )
1058{
1059 /* >>refer hydro lifetimes Horbatsch, M. W., Horbatsch, M. and Hessels, E. A. 2005, JPhysB, 38, 1765 */
1060
1061 double tau, t0, eps2;
1062 /* mass of electron */
1063 double m = ELECTRON_MASS;
1064 /* nuclear mass */
1065 double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT;
1066 double mu = (m*M)/(M+m);
1067 long z = 1;
1068 long Z = nelem + 1 - ipISO;
1069
1070 DEBUG_ENTRY( "iso_state_lifetime()" );
1071
1072 /* this should not be used for l=0 per the Horbatsch et al. paper */
1073 ASSERT( l > 0 );
1074
1075 eps2 = 1. - ( l*l + l + 8./47. - (l+1.)/69./n ) / POW2( (double)n );
1076
1077 t0 = 3. * H_BAR * pow( (double)n, 5.) /
1078 ( 2. * POW4( (double)( z * Z ) ) * pow( FINE_STRUCTURE, 5. ) * mu * POW2( SPEEDLIGHT ) ) *
1079 POW2( (m + M)/(Z*m + z*M) );
1080
1081 tau = t0 * ( 1. - eps2 ) /
1082 ( 1. + 19./88.*( (1./eps2 - 1.) * log( 1. - eps2 ) + 1. -
1083 0.5 * eps2 - 0.025 * eps2 * eps2 ) );
1084
1085 if( ipISO == ipHE_LIKE )
1086 {
1087 /* iso_state_lifetime is not spin specific, must exclude helike triplet here. */
1088 tau /= 3.;
1089 /* this is also necessary to correct the helike lifetimes */
1090 tau *= 1.1722 * pow( (double)nelem, 0.1 );
1091 }
1092
1093 /* would probably need a new lifetime algorithm for any other iso sequences. */
1094 ASSERT( ipISO <= ipHE_LIKE );
1095 ASSERT( tau > 0. );
1096
1097 return tau;
1098}
1099
1100/* calculate cascade probabilities, branching ratios, and associated errors. */
1101void iso_cascade( long ipISO, long nelem )
1102{
1103 /* The sum of all A's coming out of a given n,
1104 * Below we assert a monotonic trend. */
1105 double *SumAPerN;
1106
1107 long int i, j, ipLo, ipHi;
1108
1109 DEBUG_ENTRY( "iso_cascade()" );
1110
1111 SumAPerN = ((double*)MALLOC( (size_t)(iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double )));
1112 memset( SumAPerN, 0, (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double ) );
1113
1114 /* Initialize some ground state stuff, easier here than in loops. */
1115 iso_sp[ipISO][nelem].CascadeProb[0][0] = 1.;
1116 if( iso_ctrl.lgRandErrGen[ipISO] )
1117 {
1118 iso_sp[ipISO][nelem].fb[0].SigmaAtot = 0.;
1119 iso_sp[ipISO][nelem].ex[0][0].SigmaCascadeProb = 0.;
1120 }
1121
1122 /***************************************************************************/
1123 /****** Cascade probabilities, Branching ratios, and associated errors *****/
1124 /***************************************************************************/
1125 for( ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1126 {
1127 double SumAs = 0.;
1128
1134 /* initialize variables. */
1135 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipHi] = 1.;
1136 iso_sp[ipISO][nelem].CascadeProb[ipHi][0] = 0.;
1137 iso_sp[ipISO][nelem].BranchRatio[ipHi][0] = 0.;
1138
1139 if( iso_ctrl.lgRandErrGen[ipISO] )
1140 {
1141 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = 0.;
1142 iso_sp[ipISO][nelem].ex[ipHi][ipHi].SigmaCascadeProb = 0.;
1143 }
1144
1145 long ipLoStart = 0;
1146 if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) )
1147 ipLoStart = 1;
1148
1149 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1150 {
1151 SumAs += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1152 }
1153
1154 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1155 {
1156 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1157 {
1158 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1159 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] = 0.;
1160 continue;
1161 }
1162
1163 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1164 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] =
1165 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() / SumAs;
1166
1167 ASSERT( iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] <= 1.0000001 );
1168
1169 SumAPerN[N_(ipHi)] += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1170
1171 /* there are some negative energy transitions, where the order
1172 * has changed, but these are not optically allowed, these are
1173 * same n, different L, forbidden transitions */
1174 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() > 0. ||
1175 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA );
1176
1177 if( iso_ctrl.lgRandErrGen[ipISO] )
1178 {
1179 ASSERT( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] >= 0. );
1180 /* Uncertainties in A's are added in quadrature, square root is taken below. */
1181 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot +=
1182 pow( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] *
1183 (double)iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(), 2. );
1184 }
1185 }
1186
1187 if( iso_ctrl.lgRandErrGen[ipISO] )
1188 {
1189 /* Uncertainties in A's are added in quadrature above, square root taken here. */
1190 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = sqrt( iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot );
1191 }
1192
1193 /* cascade probabilities */
1194 for( i=0; i<ipHi; i++ )
1195 {
1196 for( ipLo=0; ipLo<=i; ipLo++ )
1197 {
1198 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] += iso_sp[ipISO][nelem].BranchRatio[ipHi][i] * iso_sp[ipISO][nelem].CascadeProb[i][ipLo];
1199 }
1200 }
1201
1202 if( iso_ctrl.lgRandErrGen[ipISO] )
1203 {
1204 for( ipLo=0; ipLo<ipHi; ipLo++ )
1205 {
1206 double SigmaCul = 0.;
1207 for( i=ipLo; i<ipHi; i++ )
1208 {
1209 if( iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul() > iso_ctrl.SmallA )
1210 {
1211 /* Uncertainties in A's and cascade probabilities */
1212 double SigmaA = iso_sp[ipISO][nelem].ex[ipHi][i].Error[IPRAD] *
1213 iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul();
1214 SigmaCul +=
1215 pow(SigmaA*iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
1216 pow(iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*iso_sp[ipISO][nelem].BranchRatio[ipHi][i]*
1217 iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
1218 pow(iso_sp[ipISO][nelem].ex[i][ipLo].SigmaCascadeProb*iso_sp[ipISO][nelem].BranchRatio[ipHi][i], 2.);
1219 }
1220 }
1221 SigmaCul = sqrt(SigmaCul);
1222 iso_sp[ipISO][nelem].ex[ipHi][ipLo].SigmaCascadeProb = SigmaCul;
1223 }
1224 }
1225 }
1226
1227 /************************************************************************/
1228 /*** Allowed decay conversion probabilities. See Robbins68b, Table 1. ***/
1229 /************************************************************************/
1230 {
1231 enum {DEBUG_LOC=false};
1232
1233 if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) )
1234 {
1235 /* To output Bm(n,l; ipLo), set ipLo, hi_l, and hi_s accordingly. */
1236 long int hi_l,hi_s;
1237 double Bm;
1238
1239 /* these must be set for following output to make sense
1240 * as is, a dangerous bit of code - set NaN for safety */
1241 hi_s = -100000;
1242 hi_l = -100000;
1243 ipLo = -100000;
1244 /* tripS to 2^3P */
1245 //hi_l = 0, hi_s = 3, ipLo = ipHe2p3P0;
1246
1247 /* tripD to 2^3P */
1248 //hi_l = 2, hi_s = 3, ipLo = ipHe2p3P0;
1249
1250 /* tripP to 2^3S */
1251 //hi_l = 1, hi_s = 3, ipLo = ipHe2s3S;
1252
1253 ASSERT( hi_l != iso_sp[ipISO][nelem].st[ipLo].l() );
1254
1255 fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
1256 fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");
1257
1258 for( ipHi=ipHe2p3P2; ipHi<iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipHi++ )
1259 {
1260 /* Pick out excitations from metastable 2tripS to ntripP. */
1261 if( (iso_sp[ipISO][nelem].st[ipHi].l() == 1) && (iso_sp[ipISO][nelem].st[ipHi].S() == 3) )
1262 {
1263 fprintf(ioQQQ,"\n%ld\t",iso_sp[ipISO][nelem].st[ipHi].n());
1264 j = 0;
1265 Bm = 0;
1266 for( i = ipLo; i<=ipHi; i++)
1267 {
1268 if( (iso_sp[ipISO][nelem].st[i].l() == hi_l) && (iso_sp[ipISO][nelem].st[i].S() == hi_s) )
1269 {
1270 if( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
1271 {
1272 Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * ( iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P0] +
1273 iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P1] + iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P2] );
1274 }
1275 else
1276 Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * iso_sp[ipISO][nelem].BranchRatio[i][ipLo];
1277
1278 if( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
1279 {
1280 j++;
1281 if(j == 3)
1282 {
1283 fprintf(ioQQQ,"%2.4e\t",Bm);
1284 Bm = 0;
1285 }
1286 }
1287 else
1288 {
1289 fprintf(ioQQQ,"%2.4e\t",Bm);
1290 Bm = 0;
1291 }
1292 }
1293 }
1294 }
1295 }
1296 fprintf(ioQQQ,"\n\n");
1297 }
1298 }
1299
1300 /******************************************************/
1301 /*** Lifetimes should increase monotonically with ***/
1302 /*** increasing n...Make sure the A's decrease. ***/
1303 /******************************************************/
1304 for( i=2; i < iso_sp[ipISO][nelem].n_HighestResolved_max; ++i)
1305 {
1306 ASSERT( (SumAPerN[i] > SumAPerN[i+1]) || opac.lgCaseB );
1307 }
1308
1309 {
1310 enum {DEBUG_LOC=false};
1311 if( DEBUG_LOC /* && (ipISO == ipH_LIKE) && (nelem == ipHYDROGEN) */)
1312 {
1313 for( i = 2; i<= (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max); ++i)
1314 {
1315 fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]);
1316 }
1317 }
1318 }
1319
1320 free( SumAPerN );
1321
1322 return;
1323}
1324
1326/* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002 */
1327/* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001. */
1328STATIC void iso_satellite( void )
1329{
1330 DEBUG_ENTRY( "iso_satellite()" );
1331
1332 for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
1333 {
1334 for( long nelem = ipISO; nelem < LIMELM; nelem++ )
1335 {
1336 if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
1337 {
1338 for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
1339 {
1340 char chLab[5]=" ";
1341
1342 TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
1343 (*tr).Zero();
1344
1345 /* Make approximation that all levels have energy of H-like 2s level */
1346 /* Lines to 1s2s have roughly energy of parent Ly-alpha. So lines to 1snL will have an energy
1347 * smaller by the difference between nL and 2s energies. Therefore, the following has
1348 * energy of parent Ly-alpha MINUS the difference between daughter level and daughter n=2 level. */
1349 (*tr).WLAng() = (realnum)(RYDLAM/
1350 ((iso_sp[ipISO-1][nelem].fb[0].xIsoLevNIonRyd - iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd) -
1351 (iso_sp[ipISO][nelem].fb[1].xIsoLevNIonRyd- iso_sp[ipISO][nelem].fb[i].xIsoLevNIonRyd)) );
1352
1353 (*tr).EnergyWN() = 1.e8f /
1354 (*tr).WLAng();
1355
1356 /* generate label for this ion */
1357 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
1358
1359 (*tr).Emis().iRedisFun() = ipCRDW;
1360 /* this is not the usual nelem, is it atomic not C scale. */
1361 (*(*tr).Hi()).nelem() = nelem + 1;
1362 (*(*tr).Hi()).IonStg() = nelem + 1 - ipISO;
1363 fixit(); /* what should the stat weight of the upper level be? For now say 2. */
1364 (*(*tr).Hi()).g() = 2.f;
1365 // The lower level must already be initialized.
1366 ASSERT( (*(*tr).Lo()).g() == iso_sp[ipISO][nelem].st[i].g() );
1367 //(*(*tr).Lo()).g() = iso_sp[ipISO][nelem].st[i].g();
1368 (*tr).Emis().PopOpc() =
1369 (*(*tr).Lo()).Pop();
1370
1371 (*tr).Emis().pump() = 0.;
1372
1373 }
1374 }
1375 }
1376 }
1377
1378 return;
1379}
1380
1381void iso_satellite_update( long nelem )
1382{
1383 double ConBoltz, LTE_pop=SMALLFLOAT+FLT_EPSILON, factor1, ConvLTEPOP;
1384
1385 DEBUG_ENTRY( "iso_satellite_update()" );
1386
1387 for( long ipISO = ipHE_LIKE; ipISO < MIN2(NISO,nelem+1); ipISO++ )
1388 {
1389 if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
1390 {
1391 for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
1392 {
1393 double dr_rate = iso_sp[ipISO][nelem].fb[i].DielecRecomb * iso_ctrl.lgDielRecom[ipISO];
1394
1395 TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
1396 (*tr).Emis().phots() =
1397 dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
1398
1399 (*tr).Emis().xIntensity() =
1400 (*tr).Emis().phots() *
1401 ERG1CM * (*tr).EnergyWN();
1402
1403 /* We set line intensity above using a rate, but here we need a transition probability.
1404 * We can obtain this by dividing dr_rate by the population of the autoionizing level.
1405 * We assume this level is in statistical equilibrium. */
1406 factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
1407 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
1408
1409 /* term in () is stat weight of electron * ion */
1410 ConvLTEPOP = pow(factor1,1.5)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;
1411
1412 /* This Boltzmann factor is exp( +ioniz energy / Te ). For simplicity, we make
1413 * the fair approximation that all of the autoionizing levels have an energy
1414 * equal to the parents n=2. */
1415 ConBoltz = dsexp(iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd/phycon.te_ryd);
1416
1417 if( ConBoltz >= SMALLDOUBLE )
1418 {
1419 /* The energy used to calculate ConBoltz above
1420 * should be negative since this is above the continuum, but
1421 * to be safe we calculate ConBoltz with a positive energy above
1422 * and multiply by it here instead of dividing. */
1423 LTE_pop = (*(*tr).Hi()).g() * ConBoltz * ConvLTEPOP;
1424 }
1425
1426 LTE_pop = max( LTE_pop, 1e-30f );
1427
1428 /* Now the transition probability is simply dr_rate/LTE_pop. */
1429 (*tr).Emis().Aul() = (realnum)(dr_rate/LTE_pop);
1430 (*tr).Emis().Aul() =
1431 max( iso_ctrl.SmallA, (*tr).Emis().Aul() );
1432
1433 (*tr).Emis().gf() = (realnum)GetGF(
1434 (*tr).Emis().Aul(),
1435 (*tr).EnergyWN(),
1436 (*(*tr).Hi()).g());
1437
1438 (*tr).Emis().gf() =
1439 max( 1e-20f, (*tr).Emis().gf() );
1440
1441 (*(*tr).Hi()).Pop() = LTE_pop * dense.xIonDense[nelem][nelem+1-ipISO] * dense.eden;
1442
1443 (*tr).Emis().PopOpc() =
1444 (*(*tr).Lo()).Pop() -
1445 (*(*tr).Hi()).Pop() *
1446 (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
1447
1448 (*tr).Emis().opacity() =
1449 (realnum)(abscf((*tr).Emis().gf(),
1450 (*tr).EnergyWN(),
1451 (*(*tr).Lo()).g()));
1452
1453 /* a typical transition probability is of order 1e10 s-1 */
1454 double lifetime = 1e-10;
1455
1456 (*tr).Emis().dampXvel() = (realnum)(
1457 (1.f/lifetime)/PI4/(*tr).EnergyWN());
1458 }
1459 }
1460 }
1461
1462 return;
1463}
1464
1465long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed )
1466{
1467 DEBUG_ENTRY( "iso_get_total_num_levels()" );
1468
1469 long tot_num_levels;
1470
1471 /* return the number of levels up to and including nmaxResolved PLUS
1472 * the number (numCollapsed) of collapsed n-levels */
1473
1474 if( ipISO == ipH_LIKE )
1475 {
1476 tot_num_levels = (long)( nmaxResolved * 0.5 *( nmaxResolved + 1 ) ) + numCollapsed;
1477 }
1478 else if( ipISO == ipHE_LIKE )
1479 {
1480 tot_num_levels = nmaxResolved*nmaxResolved + nmaxResolved + 1 + numCollapsed;
1481 }
1482 else
1483 TotalInsanity();
1484
1485 return tot_num_levels;
1486}
1487
1488void iso_update_num_levels( long ipISO, long nelem )
1489{
1490 DEBUG_ENTRY( "iso_update_num_levels()" );
1491
1492 /* This is the minimum resolved nmax. */
1493 ASSERT( iso_sp[ipISO][nelem].n_HighestResolved_max >= 3 );
1494
1495 iso_sp[ipISO][nelem].numLevels_max =
1496 iso_get_total_num_levels( ipISO, iso_sp[ipISO][nelem].n_HighestResolved_max, iso_sp[ipISO][nelem].nCollapsed_max );
1497
1498 if( iso_sp[ipISO][nelem].numLevels_max > iso_sp[ipISO][nelem].numLevels_malloc )
1499 {
1500 fprintf( ioQQQ, "The number of levels for ipISO %li, nelem %li, has been increased since the initial coreload.\n",
1501 ipISO, nelem );
1502 fprintf( ioQQQ, "This cannot be done.\n" );
1503 cdEXIT(EXIT_FAILURE);
1504 }
1505
1506 /* set local copies to the max values */
1507 iso_sp[ipISO][nelem].numLevels_local = iso_sp[ipISO][nelem].numLevels_max;
1508 iso_sp[ipISO][nelem].nCollapsed_local = iso_sp[ipISO][nelem].nCollapsed_max;
1509 iso_sp[ipISO][nelem].n_HighestResolved_local = iso_sp[ipISO][nelem].n_HighestResolved_max;
1510
1511 /* find the largest number of levels in any element in all iso sequences
1512 * we will allocate one matrix for ionization solver, and just use a piece of that memory
1513 * for smaller models. */
1514 max_num_levels = MAX2( max_num_levels, iso_sp[ipISO][nelem].numLevels_max);
1515
1516 return;
1517}
1518
1519void iso_collapsed_bnl_set( long ipISO, long nelem )
1520{
1521
1522 DEBUG_ENTRY( "iso_collapsed_bnl_set()" );
1523
1524 double bnl_array[4][3][4][10] = {
1525 {
1526 /* H */
1527 {
1528 {6.13E-01, 2.56E-01, 1.51E-01, 2.74E-01, 3.98E-01, 4.98E-01, 5.71E-01, 6.33E-01, 7.28E-01, 9.59E-01},
1529 {1.31E+00, 5.17E-01, 2.76E-01, 4.47E-01, 5.87E-01, 6.82E-01, 7.44E-01, 8.05E-01, 9.30E-01, 1.27E+00},
1530 {1.94E+00, 7.32E-01, 3.63E-01, 5.48E-01, 6.83E-01, 7.66E-01, 8.19E-01, 8.80E-01, 1.02E+00, 1.43E+00},
1531 {2.53E+00, 9.15E-01, 4.28E-01, 6.16E-01, 7.42E-01, 8.13E-01, 8.60E-01, 9.22E-01, 1.08E+00, 1.56E+00}
1532 },
1533 {
1534 {5.63E-01, 2.65E-01, 1.55E-01, 2.76E-01, 3.91E-01, 4.75E-01, 5.24E-01, 5.45E-01, 5.51E-01, 5.53E-01},
1535 {1.21E+00, 5.30E-01, 2.81E-01, 4.48E-01, 5.80E-01, 6.62E-01, 7.05E-01, 7.24E-01, 7.36E-01, 7.46E-01},
1536 {1.81E+00, 7.46E-01, 3.68E-01, 5.49E-01, 6.78E-01, 7.51E-01, 7.88E-01, 8.09E-01, 8.26E-01, 8.43E-01},
1537 {2.38E+00, 9.27E-01, 4.33E-01, 6.17E-01, 7.38E-01, 8.05E-01, 8.40E-01, 8.65E-01, 8.92E-01, 9.22E-01}
1538 },
1539 {
1540 {2.97E-01, 2.76E-01, 2.41E-01, 3.04E-01, 3.66E-01, 4.10E-01, 4.35E-01, 4.48E-01, 4.52E-01, 4.53E-01},
1541 {5.63E-01, 5.04E-01, 3.92E-01, 4.67E-01, 5.39E-01, 5.85E-01, 6.10E-01, 6.20E-01, 6.23E-01, 6.23E-01},
1542 {7.93E-01, 6.90E-01, 4.94E-01, 5.65E-01, 6.36E-01, 6.79E-01, 7.00E-01, 7.09E-01, 7.11E-01, 7.11E-01},
1543 {1.04E+00, 8.66E-01, 5.62E-01, 6.31E-01, 7.01E-01, 7.43E-01, 7.63E-01, 7.71E-01, 7.73E-01, 7.73E-01}
1544 }
1545 },
1546 {
1547 /* He+ */
1548 {
1549 {6.70E-02, 2.93E-02, 1.94E-02, 4.20E-02, 7.40E-02, 1.12E-01, 1.51E-01, 1.86E-01, 2.26E-01, 3.84E-01},
1550 {2.39E-01, 1.03E-01, 6.52E-02, 1.31E-01, 2.11E-01, 2.91E-01, 3.61E-01, 4.17E-01, 4.85E-01, 8.00E-01},
1551 {4.26E-01, 1.80E-01, 1.10E-01, 2.09E-01, 3.18E-01, 4.15E-01, 4.93E-01, 5.54E-01, 6.34E-01, 1.04E+00},
1552 {6.11E-01, 2.55E-01, 1.51E-01, 2.74E-01, 3.99E-01, 5.02E-01, 5.80E-01, 6.41E-01, 7.30E-01, 1.21E+00}
1553 },
1554 {
1555 {6.79E-02, 3.00E-02, 2.00E-02, 4.30E-02, 7.48E-02, 1.11E-01, 1.44E-01, 1.70E-01, 1.87E-01, 1.96E-01},
1556 {2.40E-01, 1.04E-01, 6.62E-02, 1.32E-01, 2.11E-01, 2.87E-01, 3.51E-01, 3.98E-01, 4.32E-01, 4.58E-01},
1557 {4.26E-01, 1.81E-01, 1.11E-01, 2.10E-01, 3.17E-01, 4.12E-01, 4.89E-01, 5.53E-01, 6.14E-01, 6.84E-01},
1558 {6.12E-01, 2.55E-01, 1.51E-01, 2.73E-01, 3.97E-01, 4.98E-01, 5.77E-01, 6.51E-01, 7.82E-01, 1.18E+00}
1559 },
1560 {
1561 {4.98E-02, 3.47E-02, 2.31E-02, 4.54E-02, 7.14E-02, 9.37E-02, 1.08E-01, 1.13E-01, 1.13E-01, 1.11E-01},
1562 {1.75E-01, 1.16E-01, 7.36E-02, 1.36E-01, 2.01E-01, 2.50E-01, 2.76E-01, 2.84E-01, 2.81E-01, 2.77E-01},
1563 {3.38E-01, 1.97E-01, 1.18E-01, 2.13E-01, 3.06E-01, 3.72E-01, 4.06E-01, 4.15E-01, 4.10E-01, 4.04E-01},
1564 {6.01E-01, 2.60E-01, 1.53E-01, 2.76E-01, 3.95E-01, 4.87E-01, 5.45E-01, 5.76E-01, 5.93E-01, 6.05E-01}
1565 }
1566 },
1567 {
1568 /* He singlets */
1569 {
1570 {1.77E-01, 3.59E-01, 1.54E-01, 2.75E-01, 3.98E-01, 4.94E-01, 5.51E-01, 5.68E-01, 5.46E-01, 4.97E-01},
1571 {4.09E-01, 7.23E-01, 2.83E-01, 4.48E-01, 5.89E-01, 6.78E-01, 7.22E-01, 7.30E-01, 7.07E-01, 6.65E-01},
1572 {6.40E-01, 1.02E+00, 3.74E-01, 5.49E-01, 6.85E-01, 7.63E-01, 7.98E-01, 8.03E-01, 7.84E-01, 7.53E-01},
1573 {8.70E-01, 1.28E+00, 4.42E-01, 6.17E-01, 7.44E-01, 8.13E-01, 8.42E-01, 8.46E-01, 8.34E-01, 8.13E-01}
1574 },
1575 {
1576 {1.78E-01, 3.62E-01, 1.55E-01, 2.73E-01, 3.91E-01, 4.73E-01, 5.10E-01, 5.04E-01, 4.70E-01, 4.32E-01},
1577 {4.08E-01, 7.26E-01, 2.83E-01, 4.45E-01, 5.79E-01, 6.54E-01, 6.78E-01, 6.64E-01, 6.30E-01, 5.98E-01},
1578 {6.37E-01, 1.03E+00, 3.73E-01, 5.46E-01, 6.75E-01, 7.40E-01, 7.57E-01, 7.43E-01, 7.15E-01, 6.92E-01},
1579 {8.65E-01, 1.28E+00, 4.41E-01, 6.14E-01, 7.35E-01, 7.92E-01, 8.05E-01, 7.95E-01, 7.74E-01, 7.59E-01}
1580 },
1581 {
1582 {2.07E-01, 3.73E-01, 1.73E-01, 2.85E-01, 4.03E-01, 4.76E-01, 5.06E-01, 5.03E-01, 4.84E-01, 4.63E-01},
1583 {4.32E-01, 7.13E-01, 3.06E-01, 4.54E-01, 5.81E-01, 6.44E-01, 6.59E-01, 6.49E-01, 6.28E-01, 6.11E-01},
1584 {6.40E-01, 9.85E-01, 3.98E-01, 5.53E-01, 6.74E-01, 7.27E-01, 7.36E-01, 7.26E-01, 7.10E-01, 6.98E-01},
1585 {8.38E-01, 1.21E+00, 4.67E-01, 6.20E-01, 7.34E-01, 7.79E-01, 7.87E-01, 7.79E-01, 7.69E-01, 7.63E-01}
1586 }
1587 },
1588 {
1589 /* He triplets */
1590 {
1591 {9.31E-02, 3.96E-01, 1.36E-01, 2.74E-01, 3.99E-01, 4.95E-01, 5.52E-01, 5.70E-01, 5.48E-01, 4.96E-01},
1592 {2.25E-01, 8.46E-01, 2.49E-01, 4.46E-01, 5.89E-01, 6.79E-01, 7.23E-01, 7.31E-01, 7.08E-01, 6.64E-01},
1593 {3.59E-01, 1.24E+00, 3.30E-01, 5.47E-01, 6.85E-01, 7.63E-01, 7.98E-01, 8.04E-01, 7.85E-01, 7.53E-01},
1594 {4.93E-01, 1.60E+00, 3.91E-01, 6.15E-01, 7.44E-01, 8.13E-01, 8.42E-01, 8.47E-01, 8.35E-01, 8.12E-01}
1595 },
1596 {
1597 {9.32E-02, 3.99E-01, 1.35E-01, 2.72E-01, 3.91E-01, 4.75E-01, 5.14E-01, 5.09E-01, 4.73E-01, 4.31E-01},
1598 {2.25E-01, 8.49E-01, 2.49E-01, 4.44E-01, 5.79E-01, 6.56E-01, 6.81E-01, 6.68E-01, 6.31E-01, 5.96E-01},
1599 {3.58E-01, 1.25E+00, 3.29E-01, 5.44E-01, 6.76E-01, 7.42E-01, 7.60E-01, 7.46E-01, 7.16E-01, 6.91E-01},
1600 {4.92E-01, 1.60E+00, 3.90E-01, 6.12E-01, 7.36E-01, 7.93E-01, 8.07E-01, 7.97E-01, 7.74E-01, 7.58E-01}
1601 },
1602 {
1603 {1.13E-01, 4.21E-01, 1.47E-01, 2.83E-01, 4.04E-01, 4.80E-01, 5.13E-01, 5.12E-01, 4.93E-01, 4.71E-01},
1604 {2.52E-01, 8.56E-01, 2.61E-01, 4.50E-01, 5.82E-01, 6.48E-01, 6.66E-01, 6.56E-01, 6.35E-01, 6.16E-01},
1605 {3.85E-01, 1.23E+00, 3.41E-01, 5.49E-01, 6.75E-01, 7.30E-01, 7.41E-01, 7.31E-01, 7.15E-01, 7.02E-01},
1606 {5.14E-01, 1.56E+00, 4.01E-01, 6.15E-01, 7.34E-01, 7.82E-01, 7.90E-01, 7.83E-01, 7.72E-01, 7.65E-01}
1607 }
1608 }
1609 };
1610
1611 double temps[4] = {5000., 10000., 15000., 20000. };
1612 double log_dens[3] = {2., 4., 6.};
1613 long ipTe, ipDens;
1614
1615 ASSERT( nelem <= 1 );
1616
1617 /* find temperature in tabulated values. */
1618 ipTe = hunt_bisect( temps, 4, phycon.te );
1619 ipDens = hunt_bisect( log_dens, 3, log10(dense.eden) );
1620
1621 ASSERT( (ipTe >=0) && (ipTe < 3) );
1622 ASSERT( (ipDens >=0) && (ipDens < 2) );
1623
1624 for( long nHi=iso_sp[ipISO][nelem].n_HighestResolved_max+1; nHi<=iso_sp[ipISO][nelem].n_HighestResolved_max+iso_sp[ipISO][nelem].nCollapsed_max; nHi++ )
1625 {
1626 for( long lHi=0; lHi<nHi; lHi++ )
1627 {
1628 for( long sHi=1; sHi<4; sHi++ )
1629 {
1630 if( ipISO == ipH_LIKE && sHi != 2 )
1631 continue;
1632 else if( ipISO == ipHE_LIKE && sHi != 1 && sHi != 3 )
1633 continue;
1634
1635 double bnl_at_lo_den, bnl_at_hi_den, bnl;
1636 double bnl_max, bnl_min, temp, dens;
1637
1638 long ipL = MIN2(9,lHi);
1639 long ip1;
1640
1641 if( nelem==ipHYDROGEN )
1642 ip1 = 0;
1643 else if( nelem==ipHELIUM )
1644 {
1645 if( ipISO==ipH_LIKE )
1646 ip1 = 1;
1647 else if( ipISO==ipHE_LIKE )
1648 {
1649 if( sHi==1 )
1650 ip1 = 2;
1651 else if( sHi==3 )
1652 ip1 = 3;
1653 else
1654 TotalInsanity();
1655 }
1656 else
1657 TotalInsanity();
1658 }
1659 else
1660 TotalInsanity();
1661
1662 temp = MAX2( temps[0], phycon.te );
1663 temp = MIN2( temps[3], temp );
1664
1665 dens = MAX2( log_dens[0], log10(dense.eden) );
1666 dens = MIN2( log_dens[2], dens );
1667
1668 /* Calculate the answer...must interpolate on two variables.
1669 * First interpolate on T, at both the lower and upper densities.
1670 * Then interpolate between these results for the right density. */
1671
1672 if( temp < temps[0] && dens < log_dens[0] )
1673 bnl = bnl_array[ip1][0][0][ipL];
1674 else if( temp < temps[0] && dens >= log_dens[2] )
1675 bnl = bnl_array[ip1][2][0][ipL];
1676 else if( temp >= temps[3] && dens < log_dens[0] )
1677 bnl = bnl_array[ip1][0][3][ipL];
1678 else if( temp >= temps[3] && dens >= log_dens[2] )
1679 bnl = bnl_array[ip1][2][3][ipL];
1680 else
1681 {
1682 bnl_at_lo_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
1683 (bnl_array[ip1][ipDens][ipTe+1][ipL] - bnl_array[ip1][ipDens][ipTe][ipL]) + bnl_array[ip1][ipDens][ipTe][ipL];
1684
1685 bnl_at_hi_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
1686 (bnl_array[ip1][ipDens+1][ipTe+1][ipL] - bnl_array[ip1][ipDens+1][ipTe][ipL]) + bnl_array[ip1][ipDens+1][ipTe][ipL];
1687
1688 bnl = ( dens - log_dens[ipDens]) / (log_dens[ipDens+1] - log_dens[ipDens]) *
1689 (bnl_at_hi_den - bnl_at_lo_den) + bnl_at_lo_den;
1690 }
1691
1693 {
1694 bnl_max = MAX4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
1695 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
1696 ASSERT( bnl <= bnl_max );
1697
1698 bnl_min = MIN4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
1699 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
1700 ASSERT( bnl >= bnl_min );
1701 }
1702
1703 iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] = bnl;
1704
1705 ASSERT( iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] > 0. );
1706 }
1707 }
1708 }
1709
1710 return;
1711}
1712
1713
1714void iso_collapsed_bnl_print( long ipISO, long nelem )
1715{
1716 DEBUG_ENTRY( "iso_collapsed_bnl_print()" );
1717
1718 for( long is = 1; is<=3; ++is)
1719 {
1720 if( ipISO == ipH_LIKE && is != 2 )
1721 continue;
1722 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
1723 continue;
1724
1725 char chSpin[3][9]= {"singlets", "doublets", "triplets"};
1726
1727 /* give element number and spin */
1728 fprintf(ioQQQ," %s %s %s bnl\n",
1729 iso_ctrl.chISO[ipISO],
1730 elementnames.chElementSym[nelem],
1731 chSpin[is-1]);
1732
1733 /* header with the l states */
1734 fprintf(ioQQQ," n\\l=> ");
1735 for( long i =0; i < iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++i)
1736 {
1737 fprintf(ioQQQ,"%2ld ",i);
1738 }
1739 fprintf(ioQQQ,"\n");
1740
1741 /* loop over prin quant numbers, one per line, with l across */
1742 for( long in = 1; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in)
1743 {
1744 if( is==3 && in==1 )
1745 continue;
1746
1747 fprintf(ioQQQ," %2ld ",in);
1748
1749 for( long il = 0; il < in; ++il)
1750 {
1751 fprintf( ioQQQ, "%9.3e ", iso_sp[ipISO][nelem].bnl_effective[in][il][is] );
1752 }
1753 fprintf(ioQQQ,"\n");
1754 }
1755 }
1756
1757 return;
1758}
1759
1760void iso_collapsed_Aul_update( long ipISO, long nelem )
1761{
1762 DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
1763
1764 long ipFirstCollapsed = iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max;
1765
1766 for( long ipLo=0; ipLo<ipFirstCollapsed; ipLo++ )
1767 {
1768 long spin = iso_sp[ipISO][nelem].st[ipLo].S();
1769
1770 /* calculate effective Aul's from collapsed levels */
1771 for( long nHi=iso_sp[ipISO][nelem].n_HighestResolved_max+1; nHi<=iso_sp[ipISO][nelem].n_HighestResolved_max+iso_sp[ipISO][nelem].nCollapsed_max; nHi++ )
1772 {
1773 realnum Auls[2] = {
1774 iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][0],
1775 iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] };
1776
1777 realnum EffectiveAul =
1778 Auls[0]*spin*(2.f*(L_(ipLo)+1.f)+1.f)*(realnum)iso_sp[ipISO][nelem].bnl_effective[nHi][ L_(ipLo)+1 ][spin];
1779
1780 /* this is for n,L-1 -> n',L
1781 * make sure L-1 exists. */
1782 if( L_(ipLo) > 0 )
1783 {
1784 EffectiveAul +=
1785 Auls[1]*spin*(2.f*(L_(ipLo)-1.f)+1.f)*(realnum)iso_sp[ipISO][nelem].bnl_effective[nHi][ L_(ipLo)-1 ][spin];
1786 }
1787
1788 if( ipISO==ipH_LIKE )
1789 EffectiveAul /= (2.f*nHi*nHi);
1790 else if( ipISO==ipHE_LIKE )
1791 EffectiveAul /= (4.f*nHi*nHi);
1792 else
1793 TotalInsanity();
1794
1795 long ipHi = iso_sp[ipISO][nelem].QuantumNumbers2Index[nHi][ L_(ipLo)+1 ][spin];
1796
1797 /* FINALLY, put the effective A in the proper Emis structure. */
1798 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = EffectiveAul;
1799
1800 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0. );
1801 }
1802 }
1803
1804 return;
1805}
1806
1807void iso_collapsed_lifetimes_update( long ipISO, long nelem )
1808{
1809 DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
1810
1811 for( long ipHi=iso_sp[ipISO][nelem].numLevels_max- iso_sp[ipISO][nelem].nCollapsed_max; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1812 {
1813 iso_sp[ipISO][nelem].st[ipHi].lifetime() = SMALLFLOAT;
1814
1815 for( long ipLo=0; ipLo < ipHi; ipLo++ )
1816 {
1817 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1818 continue;
1819
1820 iso_sp[ipISO][nelem].st[ipHi].lifetime() += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1821 }
1822
1823 /* sum of A's was just stuffed, now invert for lifetime. */
1824 iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./iso_sp[ipISO][nelem].st[ipHi].lifetime();
1825
1826 for( long ipLo=0; ipLo < ipHi; ipLo++ )
1827 {
1828 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
1829 continue;
1830
1831 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1832 continue;
1833
1834 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
1835 (1.f/iso_sp[ipISO][nelem].st[ipHi].lifetime())/
1836 PI4/iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
1837
1838 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
1839 }
1840 }
1841
1842 return;
1843}
1844
1845#if 0
1846STATIC void Prt_AGN_table( void )
1847{
1848 /* the designation of the levels, chLevel[n][string] */
1849 multi_arr<char,2> chLevel(max_num_levels,10);
1850
1851 /* create spectroscopic designation of labels */
1852 for( long ipLo=0; ipLo < iso_sp[ipISO][ipISO].numLevels_max-iso_sp[ipISO][ipISO].nCollapsed_max; ++ipLo )
1853 {
1854 long nelem = ipISO;
1855 sprintf( &chLevel.front(ipLo) , "%li %li%c", N_(ipLo), S_(ipLo), chL[MIN2(20,L_(ipLo))] );
1856 }
1857
1858 /* option to print cs data for AGN */
1859 /* create spectroscopic designation of labels */
1860 {
1861 /* option to print particulars of some line when called */
1862 enum {AGN=false};
1863 if( AGN )
1864 {
1865# define NTEMP 6
1866 double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
1867 double telog[NTEMP] ,
1868 cs ,
1869 ratecoef;
1870 long nelem = ipHELIUM;
1871 fprintf( ioQQQ,"trans");
1872 for( long i=0; i < NTEMP; ++i )
1873 {
1874 telog[i] = log10( te[i] );
1875 fprintf( ioQQQ,"\t%.3e",te[i]);
1876 }
1877 for( long i=0; i < NTEMP; ++i )
1878 {
1879 fprintf( ioQQQ,"\t%.3e",te[i]);
1880 }
1881 fprintf(ioQQQ,"\n");
1882
1883 for( long ipHi=ipHe2s3S; ipHi< iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max; ++ipHi )
1884 {
1885 for( long ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
1886 {
1887
1888 /* deltaN = 0 transitions may be wrong because
1889 * COLL_CONST below is only correct for electron colliders */
1890 if( N_(ipHi) == N_(ipLo) )
1891 continue;
1892
1893 /* print the designations of the lower and upper levels */
1894 fprintf( ioQQQ,"%s - %s",
1895 &chLevel.front(ipLo) , &chLevel.front(ipHi) );
1896
1897 /* print the interpolated collision strengths */
1898 for( long i=0; i < NTEMP; ++i )
1899 {
1900 phycon.alogte = telog[i];
1901 /* print cs */
1902 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
1903 fprintf(ioQQQ,"\t%.2e", cs );
1904 }
1905
1906 /* print the rate coefficients */
1907 for( long i=0; i < NTEMP; ++i )
1908 {
1909 phycon.alogte = telog[i];
1910 phycon.te = pow(10.,telog[i] );
1911 tfidle(false);
1912 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
1913 /* collisional deexcitation rate */
1914 ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso_sp[ipHE_LIKE][nelem].st[ipLo].g() *
1915 sexp( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyK() / phycon.te );
1916 fprintf(ioQQQ,"\t%.2e", ratecoef );
1917 }
1918 fprintf(ioQQQ,"\n");
1919 }
1920 }
1921 cdEXIT(EXIT_FAILURE);
1922 }
1923 }
1924
1925 return;
1926}
1927#endif
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
STATIC
#define STATIC
Definition: cddefines.h:97
MIN4
#define MIN4(a, b, c, d)
Definition: cddefines.h:771
MIN2
#define MIN2
Definition: cddefines.h:761
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
LIMELM
const int LIMELM
Definition: cddefines.h:258
POW2
#define POW2
Definition: cddefines.h:929
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
NISO
const int NISO
Definition: cddefines.h:261
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
realnum
float realnum
Definition: cddefines.h:103
MAX2
#define MAX2
Definition: cddefines.h:782
POW4
#define POW4
Definition: cddefines.h:943
max
long max(int a, long b)
Definition: cddefines.h:775
ipCRDW
const int ipCRDW
Definition: cddefines.h:294
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
MAX4
#define MAX4(a, b, c, d)
Definition: cddefines.h:792
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
dsexp
double dsexp(double x)
Definition: service.cpp:953
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipLY_A
const int ipLY_A
Definition: cddefines.h:296
fixit
void fixit(void)
Definition: service.cpp:991
lgHydroMalloc
bool lgHydroMalloc
Definition: cdinit.cpp:61
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:403
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
EmissionProxy::opacity
realnum & opacity() const
Definition: emission.h:593
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:175
TransitionProxy::ipEmis
int & ipEmis() const
Definition: transition.h:416
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:664
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
molecule
Definition: mole.h:132
molecule::index
int index
Definition: mole.h:169
multi_arr
Definition: container_classes.h:942
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1080
multi_arr::invalidate
void invalidate()
Definition: container_classes.h:1057
multi_arr::clone
const multi_geom< d, ALLOC > & clone() const
Definition: container_classes.h:1784
multi_arr::alloc
void alloc()
Definition: container_classes.h:1116
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat.h:329
t_iso_sp
Definition: iso.h:442
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:567
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
t_iso_sp::ex
multi_arr< extra_tr, 2 > ex
Definition: iso.h:449
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
t_iso_sp::BranchRatio
multi_arr< double, 2 > BranchRatio
Definition: iso.h:451
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
t_iso_sp::numLevels_malloc
long int numLevels_malloc
Definition: iso.h:502
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_iso_sp::n_HighestResolved_local
long int n_HighestResolved_local
Definition: iso.h:507
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:487
t_iso_sp::CascadeProb
multi_arr< double, 2 > CascadeProb
Definition: iso.h:450
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_iso_sp::st
qList st
Definition: iso.h:453
t_iso_sp::ipTrans
multi_arr< long, 2 > ipTrans
Definition: iso.h:448
t_iso_sp::bnl_effective
multi_arr< double, 3 > bnl_effective
Definition: iso.h:566
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
t_iso_sp::tr
TransitionList * tr
Definition: iso.h:454
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:365
t_isoCTRL::lgFSM
bool lgFSM[NISO]
Definition: iso.h:399
t_isoCTRL::ipLyaRedist
int ipLyaRedist[NISO]
Definition: iso.h:374
t_isoCTRL::ipSubRedist
int ipSubRedist[NISO]
Definition: iso.h:374
t_isoCTRL::nLyaLevel
int nLyaLevel[NISO]
Definition: iso.h:377
t_isoCTRL::nLyman_malloc
long int nLyman_malloc[NISO]
Definition: iso.h:336
t_isoCTRL::ipResoRedist
int ipResoRedist[NISO]
Definition: iso.h:374
t_isoCTRL::stat_ion
realnum stat_ion[NISO]
Definition: iso.h:362
t_isoCTRL::chISO
const char * chISO[NISO]
Definition: iso.h:330
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:334
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
t_isoCTRL::lgRandErrGen
bool lgRandErrGen[NISO]
Definition: iso.h:403
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
ipELECTRON
@ ipELECTRON
Definition: collision.h:9
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:195
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
dense
t_dense dense
Definition: dense.cpp:24
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
HeCollidSetup
void HeCollidSetup(void)
Definition: helike_cs.cpp:90
helike_energy
double helike_energy(long nelem, long ipLev)
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1059
HeCSInterp
realnum HeCSInterp(long int nelem, long int ipHi, long int ipLo, long int Collider)
Definition: helike_cs.cpp:227
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:894
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1132
helike_einsta.h
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
hydro_transprob
realnum hydro_transprob(long nelem, long ipHi, long ipLo)
Definition: hydroeinsta.cpp:45
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
max_num_levels
long int max_num_levels
Definition: iso.cpp:10
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
ipH3p
const int ipH3p
Definition: iso.h:31
IPRAD
#define IPRAD
Definition: iso.h:86
ipH1s
const int ipH1s
Definition: iso.h:27
iso_recomb_malloc
void iso_recomb_malloc(void)
Definition: iso_radiative_recomb.cpp:765
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
N_
#define N_(A_)
Definition: iso.h:20
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
ipH3s
const int ipH3s
Definition: iso.h:30
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
iso_recomb_setup
void iso_recomb_setup(long ipISO)
Definition: iso_radiative_recomb.cpp:838
S_
#define S_(A_)
Definition: iso.h:22
ipDOUBLET
@ ipDOUBLET
Definition: iso.h:82
ipH2p
const int ipH2p
Definition: iso.h:29
ipH3d
const int ipH3d
Definition: iso.h:32
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
ipH2s
const int ipH2s
Definition: iso.h:28
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
iso_recomb_auxiliary_free
void iso_recomb_auxiliary_free(void)
Definition: iso_radiative_recomb.cpp:825
L_
#define L_(A_)
Definition: iso.h:21
iso_update_num_levels
void iso_update_num_levels(long ipISO, long nelem)
Definition: iso_create.cpp:1488
FillExtraLymanLine
STATIC void FillExtraLymanLine(const TransitionList::iterator &t, long ipISO, long nelem, long nHi)
Definition: iso_create.cpp:977
iso_collapsed_bnl_print
void iso_collapsed_bnl_print(long ipISO, long nelem)
Definition: iso_create.cpp:1714
iso_zero
STATIC void iso_zero(void)
Definition: iso_create.cpp:413
iso_allocate
STATIC void iso_allocate(void)
Definition: iso_create.cpp:457
iso_collapsed_lifetimes_update
void iso_collapsed_lifetimes_update(long ipISO, long nelem)
Definition: iso_create.cpp:1807
iso_satellite_update
void iso_satellite_update(long nelem)
Definition: iso_create.cpp:1381
iso_cascade
void iso_cascade(long ipISO, long nelem)
Definition: iso_create.cpp:1101
iso_assign_quantum_numbers
STATIC void iso_assign_quantum_numbers(void)
Definition: iso_create.cpp:700
iso_get_total_num_levels
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
Definition: iso_create.cpp:1465
iso_collapsed_Aul_update
void iso_collapsed_Aul_update(long ipISO, long nelem)
Definition: iso_create.cpp:1760
iso_collapsed_bnl_set
void iso_collapsed_bnl_set(long ipISO, long nelem)
Definition: iso_create.cpp:1519
iso_state_lifetime
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
Definition: iso_create.cpp:1057
chL
char chL[21]
Definition: iso_create.cpp:37
iso_create
void iso_create(void)
Definition: iso_create.cpp:39
iso_satellite
STATIC void iso_satellite(void)
Definition: iso_create.cpp:1328
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:141
abscf
double abscf(double gf, double enercm, double gl)
Definition: lines_service.cpp:122
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:101
lines_service.h
mole
t_mole_local mole
Definition: mole.cpp:7
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
null_mole
molecule * null_mole
Definition: mole_species.cpp:65
opac
t_opac opac
Definition: opacity.cpp:5
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
phycon
t_phycon phycon
Definition: phycon.cpp:6
H_BAR
UNUSED const double H_BAR
Definition: physconst.h:144
FINE_STRUCTURE
UNUSED const double FINE_STRUCTURE
Definition: physconst.h:216
PI4
UNUSED const double PI4
Definition: physconst.h:35
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
EVRYD
UNUSED const double EVRYD
Definition: physconst.h:189
ERG1CM
UNUSED const double ERG1CM
Definition: physconst.h:164
ATOMIC_MASS_UNIT
UNUSED const double ATOMIC_MASS_UNIT
Definition: physconst.h:88
COLL_CONST
UNUSED const double COLL_CONST
Definition: physconst.h:229
WAVNRYD
UNUSED const double WAVNRYD
Definition: physconst.h:173
HION_LTE_POP
UNUSED const double HION_LTE_POP
Definition: physconst.h:157
RYDLAM
UNUSED const double RYDLAM
Definition: physconst.h:176
nLine
static long int nLine
Definition: save_line.cpp:288
ex
static double * ex
Definition: species2.cpp:28
g
static double * g
Definition: species2.cpp:28
t_dense::eden
double eden
Definition: dense.h:190
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_hydro::HLineWidth
realnum HLineWidth
Definition: hydrogenic.h:63
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:161
t_phycon::te
double te
Definition: phycon.h:11
t_phycon::te32
double te32
Definition: phycon.h:49
t_phycon::te_ryd
double te_ryd
Definition: phycon.h:17
t_phycon::alogte
double alogte
Definition: phycon.h:82
SatelliteLines
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:38
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
AnonStates
qList AnonStates(1)
ipSatelliteLines
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:37
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
Transitions
vector< vector< TransitionList > > Transitions
Definition: taulines.cpp:33
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
tfidle
STATIC void tfidle(bool lgForceUpdate)
Definition: temp_change.cpp:150
M
static const int M
Definition: thirdparty.cpp:2815
hunt_bisect
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:270
Generated by doxygen 1.9.5