cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  return;
407 }
408 
409 /* ============================================================================== */
410 STATIC void iso_zero(void)
411 {
412  DEBUG_ENTRY( "iso_zero()" );
413 
414  hydro.HLineWidth = 0.;
415 
416  /****************************************************/
417  /********** initialize some variables **********/
418  /****************************************************/
419  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
420  {
421  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
422  {
423  if( nelem < 2 || dense.lgElmtOn[nelem] )
424  {
425  for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
426  {
427  iso_sp[ipISO][nelem].st[ipHi].Pop() = 0.;
428  iso_sp[ipISO][nelem].fb[ipHi].Reset();
429  }
430  if (ipISO <= nelem)
431  iso_sp[ipISO][nelem].st[0].Pop() =
432  dense.xIonDense[nelem][nelem-ipISO];
433  }
434 
435  if( nelem < 2 )
436  {
437  iso_collapsed_bnl_set( ipISO, nelem );
438  //iso_collapsed_bnl_print( ipISO, nelem );
439  iso_collapsed_Aul_update( ipISO, nelem );
440  iso_collapsed_lifetimes_update( ipISO, nelem );
441  }
442  }
443  }
444 
445  /* ground state of H and He is different since totally determine
446  * their own opacities */
447  iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ConOpacRatio = 1e-5;
448  iso_sp[ipH_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
449  iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
450 
451  return;
452 }
453 
454 STATIC void iso_allocate(void)
455 {
456 
457  DEBUG_ENTRY( "iso_allocate()" );
458 
459  /* the hydrogen and helium like iso-sequences */
460  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
461  {
462  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
463  {
464  /* only grab core for elements that are turned on */
465  if( nelem < 2 || dense.lgElmtOn[nelem] )
466  {
467  t_iso_sp *sp = &iso_sp[ipISO][nelem];
468  sp->numLevels_malloc = sp->numLevels_max;
469 
470  ASSERT( sp->numLevels_max > 0 );
471  ASSERT( iso_ctrl.nLyman_malloc[ipISO] == iso_ctrl.nLyman[ipISO] );
472 
473  sp->CachedAs.reserve( MAX2(1, sp->nCollapsed_max) );
474 
475  sp->ipTrans.reserve( sp->numLevels_max );
476  sp->ex.reserve( sp->numLevels_max );
477  sp->CascadeProb.reserve( sp->numLevels_max );
478  sp->BranchRatio.reserve( sp->numLevels_max );
479  //sp->st.resize( sp->numLevels_max );
480  sp->fb.resize( sp->numLevels_max );
481 
482  for( long i = 0; i < sp->nCollapsed_max; ++i )
483  {
484  sp->CachedAs.reserve( i, sp->numLevels_max - sp->nCollapsed_max );
485  for( long i1 = 0; i1 < sp->numLevels_max - sp->nCollapsed_max; ++i1 )
486  {
487  /* allocate two spaces delta L +/- 1 */
488  sp->CachedAs.reserve( i, i1, 2 );
489  }
490  }
491 
492  sp->QuantumNumbers2Index.reserve( sp->n_HighestResolved_max + sp->nCollapsed_max + 1 );
493 
494  for( long i = 1; i <= sp->n_HighestResolved_max + sp->nCollapsed_max; ++i )
495  {
496  /* Allocate proper number of angular momentum quantum number. */
497  sp->QuantumNumbers2Index.reserve( i, i );
498 
499  for( long i1 = 0; i1 < i; ++i1 )
500  {
501  /* This may have to change for other iso sequences. */
502  ASSERT( NISO == 2 );
503  /* Allocate 4 spaces for multiplicity. H-like will be accessed with "2" for doublet,
504  * He-like will be accessed via "1" for singlet or "3" for triplet. "0" will not be used. */
505  sp->QuantumNumbers2Index.reserve( i, i1, 4 );
506  }
507  }
508 
509  for( long n=1; n < sp->numLevels_max; ++n )
510  {
511  sp->ipTrans.reserve( n, n );
512  }
513 
514  for( long n=0; n < sp->numLevels_max; ++n )
515  {
516  sp->ex.reserve( n, sp->numLevels_max );
517  sp->CascadeProb.reserve( n, sp->numLevels_max );
518  sp->BranchRatio.reserve( n, sp->numLevels_max );
519  }
520 
521  sp->ipTrans.alloc();
522  sp->ex.alloc();
523  sp->CascadeProb.alloc();
524  sp->BranchRatio.alloc();
525 
526  sp->CachedAs.alloc();
527  sp->QuantumNumbers2Index.alloc();
528  sp->bnl_effective.alloc( sp->QuantumNumbers2Index.clone() );
529  sp->QuantumNumbers2Index.invalidate();
530  sp->bnl_effective.invalidate();
531  }
532  }
533  }
534 
535  ipSatelliteLines.reserve( NISO );
536  ipExtraLymanLines.reserve( NISO );
537 
538  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
539  {
540  ipSatelliteLines.reserve( ipISO, LIMELM );
541  ipExtraLymanLines.reserve( ipISO, LIMELM );
542 
543  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
544  {
545  /* only grab core for elements that are turned on */
546  if( nelem < 2 || dense.lgElmtOn[nelem] )
547  {
548  ASSERT( iso_sp[ipISO][nelem].numLevels_max > 0 );
549 
550  ipSatelliteLines.reserve( ipISO, nelem, iso_sp[ipISO][nelem].numLevels_max );
551  ipExtraLymanLines.reserve( ipISO, nelem, iso_ctrl.nLyman_malloc[ipISO] );
552  }
553  }
554  }
555 
556  ipSatelliteLines.alloc();
557  ipExtraLymanLines.alloc();
558 
559  Transitions.resize(NISO);
560  SatelliteLines.resize(NISO);
561  ExtraLymanLines.resize(NISO);
562  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
563  {
564  Transitions[ipISO].reserve(LIMELM);
565  SatelliteLines[ipISO].reserve(LIMELM);
566  ExtraLymanLines[ipISO].reserve(LIMELM);
567  for( long nelem=0; nelem < ipISO; ++nelem )
568  {
569  Transitions[ipISO].push_back(
570  TransitionList("Insanity",&AnonStates));
571  SatelliteLines[ipISO].push_back(
572  TransitionList("Insanity",&AnonStates));
573  ExtraLymanLines[ipISO].push_back(
574  TransitionList("Insanity",&AnonStates));
575  }
576  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
577  {
578  if( nelem < 2 || dense.lgElmtOn[nelem] )
579  {
580  Transitions[ipISO].push_back(
581  TransitionList("Isosequence",&iso_sp[ipISO][nelem].st));
582  SatelliteLines[ipISO].push_back(
583  TransitionList("SatelliteLines",&iso_sp[ipISO][nelem].st));
584  ExtraLymanLines[ipISO].push_back(
585  TransitionList("ExtraLymanLines",&iso_sp[ipISO][nelem].st));
586  }
587  else
588  {
589  Transitions[ipISO].push_back(
590  TransitionList("Insanity",&AnonStates));
591  SatelliteLines[ipISO].push_back(
592  TransitionList("Insanity",&AnonStates));
593  ExtraLymanLines[ipISO].push_back(
594  TransitionList("Insanity",&AnonStates));
595  }
596  }
597  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
598  {
599  /* only grab core for elements that are turned on */
600  if( nelem < 2 || dense.lgElmtOn[nelem] )
601  {
602  if( iso_ctrl.lgDielRecom[ipISO] )
603  {
604  SatelliteLines[ipISO][nelem].resize( iso_sp[ipISO][nelem].numLevels_max );
605  AllTransitions.push_back(SatelliteLines[ipISO][nelem]);
606  unsigned int nLine = 0;
607  for( long ipLo=0; ipLo<iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
608  {
609  /* Upper level is continuum, use a generic state
610  * lower level is the same as the index. */
611  ipSatelliteLines[ipISO][nelem][ipLo] = nLine;
612  SatelliteLines[ipISO][nelem][nLine].Junk();
613  long ipHi = iso_sp[ipISO][nelem].numLevels_max;
614  SatelliteLines[ipISO][nelem][nLine].setHi(ipHi);
615  SatelliteLines[ipISO][nelem][nLine].setLo(ipLo);
616  SatelliteLines[ipISO][nelem][nLine].AddLine2Stack();
617  ++nLine;
618  }
619  ASSERT(SatelliteLines[ipISO][nelem].size() == nLine);
620  }
621 
622  //iso_sp[ipISO][nelem].tr.resize( iso_sp[ipISO][nelem].ipTrans.size() );
623  //iso_sp[ipISO][nelem].tr.states() = &iso_sp[ipISO][nelem].st;
624  Transitions[ipISO][nelem].resize( iso_sp[ipISO][nelem].ipTrans.size() );
625  AllTransitions.push_back(Transitions[ipISO][nelem]);
626  unsigned int nTransition=0;
627  for( long ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
628  {
629  for( long ipLo=0; ipLo < ipHi; ipLo++ )
630  {
631  /* set ENTIRE array to impossible values, in case of bad pointer */
632  iso_sp[ipISO][nelem].ipTrans[ipHi][ipLo] = nTransition;
633  Transitions[ipISO][nelem][nTransition].Junk();
634  Transitions[ipISO][nelem][nTransition].setHi(ipHi);
635  Transitions[ipISO][nelem][nTransition].setLo(ipLo);
636  ++nTransition;
637  }
638  }
639  ASSERT(Transitions[ipISO][nelem].size() == nTransition);
640  iso_sp[ipISO][nelem].tr = &Transitions[ipISO][nelem];
641 
642  /* junk the extra Lyman lines */
643  AllTransitions.push_back(ExtraLymanLines[ipISO][nelem]);
644  ExtraLymanLines[ipISO][nelem].resize(iso_ctrl.nLyman_malloc[ipISO]-2);
645  ExtraLymanLines[ipISO][nelem].states() = &iso_sp[ipISO][nelem].st;
646  unsigned int nExtraLyman = 0;
647  for( long ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
648  {
649  ipExtraLymanLines[ipISO][nelem][ipHi] = nExtraLyman;
650  ExtraLymanLines[ipISO][nelem][nExtraLyman].Junk();
651  long ipHi_offset = iso_sp[ipISO][nelem].numLevels_max + ipHi - 2;
652  if( iso_ctrl.lgDielRecom[ipISO] )
653  ipHi_offset += 1;
654  ExtraLymanLines[ipISO][nelem][nExtraLyman].setHi(ipHi_offset);
655  /* lower level is just ground state of the ion */
656  ExtraLymanLines[ipISO][nelem][nExtraLyman].setLo(0);
657  ExtraLymanLines[ipISO][nelem][nExtraLyman].AddLine2Stack();
658  ++nExtraLyman;
659  }
660  ASSERT(ExtraLymanLines[ipISO][nelem].size() == nExtraLyman);
661  }
662  }
663  }
664 
665  // associate line and level stacks with species
666  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
667  {
668  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
669  {
670  if( dense.lgElmtOn[nelem] )
671  {
672  long ion = nelem - ipISO;
673  ASSERT( ion >= 0 && ion <= nelem );
674  char chLabel[6] = {'\0'}, chTemp[4] = {'\0'};
675  sprintf( chLabel, "%s", elementnames.chElementSym[nelem] );
676  if( chLabel[1]==' ' )
677  chLabel[1] = '\0';
678  else
679  chLabel[2] = '\0';
680  if( ion==1 )
681  sprintf( chTemp, "+" );
682  else if( ion>0 )
683  sprintf( chTemp, "+%li", ion );
684  strcat( chLabel, chTemp );
685 
686  molecule *spmole = findspecies(chLabel);
687  ASSERT( spmole != null_mole );
688  mole.species[ spmole->index ].levels = &iso_sp[ipISO][nelem].st;
689  mole.species[ spmole->index ].lines = &Transitions[ipISO][nelem];
690  }
691  }
692  }
693 
694  return;
695 }
696 
697 STATIC void iso_assign_quantum_numbers(void)
698 {
699  long int
700  ipLo,
701  level,
702  i,
703  in,
704  il,
705  is,
706  ij;
707 
708  DEBUG_ENTRY( "iso_assign_quantum_numbers()" );
709 
710  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
711  {
712  long ipISO = ipH_LIKE;
713  /* only check elements that are turned on */
714  if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
715  {
716  i = 0;
717 
718  /* 2 for doublet */
719  is = ipDOUBLET;
720 
721  /* this loop is over quantum number n */
722  for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
723  {
724  for( il = 0L; il < in; ++il )
725  {
726  iso_sp[ipISO][nelem].st[i].n() = in;
727  iso_sp[ipISO][nelem].st[i].S() = is;
728  iso_sp[ipISO][nelem].st[i].l() = il;
729  iso_sp[ipISO][nelem].st[i].j() = -1;
730  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
731  ++i;
732  }
733  }
734  /* now do the collapsed levels */
735  in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
736  for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
737  {
738  iso_sp[ipISO][nelem].st[level].n() = in;
739  iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
740  iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
741  iso_sp[ipISO][nelem].st[level].j() = -1;
742  /* Point every l to same index for collapsed levels. */
743  for( il = 0; il < in; ++il )
744  {
745  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
746  }
747  ++in;
748  }
749  --in;
750 
751  /* confirm that we did not overrun the array */
752  ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
753 
754  /* confirm that n is positive and not greater than the max n. */
755  ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
756 
757  /* Verify states and QuantumNumbers2Index agree in all cases */
758  for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
759  {
760  for( il = 0L; il < in; ++il )
761  {
762  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
763  if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
764  {
765  /* Must only check these for resolved levels...
766  * collapsed levels have pointers for l and s that will blow if used. */
767  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
768  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
769  }
770  }
771  }
772  }
773  }
774 
775  /* then do he-like */
776  for( long nelem=ipHELIUM; nelem < LIMELM; nelem++ )
777  {
778  long ipISO = ipHE_LIKE;
779  /* only check elements that are turned on */
780  if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
781  {
782  i = 0;
783 
784  /* this loop is over quantum number n */
785  for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
786  {
787  for( il = 0L; il < in; ++il )
788  {
789  for( is = 3L; is >= 1L; is -= 2 )
790  {
791  /* All levels except singlet P follow the ordering scheme: */
792  /* lower l's have lower energy */
793  /* triplets have lower energy */
794  if( (il == 1L) && (is == 1L) )
795  continue;
796  /* n = 1 has no triplet, of course. */
797  if( (in == 1L) && (is == 3L) )
798  continue;
799 
800  /* singlets */
801  if( is == 1 )
802  {
803  iso_sp[ipISO][nelem].st[i].n() = in;
804  iso_sp[ipISO][nelem].st[i].S() = is;
805  iso_sp[ipISO][nelem].st[i].l() = il;
806  /* this is not a typo, J=L for singlets. */
807  iso_sp[ipISO][nelem].st[i].j() = il;
808  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
809  ++i;
810  }
811  /* 2 triplet P is j-resolved */
812  else if( (in == 2) && (il == 1) && (is == 3) )
813  {
814  ij = 0;
815  do
816  {
817  iso_sp[ipISO][nelem].st[i].n() = in;
818  iso_sp[ipISO][nelem].st[i].S() = is;
819  iso_sp[ipISO][nelem].st[i].l() = il;
820  iso_sp[ipISO][nelem].st[i].j() = ij;
821  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
822  ++i;
823  ++ij;
824  /* repeat this for the separate j-levels within 2^3P. */
825  } while ( ij < 3 );
826  }
827  else
828  {
829  iso_sp[ipISO][nelem].st[i].n() = in;
830  iso_sp[ipISO][nelem].st[i].S() = is;
831  iso_sp[ipISO][nelem].st[i].l() = il;
832  iso_sp[ipISO][nelem].st[i].j() = -1L;
833  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
834  ++i;
835  }
836  }
837  }
838  /* Insert singlet P at the end of every sequence for a given n. */
839  if( in > 1L )
840  {
841  iso_sp[ipISO][nelem].st[i].n() = in;
842  iso_sp[ipISO][nelem].st[i].S() = 1L;
843  iso_sp[ipISO][nelem].st[i].l() = 1L;
844  iso_sp[ipISO][nelem].st[i].j() = 1L;
845  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][1][1] = i;
846  ++i;
847  }
848  }
849  /* now do the collapsed levels */
850  in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
851  for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
852  {
853  iso_sp[ipISO][nelem].st[level].n() = in;
854  iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
855  iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
856  iso_sp[ipISO][nelem].st[level].j() = -1;
857  /* Point every l and s to same index for collapsed levels. */
858  for( il = 0; il < in; ++il )
859  {
860  for( is = 1; is <= 3; is += 2 )
861  {
862  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
863  }
864  }
865  ++in;
866  }
867  --in;
868 
869  /* confirm that we did not overrun the array */
870  ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
871 
872  /* confirm that n is positive and not greater than the max n. */
873  ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
874 
875  /* Verify states and QuantumNumbers2Index agree in all cases */
876  for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
877  {
878  for( il = 0L; il < in; ++il )
879  {
880  for( is = 3L; is >= 1; is -= 2 )
881  {
882  /* Ground state is not triplicate. */
883  if( (in == 1L) && (is == 3L) )
884  continue;
885 
886  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
887  if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
888  {
889  /* Must only check these for resolved levels...
890  * collapsed levels have pointers for l and s that will blow if used. */
891  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
892  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
893  }
894  }
895  }
896  }
897  }
898  }
899 
900  for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
901  {
902  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
903  {
904  /* must always do helium even if turned off */
905  if( nelem < 2 || dense.lgElmtOn[nelem] )
906  {
907  for( ipLo=ipH1s; ipLo < iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
908  {
909  iso_sp[ipISO][nelem].st[ipLo].nelem() = (int)(nelem+1);
910  iso_sp[ipISO][nelem].st[ipLo].IonStg() = (int)(nelem+1-ipISO);
911 
912  if( iso_sp[ipISO][nelem].st[ipLo].j() >= 0 )
913  {
914  iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*iso_sp[ipISO][nelem].st[ipLo].j()+1.f;
915  }
916  else if( iso_sp[ipISO][nelem].st[ipLo].l() >= 0 )
917  {
918  iso_sp[ipISO][nelem].st[ipLo].g() = (2.f*iso_sp[ipISO][nelem].st[ipLo].l()+1.f) *
919  iso_sp[ipISO][nelem].st[ipLo].S();
920  }
921  else
922  {
923  if( ipISO == ipH_LIKE )
924  iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
925  else if( ipISO == ipHE_LIKE )
926  iso_sp[ipISO][nelem].st[ipLo].g() = 4.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
927  else
928  {
929  /* replace this with correct thing if more sequences are added. */
930  TotalInsanity();
931  }
932  }
933  char chConfiguration[11] = " ";
934  long nCharactersWritten = 0;
935 
936  ASSERT( iso_sp[ipISO][nelem].st[ipLo].n() < 1000 );
937 
938  /* include j only if defined. */
939  if( iso_sp[ipISO][nelem].st[ipLo].n() > iso_sp[ipISO][nelem].n_HighestResolved_max )
940  {
941  nCharactersWritten = sprintf( chConfiguration, "n=%3li",
942  iso_sp[ipISO][nelem].st[ipLo].n() );
943  }
944  else if( iso_sp[ipISO][nelem].st[ipLo].j() > 0 )
945  {
946  nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c_%li",
947  iso_sp[ipISO][nelem].st[ipLo].n(),
948  iso_sp[ipISO][nelem].st[ipLo].S(),
949  chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l() ) ],
950  iso_sp[ipISO][nelem].st[ipLo].j() );
951  }
952  else
953  {
954  nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c",
955  iso_sp[ipISO][nelem].st[ipLo].n(),
956  iso_sp[ipISO][nelem].st[ipLo].S(),
957  chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l()) ] );
958  }
959 
960  ASSERT( nCharactersWritten <= 10 );
961  chConfiguration[10] = '\0';
962 
963  strncpy( iso_sp[ipISO][nelem].st[ipLo].chConfig(), chConfiguration, 10 );
964  }
965  }
966  }
967  }
968  return;
969 }
970 
971 #if defined(__ICC) && defined(__i386)
972 #pragma optimization_level 1
973 #endif
974 STATIC void FillExtraLymanLine( const TransitionList::iterator& t, long ipISO, long nelem, long nHi )
975 {
976  double Enerwn, Aul;
977 
978  DEBUG_ENTRY( "FillExtraLymanLine()" );
979 
980  /* atomic number or charge and stage: */
981  (*(*t).Hi()).nelem() = (int)(nelem+1);
982  (*(*t).Hi()).IonStg() = (int)(nelem+1-ipISO);
983 
984  (*(*t).Hi()).n() = nHi;
985 
986  /* statistical weight is same as statistical weight of corresponding LyA. */
987  (*(*t).Hi()).g() = iso_sp[ipISO][nelem].st[iso_ctrl.nLyaLevel[ipISO]].g();
988 
989  /* energies */
990  Enerwn = iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd * RYD_INF * ( 1. - 1./POW2((double)nHi) );
991 
992  /* transition energy in various units:*/
993  (*t).EnergyWN() = (realnum)(Enerwn);
994  (*t).WLAng() = (realnum)(1.0e8/ Enerwn/ RefIndex(Enerwn));
995  (*(*t).Hi()).energy().set( Enerwn, "cm^-1" );
996 
997  if( ipISO == ipH_LIKE )
998  {
999  Aul = H_Einstein_A( nHi, 1, 1, 0, nelem+1 );
1000  }
1001  else
1002  {
1003  if( nelem == ipHELIUM )
1004  {
1005  /* A simple fit for the calculation of Helium lyman Aul's. */
1006  Aul = (1.508e10) / pow((double)nHi,2.975);
1007  }
1008  else
1009  {
1010  /* Fit to values given in
1011  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
1012  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
1013  /* originally astro.ph. 0201454 */
1014  Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
1015  }
1016  }
1017 
1018  (*t).Emis().Aul() = (realnum)Aul;
1019 
1020  (*(*t).Hi()).lifetime() = iso_state_lifetime( ipISO, nelem, nHi, 1 );
1021 
1022  (*t).Emis().dampXvel() = (realnum)( 1.f / (*(*t).Hi()).lifetime() / PI4 / (*t).EnergyWN() );
1023 
1024  (*t).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
1025 
1026  (*t).Emis().gf() = (realnum)(GetGF((*t).Emis().Aul(), (*t).EnergyWN(), (*(*t).Hi()).g()));
1027 
1028  /* derive the abs coef, call to function is Emis().gf(), wl (A), g_low */
1029  (*t).Emis().opacity() = (realnum)(abscf((*t).Emis().gf(), (*t).EnergyWN(), (*(*t).Lo()).g()));
1030 
1031  /* create array indices that will blow up */
1032  (*t).ipCont() = INT_MIN;
1033  (*t).Emis().ipFine() = INT_MIN;
1034 
1035  {
1036  /* option to print particulars of some line when called
1037  * a prettier print statement is near where chSpin is defined below
1038  * search for "pretty print" */
1039  enum {DEBUG_LOC=false};
1040  if( DEBUG_LOC )
1041  {
1042  fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
1043  nelem+1,
1044  nHi,
1045  (*t).Emis().Aul() ,
1046  (*t).Emis().opacity()
1047  );
1048  }
1049  }
1050  return;
1051 }
1052 
1053 /* calculate radiative lifetime of an individual iso state */
1054 double iso_state_lifetime( long ipISO, long nelem, long n, long l )
1055 {
1056  /* >>refer hydro lifetimes Horbatsch, M. W., Horbatsch, M. and Hessels, E. A. 2005, JPhysB, 38, 1765 */
1057 
1058  double tau, t0, eps2;
1059  /* mass of electron */
1060  double m = ELECTRON_MASS;
1061  /* nuclear mass */
1062  double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT;
1063  double mu = (m*M)/(M+m);
1064  long z = 1;
1065  long Z = nelem + 1 - ipISO;
1066 
1067  DEBUG_ENTRY( "iso_state_lifetime()" );
1068 
1069  /* this should not be used for l=0 per the Horbatsch et al. paper */
1070  ASSERT( l > 0 );
1071 
1072  eps2 = 1. - ( l*l + l + 8./47. - (l+1.)/69./n ) / POW2( (double)n );
1073 
1074  t0 = 3. * H_BAR * pow( (double)n, 5.) /
1075  ( 2. * POW4( (double)( z * Z ) ) * pow( FINE_STRUCTURE, 5. ) * mu * POW2( SPEEDLIGHT ) ) *
1076  POW2( (m + M)/(Z*m + z*M) );
1077 
1078  tau = t0 * ( 1. - eps2 ) /
1079  ( 1. + 19./88.*( (1./eps2 - 1.) * log( 1. - eps2 ) + 1. -
1080  0.5 * eps2 - 0.025 * eps2 * eps2 ) );
1081 
1082  if( ipISO == ipHE_LIKE )
1083  {
1084  /* iso_state_lifetime is not spin specific, must exclude helike triplet here. */
1085  tau /= 3.;
1086  /* this is also necessary to correct the helike lifetimes */
1087  tau *= 1.1722 * pow( (double)nelem, 0.1 );
1088  }
1089 
1090  /* would probably need a new lifetime algorithm for any other iso sequences. */
1091  ASSERT( ipISO <= ipHE_LIKE );
1092  ASSERT( tau > 0. );
1093 
1094  return tau;
1095 }
1096 
1097 /* calculate cascade probabilities, branching ratios, and associated errors. */
1098 void iso_cascade( long ipISO, long nelem )
1099 {
1100  /* The sum of all A's coming out of a given n,
1101  * Below we assert a monotonic trend. */
1102  double *SumAPerN;
1103 
1104  long int i, j, ipLo, ipHi;
1105 
1106  DEBUG_ENTRY( "iso_cascade()" );
1107 
1108  SumAPerN = ((double*)MALLOC( (size_t)(iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double )));
1109  memset( SumAPerN, 0, (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double ) );
1110 
1111  /* Initialize some ground state stuff, easier here than in loops. */
1112  iso_sp[ipISO][nelem].CascadeProb[0][0] = 1.;
1113  if( iso_ctrl.lgRandErrGen[ipISO] )
1114  {
1115  iso_sp[ipISO][nelem].fb[0].SigmaAtot = 0.;
1116  iso_sp[ipISO][nelem].ex[0][0].SigmaCascadeProb = 0.;
1117  }
1118 
1119  /***************************************************************************/
1120  /****** Cascade probabilities, Branching ratios, and associated errors *****/
1121  /***************************************************************************/
1122  for( ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1123  {
1124  double SumAs = 0.;
1125 
1131  /* initialize variables. */
1132  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipHi] = 1.;
1133  iso_sp[ipISO][nelem].CascadeProb[ipHi][0] = 0.;
1134  iso_sp[ipISO][nelem].BranchRatio[ipHi][0] = 0.;
1135 
1136  if( iso_ctrl.lgRandErrGen[ipISO] )
1137  {
1138  iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = 0.;
1139  iso_sp[ipISO][nelem].ex[ipHi][ipHi].SigmaCascadeProb = 0.;
1140  }
1141 
1142  long ipLoStart = 0;
1143  if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) )
1144  ipLoStart = 1;
1145 
1146  for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1147  {
1148  SumAs += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1149  }
1150 
1151  for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1152  {
1153  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1154  {
1155  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1156  iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] = 0.;
1157  continue;
1158  }
1159 
1160  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1161  iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] =
1162  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() / SumAs;
1163 
1164  ASSERT( iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] <= 1.0000001 );
1165 
1166  SumAPerN[N_(ipHi)] += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1167 
1168  /* there are some negative energy transitions, where the order
1169  * has changed, but these are not optically allowed, these are
1170  * same n, different L, forbidden transitions */
1171  ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() > 0. ||
1172  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA );
1173 
1174  if( iso_ctrl.lgRandErrGen[ipISO] )
1175  {
1176  ASSERT( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] >= 0. );
1177  /* Uncertainties in A's are added in quadrature, square root is taken below. */
1178  iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot +=
1179  pow( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] *
1180  (double)iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(), 2. );
1181  }
1182  }
1183 
1184  if( iso_ctrl.lgRandErrGen[ipISO] )
1185  {
1186  /* Uncertainties in A's are added in quadrature above, square root taken here. */
1187  iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = sqrt( iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot );
1188  }
1189 
1190  /* cascade probabilities */
1191  for( i=0; i<ipHi; i++ )
1192  {
1193  for( ipLo=0; ipLo<=i; ipLo++ )
1194  {
1195  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] += iso_sp[ipISO][nelem].BranchRatio[ipHi][i] * iso_sp[ipISO][nelem].CascadeProb[i][ipLo];
1196  }
1197  }
1198 
1199  if( iso_ctrl.lgRandErrGen[ipISO] )
1200  {
1201  for( ipLo=0; ipLo<ipHi; ipLo++ )
1202  {
1203  double SigmaCul = 0.;
1204  for( i=ipLo; i<ipHi; i++ )
1205  {
1206  if( iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul() > iso_ctrl.SmallA )
1207  {
1208  /* Uncertainties in A's and cascade probabilities */
1209  double SigmaA = iso_sp[ipISO][nelem].ex[ipHi][i].Error[IPRAD] *
1210  iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul();
1211  SigmaCul +=
1212  pow(SigmaA*iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
1213  pow(iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*iso_sp[ipISO][nelem].BranchRatio[ipHi][i]*
1214  iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
1215  pow(iso_sp[ipISO][nelem].ex[i][ipLo].SigmaCascadeProb*iso_sp[ipISO][nelem].BranchRatio[ipHi][i], 2.);
1216  }
1217  }
1218  SigmaCul = sqrt(SigmaCul);
1219  iso_sp[ipISO][nelem].ex[ipHi][ipLo].SigmaCascadeProb = SigmaCul;
1220  }
1221  }
1222  }
1223 
1224  /************************************************************************/
1225  /*** Allowed decay conversion probabilities. See Robbins68b, Table 1. ***/
1226  /************************************************************************/
1227  {
1228  enum {DEBUG_LOC=false};
1229 
1230  if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) )
1231  {
1232  /* To output Bm(n,l; ipLo), set ipLo, hi_l, and hi_s accordingly. */
1233  long int hi_l,hi_s;
1234  double Bm;
1235 
1236  /* these must be set for following output to make sense
1237  * as is, a dangerous bit of code - set NaN for safety */
1238  hi_s = -100000;
1239  hi_l = -100000;
1240  ipLo = -100000;
1241  /* tripS to 2^3P */
1242  //hi_l = 0, hi_s = 3, ipLo = ipHe2p3P0;
1243 
1244  /* tripD to 2^3P */
1245  //hi_l = 2, hi_s = 3, ipLo = ipHe2p3P0;
1246 
1247  /* tripP to 2^3S */
1248  //hi_l = 1, hi_s = 3, ipLo = ipHe2s3S;
1249 
1250  ASSERT( hi_l != iso_sp[ipISO][nelem].st[ipLo].l() );
1251 
1252  fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
1253  fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");
1254 
1255  for( ipHi=ipHe2p3P2; ipHi<iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipHi++ )
1256  {
1257  /* Pick out excitations from metastable 2tripS to ntripP. */
1258  if( (iso_sp[ipISO][nelem].st[ipHi].l() == 1) && (iso_sp[ipISO][nelem].st[ipHi].S() == 3) )
1259  {
1260  fprintf(ioQQQ,"\n%ld\t",iso_sp[ipISO][nelem].st[ipHi].n());
1261  j = 0;
1262  Bm = 0;
1263  for( i = ipLo; i<=ipHi; i++)
1264  {
1265  if( (iso_sp[ipISO][nelem].st[i].l() == hi_l) && (iso_sp[ipISO][nelem].st[i].S() == hi_s) )
1266  {
1267  if( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
1268  {
1269  Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * ( iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P0] +
1270  iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P1] + iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P2] );
1271  }
1272  else
1273  Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * iso_sp[ipISO][nelem].BranchRatio[i][ipLo];
1274 
1275  if( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
1276  {
1277  j++;
1278  if(j == 3)
1279  {
1280  fprintf(ioQQQ,"%2.4e\t",Bm);
1281  Bm = 0;
1282  }
1283  }
1284  else
1285  {
1286  fprintf(ioQQQ,"%2.4e\t",Bm);
1287  Bm = 0;
1288  }
1289  }
1290  }
1291  }
1292  }
1293  fprintf(ioQQQ,"\n\n");
1294  }
1295  }
1296 
1297  /******************************************************/
1298  /*** Lifetimes should increase monotonically with ***/
1299  /*** increasing n...Make sure the A's decrease. ***/
1300  /******************************************************/
1301  for( i=2; i < iso_sp[ipISO][nelem].n_HighestResolved_max; ++i)
1302  {
1303  ASSERT( (SumAPerN[i] > SumAPerN[i+1]) || opac.lgCaseB );
1304  }
1305 
1306  {
1307  enum {DEBUG_LOC=false};
1308  if( DEBUG_LOC /* && (ipISO == ipH_LIKE) && (nelem == ipHYDROGEN) */)
1309  {
1310  for( i = 2; i<= (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max); ++i)
1311  {
1312  fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]);
1313  }
1314  }
1315  }
1316 
1317  free( SumAPerN );
1318 
1319  return;
1320 }
1321 
1323 /* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002 */
1324 /* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001. */
1325 STATIC void iso_satellite( void )
1326 {
1327  DEBUG_ENTRY( "iso_satellite()" );
1328 
1329  for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
1330  {
1331  for( long nelem = ipISO; nelem < LIMELM; nelem++ )
1332  {
1333  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
1334  {
1335  for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
1336  {
1337  char chLab[5]=" ";
1338 
1339  TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
1340  (*tr).Zero();
1341 
1342  /* Make approximation that all levels have energy of H-like 2s level */
1343  /* Lines to 1s2s have roughly energy of parent Ly-alpha. So lines to 1snL will have an energy
1344  * smaller by the difference between nL and 2s energies. Therefore, the following has
1345  * energy of parent Ly-alpha MINUS the difference between daughter level and daughter n=2 level. */
1346  (*tr).WLAng() = (realnum)(RYDLAM/
1347  ((iso_sp[ipISO-1][nelem].fb[0].xIsoLevNIonRyd - iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd) -
1348  (iso_sp[ipISO][nelem].fb[1].xIsoLevNIonRyd- iso_sp[ipISO][nelem].fb[i].xIsoLevNIonRyd)) );
1349 
1350  (*tr).EnergyWN() = 1.e8f /
1351  (*tr).WLAng();
1352 
1353  /* generate label for this ion */
1354  sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
1355 
1356  (*tr).Emis().iRedisFun() = ipCRDW;
1357  /* this is not the usual nelem, is it atomic not C scale. */
1358  (*(*tr).Hi()).nelem() = nelem + 1;
1359  (*(*tr).Hi()).IonStg() = nelem + 1 - ipISO;
1360  fixit(); /* what should the stat weight of the upper level be? For now say 2. */
1361  (*(*tr).Hi()).g() = 2.f;
1362  // The lower level must already be initialized.
1363  ASSERT( (*(*tr).Lo()).g() == iso_sp[ipISO][nelem].st[i].g() );
1364  //(*(*tr).Lo()).g() = iso_sp[ipISO][nelem].st[i].g();
1365  (*tr).Emis().PopOpc() =
1366  (*(*tr).Lo()).Pop();
1367 
1368  (*tr).Emis().pump() = 0.;
1369 
1370  }
1371  }
1372  }
1373  }
1374 
1375  return;
1376 }
1377 
1378 void iso_satellite_update( long nelem )
1379 {
1380  double ConBoltz, LTE_pop=SMALLFLOAT+FLT_EPSILON, factor1, ConvLTEPOP;
1381 
1382  DEBUG_ENTRY( "iso_satellite_update()" );
1383 
1384  for( long ipISO = ipHE_LIKE; ipISO < MIN2(NISO,nelem+1); ipISO++ )
1385  {
1386  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
1387  {
1388  for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
1389  {
1390  double dr_rate = iso_sp[ipISO][nelem].fb[i].DielecRecomb * iso_ctrl.lgDielRecom[ipISO];
1391 
1392  TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
1393  (*tr).Emis().phots() =
1394  dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
1395 
1396  (*tr).Emis().xIntensity() =
1397  (*tr).Emis().phots() *
1398  ERG1CM * (*tr).EnergyWN();
1399 
1400  /* We set line intensity above using a rate, but here we need a transition probability.
1401  * We can obtain this by dividing dr_rate by the population of the autoionizing level.
1402  * We assume this level is in statistical equilibrium. */
1403  factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
1404  (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
1405 
1406  /* term in () is stat weight of electron * ion */
1407  ConvLTEPOP = pow(factor1,1.5)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;
1408 
1409  /* This Boltzmann factor is exp( +ioniz energy / Te ). For simplicity, we make
1410  * the fair approximation that all of the autoionizing levels have an energy
1411  * equal to the parents n=2. */
1412  ConBoltz = dsexp(iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd/phycon.te_ryd);
1413 
1414  if( ConBoltz >= SMALLDOUBLE )
1415  {
1416  /* The energy used to calculate ConBoltz above
1417  * should be negative since this is above the continuum, but
1418  * to be safe we calculate ConBoltz with a positive energy above
1419  * and multiply by it here instead of dividing. */
1420  LTE_pop = (*(*tr).Hi()).g() * ConBoltz * ConvLTEPOP;
1421  }
1422 
1423  LTE_pop = max( LTE_pop, 1e-30f );
1424 
1425  /* Now the transition probability is simply dr_rate/LTE_pop. */
1426  (*tr).Emis().Aul() = (realnum)(dr_rate/LTE_pop);
1427  (*tr).Emis().Aul() =
1428  max( iso_ctrl.SmallA, (*tr).Emis().Aul() );
1429 
1430  (*tr).Emis().gf() = (realnum)GetGF(
1431  (*tr).Emis().Aul(),
1432  (*tr).EnergyWN(),
1433  (*(*tr).Hi()).g());
1434 
1435  (*tr).Emis().gf() =
1436  max( 1e-20f, (*tr).Emis().gf() );
1437 
1438  (*(*tr).Hi()).Pop() = LTE_pop * dense.xIonDense[nelem][nelem+1-ipISO] * dense.eden;
1439 
1440  (*tr).Emis().PopOpc() =
1441  (*(*tr).Lo()).Pop() -
1442  (*(*tr).Hi()).Pop() *
1443  (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
1444 
1445  (*tr).Emis().opacity() =
1446  (realnum)(abscf((*tr).Emis().gf(),
1447  (*tr).EnergyWN(),
1448  (*(*tr).Lo()).g()));
1449 
1450  /* a typical transition probability is of order 1e10 s-1 */
1451  double lifetime = 1e-10;
1452 
1453  (*tr).Emis().dampXvel() = (realnum)(
1454  (1.f/lifetime)/PI4/(*tr).EnergyWN());
1455  }
1456  }
1457  }
1458 
1459  return;
1460 }
1461 
1462 long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed )
1463 {
1464  DEBUG_ENTRY( "iso_get_total_num_levels()" );
1465 
1466  long tot_num_levels;
1467 
1468  /* return the number of levels up to and including nmaxResolved PLUS
1469  * the number (numCollapsed) of collapsed n-levels */
1470 
1471  if( ipISO == ipH_LIKE )
1472  {
1473  tot_num_levels = (long)( nmaxResolved * 0.5 *( nmaxResolved + 1 ) ) + numCollapsed;
1474  }
1475  else if( ipISO == ipHE_LIKE )
1476  {
1477  tot_num_levels = nmaxResolved*nmaxResolved + nmaxResolved + 1 + numCollapsed;
1478  }
1479  else
1480  TotalInsanity();
1481 
1482  return tot_num_levels;
1483 }
1484 
1485 void iso_update_num_levels( long ipISO, long nelem )
1486 {
1487  DEBUG_ENTRY( "iso_update_num_levels()" );
1488 
1489  /* This is the minimum resolved nmax. */
1490  ASSERT( iso_sp[ipISO][nelem].n_HighestResolved_max >= 3 );
1491 
1492  iso_sp[ipISO][nelem].numLevels_max =
1493  iso_get_total_num_levels( ipISO, iso_sp[ipISO][nelem].n_HighestResolved_max, iso_sp[ipISO][nelem].nCollapsed_max );
1494 
1495  if( iso_sp[ipISO][nelem].numLevels_max > iso_sp[ipISO][nelem].numLevels_malloc )
1496  {
1497  fprintf( ioQQQ, "The number of levels for ipISO %li, nelem %li, has been increased since the initial coreload.\n",
1498  ipISO, nelem );
1499  fprintf( ioQQQ, "This cannot be done.\n" );
1500  cdEXIT(EXIT_FAILURE);
1501  }
1502 
1503  /* set local copies to the max values */
1504  iso_sp[ipISO][nelem].numLevels_local = iso_sp[ipISO][nelem].numLevels_max;
1505  iso_sp[ipISO][nelem].nCollapsed_local = iso_sp[ipISO][nelem].nCollapsed_max;
1506  iso_sp[ipISO][nelem].n_HighestResolved_local = iso_sp[ipISO][nelem].n_HighestResolved_max;
1507 
1508  /* find the largest number of levels in any element in all iso sequences
1509  * we will allocate one matrix for ionization solver, and just use a piece of that memory
1510  * for smaller models. */
1511  max_num_levels = MAX2( max_num_levels, iso_sp[ipISO][nelem].numLevels_max);
1512 
1513  return;
1514 }
1515 
1516 void iso_collapsed_bnl_set( long ipISO, long nelem )
1517 {
1518 
1519  DEBUG_ENTRY( "iso_collapsed_bnl_set()" );
1520 
1521  double bnl_array[4][3][4][10] = {
1522  {
1523  /* H */
1524  {
1525  {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},
1526  {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},
1527  {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},
1528  {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}
1529  },
1530  {
1531  {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},
1532  {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},
1533  {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},
1534  {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}
1535  },
1536  {
1537  {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},
1538  {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},
1539  {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},
1540  {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}
1541  }
1542  },
1543  {
1544  /* He+ */
1545  {
1546  {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},
1547  {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},
1548  {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},
1549  {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}
1550  },
1551  {
1552  {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},
1553  {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},
1554  {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},
1555  {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}
1556  },
1557  {
1558  {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},
1559  {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},
1560  {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},
1561  {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}
1562  }
1563  },
1564  {
1565  /* He singlets */
1566  {
1567  {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},
1568  {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},
1569  {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},
1570  {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}
1571  },
1572  {
1573  {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},
1574  {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},
1575  {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},
1576  {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}
1577  },
1578  {
1579  {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},
1580  {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},
1581  {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},
1582  {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}
1583  }
1584  },
1585  {
1586  /* He triplets */
1587  {
1588  {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},
1589  {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},
1590  {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},
1591  {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}
1592  },
1593  {
1594  {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},
1595  {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},
1596  {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},
1597  {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}
1598  },
1599  {
1600  {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},
1601  {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},
1602  {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},
1603  {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}
1604  }
1605  }
1606  };
1607 
1608  double temps[4] = {5000., 10000., 15000., 20000. };
1609  double log_dens[3] = {2., 4., 6.};
1610  long ipTe, ipDens;
1611 
1612  ASSERT( nelem <= 1 );
1613 
1614  /* find temperature in tabulated values. */
1615  ipTe = hunt_bisect( temps, 4, phycon.te );
1616  ipDens = hunt_bisect( log_dens, 3, log10(dense.eden) );
1617 
1618  ASSERT( (ipTe >=0) && (ipTe < 3) );
1619  ASSERT( (ipDens >=0) && (ipDens < 2) );
1620 
1621  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++ )
1622  {
1623  for( long lHi=0; lHi<nHi; lHi++ )
1624  {
1625  for( long sHi=1; sHi<4; sHi++ )
1626  {
1627  if( ipISO == ipH_LIKE && sHi != 2 )
1628  continue;
1629  else if( ipISO == ipHE_LIKE && sHi != 1 && sHi != 3 )
1630  continue;
1631 
1632  double bnl_at_lo_den, bnl_at_hi_den, bnl;
1633  double bnl_max, bnl_min, temp, dens;
1634 
1635  long ipL = MIN2(9,lHi);
1636  long ip1;
1637 
1638  if( nelem==ipHYDROGEN )
1639  ip1 = 0;
1640  else if( nelem==ipHELIUM )
1641  {
1642  if( ipISO==ipH_LIKE )
1643  ip1 = 1;
1644  else if( ipISO==ipHE_LIKE )
1645  {
1646  if( sHi==1 )
1647  ip1 = 2;
1648  else if( sHi==3 )
1649  ip1 = 3;
1650  else
1651  TotalInsanity();
1652  }
1653  else
1654  TotalInsanity();
1655  }
1656  else
1657  TotalInsanity();
1658 
1659  temp = MAX2( temps[0], phycon.te );
1660  temp = MIN2( temps[3], temp );
1661 
1662  dens = MAX2( log_dens[0], log10(dense.eden) );
1663  dens = MIN2( log_dens[2], dens );
1664 
1665  /* Calculate the answer...must interpolate on two variables.
1666  * First interpolate on T, at both the lower and upper densities.
1667  * Then interpolate between these results for the right density. */
1668 
1669  if( temp < temps[0] && dens < log_dens[0] )
1670  bnl = bnl_array[ip1][0][0][ipL];
1671  else if( temp < temps[0] && dens >= log_dens[2] )
1672  bnl = bnl_array[ip1][2][0][ipL];
1673  else if( temp >= temps[3] && dens < log_dens[0] )
1674  bnl = bnl_array[ip1][0][3][ipL];
1675  else if( temp >= temps[3] && dens >= log_dens[2] )
1676  bnl = bnl_array[ip1][2][3][ipL];
1677  else
1678  {
1679  bnl_at_lo_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
1680  (bnl_array[ip1][ipDens][ipTe+1][ipL] - bnl_array[ip1][ipDens][ipTe][ipL]) + bnl_array[ip1][ipDens][ipTe][ipL];
1681 
1682  bnl_at_hi_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
1683  (bnl_array[ip1][ipDens+1][ipTe+1][ipL] - bnl_array[ip1][ipDens+1][ipTe][ipL]) + bnl_array[ip1][ipDens+1][ipTe][ipL];
1684 
1685  bnl = ( dens - log_dens[ipDens]) / (log_dens[ipDens+1] - log_dens[ipDens]) *
1686  (bnl_at_hi_den - bnl_at_lo_den) + bnl_at_lo_den;
1687  }
1688 
1690  {
1691  bnl_max = MAX4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
1692  bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
1693  ASSERT( bnl <= bnl_max );
1694 
1695  bnl_min = MIN4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
1696  bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
1697  ASSERT( bnl >= bnl_min );
1698  }
1699 
1700  iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] = bnl;
1701 
1702  ASSERT( iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] > 0. );
1703  }
1704  }
1705  }
1706 
1707  return;
1708 }
1709 
1710 
1711 void iso_collapsed_bnl_print( long ipISO, long nelem )
1712 {
1713  DEBUG_ENTRY( "iso_collapsed_bnl_print()" );
1714 
1715  for( long is = 1; is<=3; ++is)
1716  {
1717  if( ipISO == ipH_LIKE && is != 2 )
1718  continue;
1719  else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
1720  continue;
1721 
1722  char chSpin[3][9]= {"singlets", "doublets", "triplets"};
1723 
1724  /* give element number and spin */
1725  fprintf(ioQQQ," %s %s %s bnl\n",
1726  iso_ctrl.chISO[ipISO],
1727  elementnames.chElementSym[nelem],
1728  chSpin[is-1]);
1729 
1730  /* header with the l states */
1731  fprintf(ioQQQ," n\\l=> ");
1732  for( long i =0; i < iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++i)
1733  {
1734  fprintf(ioQQQ,"%2ld ",i);
1735  }
1736  fprintf(ioQQQ,"\n");
1737 
1738  /* loop over prin quant numbers, one per line, with l across */
1739  for( long in = 1; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in)
1740  {
1741  if( is==3 && in==1 )
1742  continue;
1743 
1744  fprintf(ioQQQ," %2ld ",in);
1745 
1746  for( long il = 0; il < in; ++il)
1747  {
1748  fprintf( ioQQQ, "%9.3e ", iso_sp[ipISO][nelem].bnl_effective[in][il][is] );
1749  }
1750  fprintf(ioQQQ,"\n");
1751  }
1752  }
1753 
1754  return;
1755 }
1756 
1757 void iso_collapsed_Aul_update( long ipISO, long nelem )
1758 {
1759  DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
1760 
1761  long ipFirstCollapsed = iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max;
1762 
1763  for( long ipLo=0; ipLo<ipFirstCollapsed; ipLo++ )
1764  {
1765  long spin = iso_sp[ipISO][nelem].st[ipLo].S();
1766 
1767  /* calculate effective Aul's from collapsed levels */
1768  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++ )
1769  {
1770  realnum Auls[2] = {
1771  iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][0],
1772  iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] };
1773 
1774  realnum EffectiveAul =
1775  Auls[0]*spin*(2.f*(L_(ipLo)+1.f)+1.f)*(realnum)iso_sp[ipISO][nelem].bnl_effective[nHi][ L_(ipLo)+1 ][spin];
1776 
1777  /* this is for n,L-1 -> n',L
1778  * make sure L-1 exists. */
1779  if( L_(ipLo) > 0 )
1780  {
1781  EffectiveAul +=
1782  Auls[1]*spin*(2.f*(L_(ipLo)-1.f)+1.f)*(realnum)iso_sp[ipISO][nelem].bnl_effective[nHi][ L_(ipLo)-1 ][spin];
1783  }
1784 
1785  if( ipISO==ipH_LIKE )
1786  EffectiveAul /= (2.f*nHi*nHi);
1787  else if( ipISO==ipHE_LIKE )
1788  EffectiveAul /= (4.f*nHi*nHi);
1789  else
1790  TotalInsanity();
1791 
1792  long ipHi = iso_sp[ipISO][nelem].QuantumNumbers2Index[nHi][ L_(ipLo)+1 ][spin];
1793 
1794  /* FINALLY, put the effective A in the proper Emis structure. */
1795  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = EffectiveAul;
1796 
1797  ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0. );
1798  }
1799  }
1800 
1801  return;
1802 }
1803 
1804 void iso_collapsed_lifetimes_update( long ipISO, long nelem )
1805 {
1806  DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
1807 
1808  for( long ipHi=iso_sp[ipISO][nelem].numLevels_max- iso_sp[ipISO][nelem].nCollapsed_max; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1809  {
1810  iso_sp[ipISO][nelem].st[ipHi].lifetime() = SMALLFLOAT;
1811 
1812  for( long ipLo=0; ipLo < ipHi; ipLo++ )
1813  {
1814  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1815  continue;
1816 
1817  iso_sp[ipISO][nelem].st[ipHi].lifetime() += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1818  }
1819 
1820  /* sum of A's was just stuffed, now invert for lifetime. */
1821  iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./iso_sp[ipISO][nelem].st[ipHi].lifetime();
1822 
1823  for( long ipLo=0; ipLo < ipHi; ipLo++ )
1824  {
1825  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
1826  continue;
1827 
1828  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1829  continue;
1830 
1831  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
1832  (1.f/iso_sp[ipISO][nelem].st[ipHi].lifetime())/
1833  PI4/iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
1834 
1835  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
1836  }
1837  }
1838 
1839  return;
1840 }
1841 
1842 #if 0
1843 STATIC void Prt_AGN_table( void )
1844 {
1845  /* the designation of the levels, chLevel[n][string] */
1846  multi_arr<char,2> chLevel(max_num_levels,10);
1847 
1848  /* create spectroscopic designation of labels */
1849  for( long ipLo=0; ipLo < iso_sp[ipISO][ipISO].numLevels_max-iso_sp[ipISO][ipISO].nCollapsed_max; ++ipLo )
1850  {
1851  long nelem = ipISO;
1852  sprintf( &chLevel.front(ipLo) , "%li %li%c", N_(ipLo), S_(ipLo), chL[MIN2(20,L_(ipLo))] );
1853  }
1854 
1855  /* option to print cs data for AGN */
1856  /* create spectroscopic designation of labels */
1857  {
1858  /* option to print particulars of some line when called */
1859  enum {AGN=false};
1860  if( AGN )
1861  {
1862 # define NTEMP 6
1863  double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
1864  double telog[NTEMP] ,
1865  cs ,
1866  ratecoef;
1867  long nelem = ipHELIUM;
1868  fprintf( ioQQQ,"trans");
1869  for( long i=0; i < NTEMP; ++i )
1870  {
1871  telog[i] = log10( te[i] );
1872  fprintf( ioQQQ,"\t%.3e",te[i]);
1873  }
1874  for( long i=0; i < NTEMP; ++i )
1875  {
1876  fprintf( ioQQQ,"\t%.3e",te[i]);
1877  }
1878  fprintf(ioQQQ,"\n");
1879 
1880  for( long ipHi=ipHe2s3S; ipHi< iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max; ++ipHi )
1881  {
1882  for( long ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
1883  {
1884 
1885  /* deltaN = 0 transitions may be wrong because
1886  * COLL_CONST below is only correct for electron colliders */
1887  if( N_(ipHi) == N_(ipLo) )
1888  continue;
1889 
1890  /* print the designations of the lower and upper levels */
1891  fprintf( ioQQQ,"%s - %s",
1892  &chLevel.front(ipLo) , &chLevel.front(ipHi) );
1893 
1894  /* print the interpolated collision strengths */
1895  for( long i=0; i < NTEMP; ++i )
1896  {
1897  phycon.alogte = telog[i];
1898  /* print cs */
1899  cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
1900  fprintf(ioQQQ,"\t%.2e", cs );
1901  }
1902 
1903  /* print the rate coefficients */
1904  for( long i=0; i < NTEMP; ++i )
1905  {
1906  phycon.alogte = telog[i];
1907  phycon.te = pow(10.,telog[i] );
1908  tfidle(false);
1909  cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
1910  /* collisional deexcitation rate */
1911  ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso_sp[ipHE_LIKE][nelem].st[ipLo].g() *
1912  sexp( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyK() / phycon.te );
1913  fprintf(ioQQQ,"\t%.2e", ratecoef );
1914  }
1915  fprintf(ioQQQ,"\n");
1916  }
1917  }
1918  cdEXIT(EXIT_FAILURE);
1919  }
1920  }
1921 
1922  return;
1923 }
1924 #endif
MIN2
#define MIN2
Definition: cddefines.h:765
MAX4
#define MAX4(a, b, c, d)
Definition: cddefines.h:796
iso_collapsed_bnl_print
void iso_collapsed_bnl_print(long ipISO, long nelem)
Definition: iso_create.cpp:1711
t_iso_sp::numLevels_malloc
long int numLevels_malloc
Definition: iso.h:502
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:403
EmissionProxy::opacity
realnum & opacity() const
Definition: emission.h:593
null_mole
molecule * null_mole
Definition: mole_species.cpp:64
t_iso_sp::st
qList st
Definition: iso.h:453
MIN4
#define MIN4(a, b, c, d)
Definition: cddefines.h:775
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
iso_assign_quantum_numbers
STATIC void iso_assign_quantum_numbers(void)
Definition: iso_create.cpp:697
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
multi_arr::alloc
void alloc()
Definition: container_classes.h:1116
opac
t_opac opac
Definition: opacity.cpp:5
ipSatelliteLines
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:37
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:567
iso_satellite
STATIC void iso_satellite(void)
Definition: iso_create.cpp:1325
abscf
double abscf(double gf, double enercm, double gl)
Definition: lines_service.cpp:122
H_BAR
const double H_BAR
Definition: physconst.h:144
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat.h:329
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:178
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
NISO
const int NISO
Definition: cddefines.h:265
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
t_isoCTRL::lgFSM
bool lgFSM[NISO]
Definition: iso.h:399
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
t_hydro::HLineWidth
realnum HLineWidth
Definition: hydrogenic.h:63
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
MAX2
#define MAX2
Definition: cddefines.h:786
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:182
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:487
iso_get_total_num_levels
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
Definition: iso_create.cpp:1462
t_isoCTRL::ipResoRedist
int ipResoRedist[NISO]
Definition: iso.h:374
phycon
t_phycon phycon
Definition: phycon.cpp:6
dense
t_dense dense
Definition: dense.cpp:24
helike_einsta.h
HeCollidSetup
void HeCollidSetup(void)
Definition: helike_cs.cpp:90
COLL_CONST
const double COLL_CONST
Definition: physconst.h:229
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:141
g
static double * g
Definition: species2.cpp:28
RYDLAM
const double RYDLAM
Definition: physconst.h:176
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_iso_sp::ipTrans
multi_arr< long, 2 > ipTrans
Definition: iso.h:448
t_isoCTRL::lgRandErrGen
bool lgRandErrGen[NISO]
Definition: iso.h:403
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
Transitions
vector< vector< TransitionList > > Transitions
Definition: taulines.cpp:33
molecule
Definition: mole.h:132
dsexp
double dsexp(double x)
Definition: service.cpp:953
SPEEDLIGHT
const double SPEEDLIGHT
Definition: physconst.h:100
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
IPRAD
#define IPRAD
Definition: iso.h:86
hunt_bisect
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:270
tfidle
STATIC void tfidle(bool lgForceUpdate)
Definition: temp_change.cpp:150
FINE_STRUCTURE
const double FINE_STRUCTURE
Definition: physconst.h:216
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:179
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
PI4
const double PI4
Definition: physconst.h:35
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
iso_collapsed_lifetimes_update
void iso_collapsed_lifetimes_update(long ipISO, long nelem)
Definition: iso_create.cpp:1804
max_num_levels
long int max_num_levels
Definition: iso.cpp:10
HeCSInterp
realnum HeCSInterp(long int nelem, long int ipHi, long int ipLo, long int Collider)
Definition: helike_cs.cpp:227
FillExtraLymanLine
STATIC void FillExtraLymanLine(const TransitionList::iterator &t, long ipISO, long nelem, long nHi)
Definition: iso_create.cpp:974
t_iso_sp::n_HighestResolved_local
long int n_HighestResolved_local
Definition: iso.h:507
nLine
static long int nLine
Definition: save_line.cpp:288
multi_arr::clone
const multi_geom< d, ALLOC > & clone() const
Definition: container_classes.h:1784
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:505
lines_service.h
t_iso_sp::BranchRatio
multi_arr< double, 2 > BranchRatio
Definition: iso.h:451
EVRYD
const double EVRYD
Definition: physconst.h:189
helike_energy
double helike_energy(long nelem, long ipLev)
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
iso_update_num_levels
void iso_update_num_levels(long ipISO, long nelem)
Definition: iso_create.cpp:1485
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
L_
#define L_(A_)
Definition: iso.h:21
chL
char chL[21]
Definition: iso_create.cpp:37
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
ERG1CM
const double ERG1CM
Definition: physconst.h:164
iso_collapsed_bnl_set
void iso_collapsed_bnl_set(long ipISO, long nelem)
Definition: iso_create.cpp:1516
POW2
#define POW2
Definition: cddefines.h:933
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:334
ipH1s
const int ipH1s
Definition: iso.h:27
multi_arr
Definition: container_classes.h:941
RYD_INF
const double RYD_INF
Definition: physconst.h:115
STATIC
#define STATIC
Definition: cddefines.h:101
WAVNRYD
const double WAVNRYD
Definition: physconst.h:173
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
iso_satellite_update
void iso_satellite_update(long nelem)
Definition: iso_create.cpp:1378
N_
#define N_(A_)
Definition: iso.h:20
t_isoCTRL::chISO
const char * chISO[NISO]
Definition: iso.h:330
TransitionProxy::ipEmis
int & ipEmis() const
Definition: transition.h:416
M
static const int M
Definition: thirdparty.cpp:2815
mole
t_mole_local mole
Definition: mole.cpp:7
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:813
multi_arr::invalidate
void invalidate()
Definition: container_classes.h:1057
t_iso_sp::tr
TransitionList * tr
Definition: iso.h:454
iso_zero
STATIC void iso_zero(void)
Definition: iso_create.cpp:410
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:161
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:365
realnum
float realnum
Definition: cddefines.h:107
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:355
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:144
ipLY_A
const int ipLY_A
Definition: cddefines.h:300
max
long max(int a, long b)
Definition: cddefines.h:779
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
HION_LTE_POP
const double HION_LTE_POP
Definition: physconst.h:157
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:438
molecule::index
int index
Definition: mole.h:169
S_
#define S_(A_)
Definition: iso.h:22
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1059
iso_collapsed_Aul_update
void iso_collapsed_Aul_update(long ipISO, long nelem)
Definition: iso_create.cpp:1757
ipH3s
const int ipH3s
Definition: iso.h:30
SatelliteLines
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:38
ipH3d
const int ipH3d
Definition: iso.h:32
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
t_iso_sp
Definition: iso.h:441
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
iso_allocate
STATIC void iso_allocate(void)
Definition: iso_create.cpp:454
t_isoCTRL::nLyaLevel
int nLyaLevel[NISO]
Definition: iso.h:377
hydro_bauman.h
ipH2p
const int ipH2p
Definition: iso.h:29
t_iso_sp::CascadeProb
multi_arr< double, 2 > CascadeProb
Definition: iso.h:450
t_iso_sp::ex
multi_arr< extra_tr, 2 > ex
Definition: iso.h:449
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:582
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:894
ipH2s
const int ipH2s
Definition: iso.h:28
t_isoCTRL::ipLyaRedist
int ipLyaRedist[NISO]
Definition: iso.h:374
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1080
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:101
t_iso_sp::bnl_effective
multi_arr< double, 3 > bnl_effective
Definition: iso.h:566
iso_recomb_setup
void iso_recomb_setup(long ipISO)
Definition: iso_radiative_recomb.cpp:838
ex
static double * ex
Definition: species2.cpp:28
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
iso_recomb_malloc
void iso_recomb_malloc(void)
Definition: iso_radiative_recomb.cpp:765
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
LIMELM
const int LIMELM
Definition: cddefines.h:262
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:688
ELECTRON_MASS
const double ELECTRON_MASS
Definition: physconst.h:91
iso_state_lifetime
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
Definition: iso_create.cpp:1054
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:310
iso_cascade
void iso_cascade(long ipISO, long nelem)
Definition: iso_create.cpp:1098
t_phycon::te_ryd
double te_ryd
Definition: phycon.h:17
iso_create
void iso_create(void)
Definition: iso_create.cpp:39
t_dense::eden
double eden
Definition: dense.h:190
t_isoCTRL::ipSubRedist
int ipSubRedist[NISO]
Definition: iso.h:374
t_isoCTRL::stat_ion
realnum stat_ion[NISO]
Definition: iso.h:362
t_phycon::alogte
double alogte
Definition: phycon.h:82
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1132
t_phycon::te
double te
Definition: phycon.h:11
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:309
fixit
void fixit(void)
Definition: service.cpp:991
AnonStates
qList AnonStates(1)
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
POW4
#define POW4
Definition: cddefines.h:947
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
ATOMIC_MASS_UNIT
const double ATOMIC_MASS_UNIT
Definition: physconst.h:88
ipH3p
const int ipH3p
Definition: iso.h:31
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
t_isoCTRL::nLyman_malloc
long int nLyman_malloc[NISO]
Definition: iso.h:334
t_phycon::te32
double te32
Definition: phycon.h:48
ipCRDW
const int ipCRDW
Definition: cddefines.h:298
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:664
iso_recomb_auxiliary_free
void iso_recomb_auxiliary_free(void)
Definition: iso_radiative_recomb.cpp:825
hydro_transprob
realnum hydro_transprob(long nelem, long ipHi, long ipLo)
Definition: hydroeinsta.cpp:45
ipDOUBLET
Definition: iso.h:82
Generated on Sat Oct 19 2013 10:43:11 for cloudy by   doxygen 1.8.5