cloudy  trunk
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 */
25 STATIC void iso_zero(void);
26 
27 /* allocate memory for iso sequence structures */
28 STATIC void iso_allocate(void);
29 
30 /* define levels of iso sequences and assign quantum numbers to those levels */
31 STATIC void iso_assign_quantum_numbers(void);
32 
33 STATIC void FillExtraLymanLine( const TransitionList::iterator& t, long ipISO, long nelem, long nHi );
34 
35 STATIC void iso_satellite( void );
36 
37 char chL[21]={'S','P','D','F','G','H','I','K','L','M','N','O','Q','R','T','U','V','W','X','Y','Z'};
38 
39 void 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 /* ============================================================================== */
413 STATIC 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 
457 STATIC 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 
700 STATIC 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
977 STATIC 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 */
1057 double 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. */
1101 void 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. */
1328 STATIC 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 
1381 void 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 
1465 long 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 
1488 void 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 
1519 void 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 
1714 void 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 
1760 void 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 
1807 void 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
1846 STATIC 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
ipDOUBLET
@ ipDOUBLET
Definition: iso.h:82
t_isoCTRL::nLyaLevel
int nLyaLevel[NISO]
Definition: iso.h:377
multi_arr::invalidate
void invalidate()
Definition: container_classes.h:1057
helike.h
Transitions
vector< vector< TransitionList > > Transitions
Definition: taulines.cpp:33
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
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:195
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
AnonStates
qList AnonStates(1)
t_dense::eden
double eden
Definition: dense.h:190
dense
t_dense dense
Definition: dense.cpp:24
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:161
elementnames.h
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
secondaries.h
iso_collapsed_lifetimes_update
void iso_collapsed_lifetimes_update(long ipISO, long nelem)
Definition: iso_create.cpp:1807
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:175
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
multi_arr::clone
const multi_geom< d, ALLOC > & clone() const
Definition: container_classes.h:1784
t_isoCTRL::chISO
const char * chISO[NISO]
Definition: iso.h:330
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1059
realnum
float realnum
Definition: cddefines.h:103
tfidle
STATIC void tfidle(bool lgForceUpdate)
Definition: temp_change.cpp:150
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
STATIC
#define STATIC
Definition: cddefines.h:97
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1080
molecule::index
int index
Definition: mole.h:169
hydro_bauman.h
t_iso_sp::tr
TransitionList * tr
Definition: iso.h:454
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:403
multi_arr
Definition: container_classes.h:941
chL
char chL[21]
Definition: iso_create.cpp:37
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
ipH3s
const int ipH3s
Definition: iso.h:30
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
L_
#define L_(A_)
Definition: iso.h:21
t_iso_sp::numLevels_malloc
long int numLevels_malloc
Definition: iso.h:502
thirdparty.h
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_isoCTRL::lgRandErrGen
bool lgRandErrGen[NISO]
Definition: iso.h:403
phycon
t_phycon phycon
Definition: phycon.cpp:6
iso_satellite_update
void iso_satellite_update(long nelem)
Definition: iso_create.cpp:1381
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:101
t_iso_sp::CascadeProb
multi_arr< double, 2 > CascadeProb
Definition: iso.h:450
ProxyIterator
Definition: proxy_iterator.h:58
t_isoCTRL::nLyman_malloc
long int nLyman_malloc[NISO]
Definition: iso.h:336
ex
static double * ex
Definition: species2.cpp:28
iso_create
void iso_create(void)
Definition: iso_create.cpp:39
lines_service.h
t_iso_sp
Definition: iso.h:441
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
PI4
const UNUSED double PI4
Definition: physconst.h:35
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
iso_recomb_malloc
void iso_recomb_malloc(void)
Definition: iso_radiative_recomb.cpp:765
opac
t_opac opac
Definition: opacity.cpp:5
atmdat.h
POW2
#define POW2
Definition: cddefines.h:929
t_iso_sp::ex
multi_arr< extra_tr, 2 > ex
Definition: iso.h:449
t_isoCTRL::lgFSM
bool lgFSM[NISO]
Definition: iso.h:399
MIN2
#define MIN2
Definition: cddefines.h:761
t_isoCTRL::ipLyaRedist
int ipLyaRedist[NISO]
Definition: iso.h:374
ATOMIC_MASS_UNIT
const UNUSED double ATOMIC_MASS_UNIT
Definition: physconst.h:88
t_iso_sp::bnl_effective
multi_arr< double, 3 > bnl_effective
Definition: iso.h:566
max_num_levels
long int max_num_levels
Definition: iso.cpp:10
HION_LTE_POP
const UNUSED double HION_LTE_POP
Definition: physconst.h:157
hunt_bisect
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:270
helike_einsta.h
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
t_hydro::HLineWidth
realnum HLineWidth
Definition: hydrogenic.h:63
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
iso_collapsed_bnl_print
void iso_collapsed_bnl_print(long ipISO, long nelem)
Definition: iso_create.cpp:1714
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
dense.h
MAX4
#define MAX4(a, b, c, d)
Definition: cddefines.h:792
mole
t_mole_local mole
Definition: mole.cpp:7
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
t_iso_sp::ipTrans
multi_arr< long, 2 > ipTrans
Definition: iso.h:448
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:365
FillExtraLymanLine
STATIC void FillExtraLymanLine(const TransitionList::iterator &t, long ipISO, long nelem, long nHi)
Definition: iso_create.cpp:977
iso_assign_quantum_numbers
STATIC void iso_assign_quantum_numbers(void)
Definition: iso_create.cpp:700
cddefines.h
null_mole
molecule * null_mole
Definition: mole_species.cpp:64
WAVNRYD
const UNUSED double WAVNRYD
Definition: physconst.h:173
ipH2s
const int ipH2s
Definition: iso.h:28
iso_collapsed_Aul_update
void iso_collapsed_Aul_update(long ipISO, long nelem)
Definition: iso_create.cpp:1760
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
iso_cascade
void iso_cascade(long ipISO, long nelem)
Definition: iso_create.cpp:1101
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
H_BAR
const UNUSED double H_BAR
Definition: physconst.h:144
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
iso_update_num_levels
void iso_update_num_levels(long ipISO, long nelem)
Definition: iso_create.cpp:1488
multi_arr::alloc
void alloc()
Definition: container_classes.h:1116
IPRAD
#define IPRAD
Definition: iso.h:86
t_isoCTRL::stat_ion
realnum stat_ion[NISO]
Definition: iso.h:362
FINE_STRUCTURE
const UNUSED double FINE_STRUCTURE
Definition: physconst.h:216
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
iso_get_total_num_levels
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
Definition: iso_create.cpp:1465
iso_satellite
STATIC void iso_satellite(void)
Definition: iso_create.cpp:1328
SPEEDLIGHT
const UNUSED double SPEEDLIGHT
Definition: physconst.h:100
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:141
ELECTRON_MASS
const UNUSED double ELECTRON_MASS
Definition: physconst.h:91
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
MAX2
#define MAX2
Definition: cddefines.h:782
t_isoCTRL::ipResoRedist
int ipResoRedist[NISO]
Definition: iso.h:374
LIMELM
const int LIMELM
Definition: cddefines.h:258
RYDLAM
const UNUSED double RYDLAM
Definition: physconst.h:176
RYD_INF
const UNUSED double RYD_INF
Definition: physconst.h:115
ipELECTRON
@ ipELECTRON
Definition: collision.h:9
t_isoCTRL::ipSubRedist
int ipSubRedist[NISO]
Definition: iso.h:374
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
M
static const int M
Definition: thirdparty.cpp:2815
t_phycon::te_ryd
double te_ryd
Definition: phycon.h:17
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:487
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:334
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
abscf
double abscf(double gf, double enercm, double gl)
Definition: lines_service.cpp:122
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:567
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
ipLY_A
const int ipLY_A
Definition: cddefines.h:296
hydrogenic.h
ipH2p
const int ipH2p
Definition: iso.h:29
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:894
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:664
TransitionList
Definition: transition.h:274
ipH3p
const int ipH3p
Definition: iso.h:31
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
COLL_CONST
const UNUSED double COLL_CONST
Definition: physconst.h:229
physconst.h
nLine
static long int nLine
Definition: save_line.cpp:288
SatelliteLines
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:38
t_iso_sp::BranchRatio
multi_arr< double, 2 > BranchRatio
Definition: iso.h:451
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
MIN4
#define MIN4(a, b, c, d)
Definition: cddefines.h:771
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
fixit
void fixit(void)
Definition: service.cpp:991
t_phycon::alogte
double alogte
Definition: phycon.h:82
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
molecule
Definition: mole.h:132
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
ERG1CM
const UNUSED double ERG1CM
Definition: physconst.h:164
phycon.h
t_iso_sp::n_HighestResolved_local
long int n_HighestResolved_local
Definition: iso.h:507
t_phycon::te32
double te32
Definition: phycon.h:49
iso_recomb_auxiliary_free
void iso_recomb_auxiliary_free(void)
Definition: iso_radiative_recomb.cpp:825
ipH1s
const int ipH1s
Definition: iso.h:27
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat.h:329
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1132
helike_energy
double helike_energy(long nelem, long ipLev)
ipSatelliteLines
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:37
lgHydroMalloc
bool lgHydroMalloc
Definition: cdinit.cpp:61
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
S_
#define S_(A_)
Definition: iso.h:22
dsexp
double dsexp(double x)
Definition: service.cpp:953
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
EmissionProxy::opacity
realnum & opacity() const
Definition: emission.h:593
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
iso_state_lifetime
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
Definition: iso_create.cpp:1057
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
ipCRDW
const int ipCRDW
Definition: cddefines.h:294
t_phycon::te
double te
Definition: phycon.h:11
iso_recomb_setup
void iso_recomb_setup(long ipISO)
Definition: iso_radiative_recomb.cpp:838
NISO
const int NISO
Definition: cddefines.h:261
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
ipH3d
const int ipH3d
Definition: iso.h:32
HeCSInterp
realnum HeCSInterp(long int nelem, long int ipHi, long int ipLo, long int Collider)
Definition: helike_cs.cpp:227
hydro_transprob
realnum hydro_transprob(long nelem, long ipHi, long ipLo)
Definition: hydroeinsta.cpp:45
hydroeinsta.h
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
HeCollidSetup
void HeCollidSetup(void)
Definition: helike_cs.cpp:90
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
iso_zero
STATIC void iso_zero(void)
Definition: iso_create.cpp:413
iso_collapsed_bnl_set
void iso_collapsed_bnl_set(long ipISO, long nelem)
Definition: iso_create.cpp:1519
POW4
#define POW4
Definition: cddefines.h:943
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
iso_allocate
STATIC void iso_allocate(void)
Definition: iso_create.cpp:457
TransitionProxy::ipEmis
int & ipEmis() const
Definition: transition.h:416
N_
#define N_(A_)
Definition: iso.h:20
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17