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-2022 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_adfa.h"
8 #include "dense.h"
9 #include "helike.h"
10 #include "helike_einsta.h"
11 #include "hydro_bauman.h"
12 #include "hydrogenic.h"
13 #include "hydroeinsta.h"
14 #include "iso.h"
15 #include "opacity.h"
16 #include "phycon.h"
17 #include "taulines.h"
18 #include "mole.h"
19 #include "freebound.h"
20 #include "lines_service.h"
21 #include "prt.h"
22 
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 static vector<species> isoSpecies( NISO );
40 
41 
48 void iso_setRedisFun (long ipISO, long nelem, long ipLo, long ipHi)
49 {
50  if( ipLo == 0 && ipHi == iso_ctrl.nLyaLevel[ipISO] )
51  {
52  long redis = iso_ctrl.ipLyaRedist[ipISO];
53  // H LyA has a special redistribution function
54  if( ipISO==ipH_LIKE && nelem==ipHYDROGEN )
55  redis = ipLY_A;
56  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = redis;
57  }
58  else if( ipLo == 0 )
59  {
60  /* these are rest of Lyman lines,
61  * complete redistribution, doppler core only, K2 core, default ipCRD */
62  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
63  }
64  else
65  {
66  /* all lines coming from excited states, default is complete
67  * redis with wings, ipCRDW*/
68  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = iso_ctrl.ipSubRedist[ipISO];
69  }
70 
71  return;
72 }
73 
74 
75 
82 void iso_setOpacity (long ipISO, long nelem, long ipLo, long ipHi)
83 {
84  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA ||
85  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0.)
86  {
87  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() = 0.;
88  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() = 0.;
89  }
90  else
91  {
92  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() =
93  (realnum)(GetGF(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(),
94  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
95  iso_sp[ipISO][nelem].st[ipHi].g()));
96  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() > 0.);
97 
98  /* derive the abs coef, call to function is gf, wl (A), g_low */
99  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() =
100  (realnum)(abscf(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf(),
101  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
102  iso_sp[ipISO][nelem].st[ipLo].g()));
103  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() > 0.);
104  }
105 
106  return;
107 }
108 
109 
110 
111 void iso_create(void)
112 {
113  long int ipHi,
114  ipLo;
115 
116  static int nCalled = 0;
117 
118  double HIonPoten;
119 
120  DEBUG_ENTRY( "iso_create()" );
121 
122  /* > 1 if not first call, then just zero arrays out */
123  if( nCalled > 0 )
124  {
125  iso_zero();
126  return;
127  }
128 
129  /* this is first call, increment the nCalled counterso never do this again */
130  ++nCalled;
131 
132  /* these are the statistical weights of the ions */
133  iso_ctrl.stat_ion[ipH_LIKE] = 1.f;
134  iso_ctrl.stat_ion[ipHE_LIKE] = 2.f;
135 
136  /* this routine allocates all the memory
137  * needed for the iso sequence structures */
138  iso_allocate();
139 
140  /* loop over iso sequences and assign quantum numbers to all levels */
141  iso_assign_quantum_numbers();
142 
143  /* this is a dummy line, junk it too. */
144  (*TauDummy).Junk();
145  (*TauDummy).AddHiState();
146  (*TauDummy).AddLoState();
147  (*TauDummy).AddLine2Stack();
148 
149  /********************************************/
150  /********** Line and level energies ********/
151  /********************************************/
152  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
153  {
154  /* main hydrogenic arrays, fill with sane values */
155  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
156  {
157  /* must always do helium even if turned off */
158  if( nelem < 2 || dense.lgElmtOn[nelem] )
159  {
160  /* the level energies and line wavelengths are derived from these
161  * Rydberg constants. Changing the constants will change the wavelength,
162  * but beware that many H I, He I, and He II wavelengths are hardwired into
163  * the code. These include lists of Case B line predictions and places
164  * where lines are retrieved. If the Rydberg constants are changed it would
165  * be wise to invest the time to derive all wavelengths self consistently
166  * with global variables and remove the explicit lists. See the source
167  * files changed in r11678 for an idea of where these lines occur.
168  */
169  if( nelem==ipHYDROGEN )
170  /* this is the NIST value for H itself */
171  HIonPoten = 0.99946650834637;
172  else if( nelem==ipHELIUM )
173  HIonPoten = 3.9996319917;
174  else
175  /* Dima's data in phfit.dat have ionization potentials in eV
176  * with four significant figures*/
177  HIonPoten = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD;
178  ASSERT(HIonPoten > 0.);
179 
180  double EnergyRydGround = 0.;
181  /* go from ground to the highest level */
182  for( ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
183  {
184  double EnergyWN, EnergyRyd;
185 
186  if( ipISO == ipH_LIKE )
187  {
188  EnergyRyd = HIonPoten/POW2((double)N_(ipHi));
189  if( ipHi==ipH2p && nelem >4 )
190  {
191  /* only Ka for Z>5 */
192  double Z = nelem+1;
193  /* >>20 sept 24, Priyanka Chakraborty's energy correction in wavenumbers
194  *updated phfit.dat according to NIST values for H-like ions
195  *Compared Ritz formula with NIST and
196  *generated correction to the polynomial */
197  double E_correction = 0.1783*pow(Z,4)-1.8313*pow(Z,3)+27.803*pow(Z,2)-208.04*Z+570.59;
198  /* convert to Ryd and apply */
199  EnergyRyd -= E_correction * WAVNRYD;
200  }
201  }
202  else if( ipISO == ipHE_LIKE )
203  {
204  EnergyRyd = helike_energy( nelem, ipHi ) * WAVNRYD;
205  }
206  else
207  {
208  /* Other iso sequences don't exist yet. */
209  TotalInsanity();
210  }
211 
212  /* >>chng 02 feb 09, change test to >= 0 since we now use 0 for 2s-2p */
213  ASSERT(EnergyRyd >= 0.);
214 
215  iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd = EnergyRyd;
216  if (ipHi == 0)
217  EnergyRydGround = EnergyRyd;
218  iso_sp[ipISO][nelem].st[ipHi].energy().set(EnergyRydGround-EnergyRyd);
219 
220  /* now loop from ground to level below ipHi */
221  for( ipLo=0; ipLo < ipHi; ipLo++ )
222  {
223  EnergyWN = RYD_INF * (iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd -
224  iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
225 
226  /* This is the minimum line energy we will allow. */
227  /* \todo 2 wire this to lowest energy of code. */
228  if( EnergyWN==0 && ipISO==ipHE_LIKE )
229  EnergyWN = 0.0001;
230 
231  if( EnergyWN < 0. )
232  EnergyWN = -1.0 * EnergyWN;
233 
234  /* transition energy in various units: */
235  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() = (realnum)EnergyWN;
236 
237  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() >= 0.);
238  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() >= 0.);
239  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyK() >= 0.);
240 
242  if( N_(ipLo)==N_(ipHi) && ipISO==ipH_LIKE )
243  {
244  iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() = 0.;
245  }
246  else
247  {
248  /* make following an air wavelength */
249  iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() =
250  (realnum)(1.0e8/
251  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()/
252  RefIndex( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()));
253  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() > 0.);
254  }
255  }
256  }
257 
258  /* fill the extra Lyman lines */
259  for( ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
260  {
261  FillExtraLymanLine( ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi], ipISO, nelem, ipHi );
262  }
263  }
264  }
265  }
266 
267  /***************************************************************/
268  /***** Set up recombination tables for later interpolation *****/
269  /***************************************************************/
270  /* NB - the above is all we need if we are compiling recombination tables. */
271  iso_recomb_malloc();
272  iso_recomb_setup( ipH_LIKE );
273  iso_recomb_setup( ipHE_LIKE );
274  iso_recomb_auxiliary_free();
275 
276  /* set up helium collision tables */
277  HeCollidSetup();
278 
279  /***********************************************************************************/
280  /********** Transition Probabilities, Redistribution Functions, Opacitites ********/
281  /***********************************************************************************/
282  enum { DEBUG_EINA_A_NNP = false };
283 
284  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
285  {
286  if( ipISO == ipH_LIKE )
287  {
288  /* do nothing here */
289  }
290  else if( ipISO == ipHE_LIKE )
291  {
292  /* This routine reads in transition probabilities from a file. */
293  HelikeTransProbSetup();
294  }
295  else
296  {
297  TotalInsanity();
298  }
299 
300  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
301  {
302  /* must always do helium even if turned off */
303  if( nelem < 2 || dense.lgElmtOn[nelem] )
304  {
305  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
306  {
307  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
308  {
309  realnum Aul;
310 
311  /* transition prob, EinstA uses current H atom indices */
312  if( ipISO == ipH_LIKE )
313  {
314  Aul = hydro_transprob( nelem, ipHi, ipLo );
315  }
316  else if( ipISO == ipHE_LIKE )
317  {
318  Aul = helike_transprob(nelem, ipHi, ipLo);
319  }
320  else
321  {
322  TotalInsanity();
323  }
324 
325  if( Aul <= iso_ctrl.SmallA )
326  iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipEmis() = -1;
327  else
328  iso_sp[ipISO][nelem].trans(ipHi,ipLo).AddLine2Stack();
329 
330  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = Aul;
331 
332  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0.);
333 
334  iso_setRedisFun(ipISO, nelem, ipLo, ipHi);
335 
336  iso_setOpacity(ipISO, nelem, ipLo, ipHi);
337 
338  if( DEBUG_EINA_A_NNP )
339  {
340  printf( "%ld\t %ld\t %ld\t %.4e\n",
341  nelem+1,
342  N_( ipHi ), N_( ipLo ),
343  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() );
344  }
345  }
346  }
347  }
348  }
349  }
350 
351  /************************************************/
352  /********** Fine Structure Mixing - FSM ********/
353  /************************************************/
354  if( iso_ctrl.lgFSM[ipHE_LIKE] )
355  {
356  /* set some special optical depth values */
357  for( long nelem=ipHE_LIKE; nelem < LIMELM; nelem++ )
358  {
359  /* must always do helium even if turned off */
360  if( nelem < 2 || dense.lgElmtOn[nelem] )
361  {
362  for( ipHi=ipHe2s3S; ipHi<iso_sp[ipHE_LIKE][nelem].numLevels_max; ipHi++ )
363  {
364  for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
365  {
366  DoFSMixing( nelem, ipLo, ipHi );
367  }
368  }
369  }
370  }
371  }
372 
373  /****************************************************/
374  /********** lifetimes and damping constants ********/
375  /****************************************************/
376  enum { DEBUG_LIFETIMES = false };
377 
378  for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
379  {
380  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
381  {
382  /* define these for H and He always */
383  if( nelem < 2 || dense.lgElmtOn[nelem] )
384  {
385  /* these are not defined and must never be used */
386  iso_sp[ipISO][nelem].st[0].lifetime() = -FLT_MAX;
387 
388  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
389  {
390  iso_sp[ipISO][nelem].st[ipHi].lifetime() = SMALLFLOAT;
391 
392  for( ipLo=0; ipLo < ipHi; ipLo++ )
393  {
394  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
395  continue;
396 
397  iso_sp[ipISO][nelem].st[ipHi].lifetime() += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
398  }
399 
400  /* sum of A's was just stuffed, now invert for lifetime. */
401  iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./iso_sp[ipISO][nelem].st[ipHi].lifetime();
402 
403  if( DEBUG_LIFETIMES )
404  {
405  printf( "%ld\t %ld\t %.4e\n",
406  nelem+1,
407  N_( ipHi ),
408  iso_sp[ipISO][nelem].st[ipHi].lifetime() );
409  }
410 
411  for( ipLo=0; ipLo < ipHi; ipLo++ )
412  {
413  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
414  continue;
415 
416  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
417  continue;
418 
419  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
420  (1.f/iso_sp[ipISO][nelem].st[ipHi].lifetime())/
421  PI4/iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
422 
423  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
424  }
425  }
426  }
427  }
428  }
429 
430  /* zero out some line information */
431  iso_zero();
432 
433  /* loop over iso sequences */
434  for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
435  {
436  for( long nelem = ipISO; nelem < LIMELM; nelem++ )
437  {
438  /* must always do helium even if turned off */
439  if( nelem == ipISO || dense.lgElmtOn[nelem] )
440  {
441  /* calculate cascade probabilities, branching ratios, and associated errors. */
442  iso_cascade( ipISO, nelem);
443  }
444  }
445  }
446 
447  iso_satellite();
448 
449  for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
450  iso_satellite_update( nelem );
451 
452  /***************************************/
453  /********** Stark Broadening **********/
454  /***************************************/
455  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
456  {
457  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
458  {
459  if( nelem < 2 || dense.lgElmtOn[nelem] )
460  {
461  for( long ipHi= 1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
462  {
463  for( long ipLo=0; ipLo < ipHi; ++ipLo )
464  {
465  iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk = 0.;
466  iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up = 0.;
467  }
468  }
469  }
470  }
471  }
472 
473  // arrays set up, do not allow number of levels to change in later sims
474  lgHydroMalloc = true;
475 
476  return;
477 }
478 
479 
480 /* ============================================================================== */
481 STATIC void iso_zero(void)
482 {
483  DEBUG_ENTRY( "iso_zero()" );
484 
485  hydro.HLineWidth = 0.;
486 
487  /****************************************************/
488  /********** initialize some variables **********/
489  /****************************************************/
490  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
491  {
492  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
493  {
494  if( nelem < 2 || dense.lgElmtOn[nelem] )
495  {
496  for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
497  {
498  iso_sp[ipISO][nelem].st[ipHi].Pop() = 0.;
499  iso_sp[ipISO][nelem].fb[ipHi].Reset();
500  }
501  if (ipISO <= nelem)
502  iso_sp[ipISO][nelem].st[0].Pop() =
503  dense.xIonDense[nelem][nelem-ipISO];
504  }
505 
506  if( nelem < 2 )
507  {
508  iso_collapsed_Aul_update( ipISO, nelem );
509  iso_collapsed_lifetimes_update( ipISO, nelem );
510  }
511  }
512  }
513 
514  /* ground state of H and He is different since totally determine
515  * their own opacities */
516  iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ConOpacRatio = 1e-5;
517  iso_sp[ipH_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
518  iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
519 
520  return;
521 }
522 
523 STATIC void iso_allocate(void)
524 {
525 
526  DEBUG_ENTRY( "iso_allocate()" );
527 
528  /* the hydrogen and helium like iso-sequences */
529  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
530  {
531  isoSpecies[ipISO].database = iso_ctrl.chISO[ipISO];
532 
533  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
534  {
535  /* only grab core for elements that are turned on */
536  if( nelem < 2 || dense.lgElmtOn[nelem] )
537  {
538  t_iso_sp *sp = &iso_sp[ipISO][nelem];
539  sp->numLevels_malloc = sp->numLevels_max;
540 
541  {
542  sp->lgPrtMatrix = false;
543  string chemicalLabel = makeChemical( nelem, nelem-ipISO );
544  if( strcmp( chemicalLabel.c_str(), prt.matrix.species ) == 0 )
545  sp->lgPrtMatrix = true;
546  }
547 
548  ASSERT( sp->numLevels_max > 0 );
549  ASSERT( iso_ctrl.nLyman_malloc[ipISO] == iso_ctrl.nLyman[ipISO] );
550 
551  sp->CachedAs.reserve( MAX2(1, sp->nCollapsed_max) );
552 
553  sp->ipTrans.reserve( sp->numLevels_max );
554  sp->ex.reserve( sp->numLevels_max );
555  sp->CascadeProb.reserve( sp->numLevels_max );
556  sp->BranchRatio.reserve( sp->numLevels_max );
557  //sp->st.resize( sp->numLevels_max );
558  sp->fb.resize( sp->numLevels_max );
559 
560  for( long i = 0; i < sp->nCollapsed_max; ++i )
561  {
562  // NB - this is the total number of levels we would have if all n were resolved.
563  long numLevels = iso_get_total_num_levels( ipISO, sp->n_HighestResolved_max + sp->nCollapsed_max, 0 );
564  sp->CachedAs.reserve( i, numLevels );
565  for( long i1 = 0; i1 < numLevels; ++i1 )
566  {
567  /* allocate two spaces delta L +/- 1 */
568  sp->CachedAs.reserve( i, i1, 2 );
569  }
570  }
571 
572  sp->QuantumNumbers2Index.reserve( sp->n_HighestResolved_max + sp->nCollapsed_max + 1 );
573 
574  for( long i = 1; i <= sp->n_HighestResolved_max + sp->nCollapsed_max; ++i )
575  {
576  /* Allocate proper number of angular momentum quantum number. */
577  sp->QuantumNumbers2Index.reserve( i, i );
578 
579  for( long i1 = 0; i1 < i; ++i1 )
580  {
581  /* This may have to change for other iso sequences. */
582  ASSERT( NISO == 2 );
583  /* Allocate 4 spaces for multiplicity. H-like will be accessed with "2" for doublet,
584  * He-like will be accessed via "1" for singlet or "3" for triplet. "0" will not be used. */
585  sp->QuantumNumbers2Index.reserve( i, i1, 4 );
586  }
587  }
588 
589  for( long n=1; n < sp->numLevels_max; ++n )
590  {
591  sp->ipTrans.reserve( n, n );
592  }
593 
594  for( long n=0; n < sp->numLevels_max; ++n )
595  {
596  sp->ex.reserve( n, sp->numLevels_max );
597  sp->CascadeProb.reserve( n, sp->numLevels_max );
598  sp->BranchRatio.reserve( n, sp->numLevels_max );
599  }
600 
601  sp->ipTrans.alloc();
602  sp->ex.alloc();
603  sp->CascadeProb.alloc();
604  sp->BranchRatio.alloc();
605 
606  sp->CachedAs.alloc();
607  sp->QuantumNumbers2Index.alloc();
608  sp->IndexIfAllResolved.alloc( sp->QuantumNumbers2Index.clone() );
609  sp->QuantumNumbers2Index.invalidate();
610  sp->IndexIfAllResolved.invalidate();
611  }
612  }
613  }
614 
615  ipSatelliteLines.reserve( NISO );
616  ipExtraLymanLines.reserve( NISO );
617 
618  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
619  {
620  ipSatelliteLines.reserve( ipISO, LIMELM );
621  ipExtraLymanLines.reserve( ipISO, LIMELM );
622 
623  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
624  {
625  /* only grab core for elements that are turned on */
626  if( nelem < 2 || dense.lgElmtOn[nelem] )
627  {
628  ASSERT( iso_sp[ipISO][nelem].numLevels_max > 0 );
629 
630  ipSatelliteLines.reserve( ipISO, nelem, iso_sp[ipISO][nelem].numLevels_max );
631  ipExtraLymanLines.reserve( ipISO, nelem, iso_ctrl.nLyman_malloc[ipISO] );
632  }
633  }
634  }
635 
636  ipSatelliteLines.alloc();
637  ipExtraLymanLines.alloc();
638 
639  Transitions.resize(NISO);
640  SatelliteLines.resize(NISO);
641  ExtraLymanLines.resize(NISO);
642  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
643  {
644  Transitions[ipISO].reserve(LIMELM);
645  SatelliteLines[ipISO].reserve(LIMELM);
646  ExtraLymanLines[ipISO].reserve(LIMELM);
647  for( long nelem=0; nelem < ipISO; ++nelem )
648  {
649  Transitions[ipISO].push_back(
650  TransitionList("Insanity",&AnonStates));
651  SatelliteLines[ipISO].push_back(
652  TransitionList("Insanity",&AnonStates));
653  ExtraLymanLines[ipISO].push_back(
654  TransitionList("Insanity",&AnonStates));
655  }
656  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
657  {
658  if( nelem < 2 || dense.lgElmtOn[nelem] )
659  {
660  Transitions[ipISO].push_back(
661  TransitionList("Isosequence",&iso_sp[ipISO][nelem].st));
662  SatelliteLines[ipISO].push_back(
663  TransitionList("SatelliteLines",&iso_sp[ipISO][nelem].st));
664  ExtraLymanLines[ipISO].push_back(
665  TransitionList("ExtraLymanLines",&iso_sp[ipISO][nelem].st));
666  }
667  else
668  {
669  Transitions[ipISO].push_back(
670  TransitionList("Insanity",&AnonStates));
671  SatelliteLines[ipISO].push_back(
672  TransitionList("Insanity",&AnonStates));
673  ExtraLymanLines[ipISO].push_back(
674  TransitionList("Insanity",&AnonStates));
675  }
676  }
677  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
678  {
679  /* only grab core for elements that are turned on */
680  if( nelem < 2 || dense.lgElmtOn[nelem] )
681  {
682  if( iso_ctrl.lgDielRecom[ipISO] )
683  {
684  SatelliteLines[ipISO][nelem].resize( iso_sp[ipISO][nelem].numLevels_max );
685  AllTransitions.push_back(SatelliteLines[ipISO][nelem]);
686  unsigned int nLine = 0;
687  for( long ipLo=0; ipLo<iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
688  {
689  /* Upper level is continuum, use a generic state
690  * lower level is the same as the index. */
691  ipSatelliteLines[ipISO][nelem][ipLo] = nLine;
692  SatelliteLines[ipISO][nelem][nLine].Junk();
693  long ipHi = iso_sp[ipISO][nelem].numLevels_max;
694  SatelliteLines[ipISO][nelem][nLine].setHi(ipHi);
695  SatelliteLines[ipISO][nelem][nLine].setLo(ipLo);
696  SatelliteLines[ipISO][nelem][nLine].AddLine2Stack();
697  ++nLine;
698  }
699  ASSERT(SatelliteLines[ipISO][nelem].size() == nLine);
700  }
701 
702  //iso_sp[ipISO][nelem].tr.resize( iso_sp[ipISO][nelem].ipTrans.size() );
703  //iso_sp[ipISO][nelem].tr.states() = &iso_sp[ipISO][nelem].st;
704  Transitions[ipISO][nelem].resize( iso_sp[ipISO][nelem].ipTrans.size() );
705  AllTransitions.push_back(Transitions[ipISO][nelem]);
706  unsigned int nTransition=0;
707  for( long ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
708  {
709  for( long ipLo=0; ipLo < ipHi; ipLo++ )
710  {
711  /* set ENTIRE array to impossible values, in case of bad pointer */
712  iso_sp[ipISO][nelem].ipTrans[ipHi][ipLo] = nTransition;
713  Transitions[ipISO][nelem][nTransition].Junk();
714  Transitions[ipISO][nelem][nTransition].setHi(ipHi);
715  Transitions[ipISO][nelem][nTransition].setLo(ipLo);
716  ++nTransition;
717  }
718  }
719  ASSERT(Transitions[ipISO][nelem].size() == nTransition);
720  iso_sp[ipISO][nelem].tr = &Transitions[ipISO][nelem];
721 
722  /* junk the extra Lyman lines */
723  AllTransitions.push_back(ExtraLymanLines[ipISO][nelem]);
724  ExtraLymanLines[ipISO][nelem].resize(iso_ctrl.nLyman_malloc[ipISO]-2);
725  ExtraLymanLines[ipISO][nelem].states() = &iso_sp[ipISO][nelem].st;
726  unsigned int nExtraLyman = 0;
727  for( long ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
728  {
729  ipExtraLymanLines[ipISO][nelem][ipHi] = nExtraLyman;
730  ExtraLymanLines[ipISO][nelem][nExtraLyman].Junk();
731  long ipHi_offset = iso_sp[ipISO][nelem].numLevels_max + ipHi - 2;
732  if( iso_ctrl.lgDielRecom[ipISO] )
733  ipHi_offset += 1;
734  ExtraLymanLines[ipISO][nelem][nExtraLyman].setHi(ipHi_offset);
735  /* lower level is just ground state of the ion */
736  ExtraLymanLines[ipISO][nelem][nExtraLyman].setLo(0);
737  ExtraLymanLines[ipISO][nelem][nExtraLyman].AddLine2Stack();
738  ++nExtraLyman;
739  }
740  ASSERT(ExtraLymanLines[ipISO][nelem].size() == nExtraLyman);
741  }
742  }
743  }
744 
745  // associate line and level stacks with species
746  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
747  {
748  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
749  {
750  if( dense.lgElmtOn[nelem] )
751  {
752  long ion = nelem - ipISO;
753  ASSERT( ion >= 0 && ion <= nelem );
754 
755  molecule *spmole = findspecies( makeChemical( nelem, ion ).c_str() );
756  ASSERT( spmole != null_mole );
757 
758  mole.species[ spmole->index ].dbase = &isoSpecies[ipISO];
759  mole.species[ spmole->index ].dbase->numLevels_local = iso_sp[ipISO][nelem].numLevels_local;
760  mole.species[ spmole->index ].dbase->numLevels_max = iso_sp[ipISO][nelem].numLevels_max;
761  mole.species[ spmole->index ].levels = &iso_sp[ipISO][nelem].st;
762  mole.species[ spmole->index ].lines = &Transitions[ipISO][nelem];
763  }
764  }
765  }
766 
767  return;
768 }
769 
770 STATIC void iso_assign_quantum_numbers(void)
771 {
772  long int
773  ipLo,
774  level,
775  i,
776  in,
777  il,
778  is,
779  ij;
780 
781  DEBUG_ENTRY( "iso_assign_quantum_numbers()" );
782 
783  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
784  {
785  long ipISO = ipH_LIKE;
786  /* only check elements that are turned on */
787  if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
788  {
789  i = 0;
790 
791  /* 2 for doublet */
792  is = ipDOUBLET;
793 
794  /* this loop is over quantum number n */
795  for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
796  {
797  for( il = 0L; il < in; ++il )
798  {
799  iso_sp[ipISO][nelem].st[i].n() = in;
800  iso_sp[ipISO][nelem].st[i].S() = is;
801  iso_sp[ipISO][nelem].st[i].l() = il;
802  iso_sp[ipISO][nelem].st[i].j() = -1;
803  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
804  iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] = i;
805  ++i;
806  }
807  }
808  /* Now find indices of levels if all were resolved */
809  {
810  long i2 = i;
811  for( in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
812  in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
813  {
814  for( il = 0L; il < in; ++il )
815  {
816  iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] = i2;
817  ++i2;
818  }
819  }
820  }
821  /* now do the collapsed levels */
822  in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
823  for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
824  {
825  iso_sp[ipISO][nelem].st[level].n() = in;
826  iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
827  iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
828  iso_sp[ipISO][nelem].st[level].j() = -1;
829  /* Point every l to same index for collapsed levels. */
830  for( il = 0; il < in; ++il )
831  {
832  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
833  }
834  ++in;
835  }
836  --in;
837 
838  /* confirm that we did not overrun the array */
839  ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
840 
841  /* confirm that n is positive and not greater than the max n. */
842  ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
843 
844  /* Verify states and QuantumNumbers2Index agree in all cases */
845  for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
846  {
847  for( il = 0L; il < in; ++il )
848  {
849  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
850  if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
851  {
852  /* Must only check these for resolved levels...
853  * collapsed levels have pointers for l and s that will blow if used. */
854  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
855  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
856  // QuantumNumbers2Index and IndexIfAllResolved must agree for resolved levels
857  ASSERT( iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] == iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] );
858  }
859  }
860  }
861  }
862  }
863 
864  /* then do he-like */
865  for( long nelem=ipHELIUM; nelem < LIMELM; nelem++ )
866  {
867  long ipISO = ipHE_LIKE;
868  /* only check elements that are turned on */
869  if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
870  {
871  i = 0;
872 
873  /* this loop is over quantum number n */
874  for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
875  {
876  for( il = 0L; il < in; ++il )
877  {
878  for( is = 3L; is >= 1L; is -= 2 )
879  {
880  /* All levels except singlet P follow the ordering scheme: */
881  /* lower l's have lower energy */
882  /* triplets have lower energy */
883  if( (il == 1L) && (is == 1L) )
884  continue;
885  /* n = 1 has no triplet, of course. */
886  if( (in == 1L) && (is == 3L) )
887  continue;
888 
889  /* singlets */
890  if( is == 1 )
891  {
892  iso_sp[ipISO][nelem].st[i].n() = in;
893  iso_sp[ipISO][nelem].st[i].S() = is;
894  iso_sp[ipISO][nelem].st[i].l() = il;
895  /* this is not a typo, J=L for singlets. */
896  iso_sp[ipISO][nelem].st[i].j() = il;
897  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
898  iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] = i;
899  ++i;
900  }
901  /* 2 triplet P is j-resolved */
902  else if( (in == 2) && (il == 1) && (is == 3) )
903  {
904  ij = 0;
905  do
906  {
907  iso_sp[ipISO][nelem].st[i].n() = in;
908  iso_sp[ipISO][nelem].st[i].S() = is;
909  iso_sp[ipISO][nelem].st[i].l() = il;
910  iso_sp[ipISO][nelem].st[i].j() = ij;
911  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
912  iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] = i;
913  ++i;
914  ++ij;
915  /* repeat this for the separate j-levels within 2^3P. */
916  } while ( ij < 3 );
917  }
918  else
919  {
920  iso_sp[ipISO][nelem].st[i].n() = in;
921  iso_sp[ipISO][nelem].st[i].S() = is;
922  iso_sp[ipISO][nelem].st[i].l() = il;
923  iso_sp[ipISO][nelem].st[i].j() = -1L;
924  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
925  iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] = i;
926  ++i;
927  }
928  }
929  }
930  /* Insert singlet P at the end of every sequence for a given n. */
931  if( in > 1L )
932  {
933  iso_sp[ipISO][nelem].st[i].n() = in;
934  iso_sp[ipISO][nelem].st[i].S() = 1L;
935  iso_sp[ipISO][nelem].st[i].l() = 1L;
936  iso_sp[ipISO][nelem].st[i].j() = 1L;
937  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][1][1] = i;
938  iso_sp[ipISO][nelem].IndexIfAllResolved[in][1][1] = i;
939  ++i;
940  }
941  }
942  /* Now find indices of levels if all were resolved */
943  {
944  long i2 = i;
945  for( in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
946  in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
947  {
948  ASSERT( in > 2L );
949  for( il = 0L; il < in; ++il )
950  {
951  for( is = 3L; is >= 1L; is -= 2 )
952  {
953  /* All levels except singlet P follow the ordering scheme: */
954  /* lower l's have lower energy */
955  /* triplets have lower energy */
956  if( (il == 1L) && (is == 1L) )
957  continue;
958 
959  iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] = i2;
960  ++i2;
961  }
962  }
963  /* Insert singlet P at the end of every sequence for a given n. */
964  if( in > 1L )
965  {
966  iso_sp[ipISO][nelem].IndexIfAllResolved[in][1][1] = i2;
967  ++i2;
968  }
969  }
970  }
971  /* now do the collapsed levels */
972  in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
973  for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
974  {
975  iso_sp[ipISO][nelem].st[level].n() = in;
976  iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
977  iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
978  iso_sp[ipISO][nelem].st[level].j() = -1;
979  /* Point every l and s to same index for collapsed levels. */
980  for( il = 0; il < in; ++il )
981  {
982  for( is = 1; is <= 3; is += 2 )
983  {
984  iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
985  }
986  }
987  ++in;
988  }
989  --in;
990 
991  /* confirm that we did not overrun the array */
992  ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
993 
994  /* confirm that n is positive and not greater than the max n. */
995  ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
996 
997  /* Verify states and QuantumNumbers2Index agree in all cases */
998  for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
999  {
1000  for( il = 0L; il < in; ++il )
1001  {
1002  for( is = 3L; is >= 1; is -= 2 )
1003  {
1004  /* Ground state is not triplicate. */
1005  if( (in == 1L) && (is == 3L) )
1006  continue;
1007 
1008  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
1009  if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
1010  {
1011  /* Must only check these for resolved levels...
1012  * collapsed levels have pointers for l and s that will blow if used. */
1013  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
1014  ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
1015  // QuantumNumbers2Index and IndexIfAllResolved must agree for resolved levels
1016  ASSERT( iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] == iso_sp[ipISO][nelem].IndexIfAllResolved[in][il][is] );
1017  }
1018  }
1019  }
1020  }
1021  }
1022  }
1023 
1024  for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
1025  {
1026  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
1027  {
1028  /* must always do helium even if turned off */
1029  if( nelem < 2 || dense.lgElmtOn[nelem] )
1030  {
1031  for( ipLo=ipH1s; ipLo < iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
1032  {
1033  iso_sp[ipISO][nelem].st[ipLo].nelem() = (int)(nelem+1);
1034  iso_sp[ipISO][nelem].st[ipLo].IonStg() = (int)(nelem+1-ipISO);
1035 
1036  if( iso_sp[ipISO][nelem].st[ipLo].j() >= 0 )
1037  {
1038  iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*iso_sp[ipISO][nelem].st[ipLo].j()+1.f;
1039  }
1040  else if( iso_sp[ipISO][nelem].st[ipLo].l() >= 0 )
1041  {
1042  iso_sp[ipISO][nelem].st[ipLo].g() = (2.f*iso_sp[ipISO][nelem].st[ipLo].l()+1.f) *
1043  iso_sp[ipISO][nelem].st[ipLo].S();
1044  }
1045  else
1046  {
1047  if( ipISO == ipH_LIKE )
1048  iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
1049  else if( ipISO == ipHE_LIKE )
1050  iso_sp[ipISO][nelem].st[ipLo].g() = 4.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
1051  else
1052  {
1053  /* replace this with correct thing if more sequences are added. */
1054  TotalInsanity();
1055  }
1056  }
1057  char chConfiguration[32];
1058  long nCharactersWritten = 0;
1059 
1060  ASSERT( iso_sp[ipISO][nelem].st[ipLo].n() < 1000 );
1061 
1062  /* Treat H-like levels as collapsed, for the purposes of
1063  * reporting comments in the output of 'save lines labels'.
1064  * For the rest, include J only if defined, and for singlets if positive. */
1065  if( ( ipISO == ipH_LIKE && ipLo > ipH2p ) ||
1066  iso_sp[ipISO][nelem].st[ipLo].n() > iso_sp[ipISO][nelem].n_HighestResolved_max )
1067  {
1068  nCharactersWritten = sprintf( chConfiguration, "n=%3li",
1069  iso_sp[ipISO][nelem].st[ipLo].n() );
1070  }
1071  else if( ( iso_sp[ipISO][nelem].st[ipLo].j() > 0 &&
1072  iso_sp[ipISO][nelem].st[ipLo].S() == ipSINGLET ) ||
1073  ( iso_sp[ipISO][nelem].st[ipLo].j() >= 0 &&
1074  iso_sp[ipISO][nelem].st[ipLo].S() == ipTRIPLET ) )
1075  {
1076  nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c_%li",
1077  iso_sp[ipISO][nelem].st[ipLo].n(),
1078  iso_sp[ipISO][nelem].st[ipLo].S(),
1079  chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l() ) ],
1080  iso_sp[ipISO][nelem].st[ipLo].j() );
1081  }
1082  else
1083  {
1084  nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c",
1085  iso_sp[ipISO][nelem].st[ipLo].n(),
1086  iso_sp[ipISO][nelem].st[ipLo].S(),
1087  chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l()) ] );
1088  }
1089 
1090  if( nCharactersWritten < 0 || nCharactersWritten > 31 )
1091  TotalInsanity();
1092 
1093  iso_sp[ipISO][nelem].st[ipLo].chConfig() = chConfiguration;
1094  }
1095  }
1096  }
1097  }
1098  return;
1099 }
1100 
1101 #if defined(__ICC) && defined(__i386)
1102 #pragma optimization_level 1
1103 #endif
1104 STATIC void FillExtraLymanLine( const TransitionList::iterator& t, long ipISO, long nelem, long nHi )
1105 {
1106  double Enerwn, Aul;
1107 
1108  DEBUG_ENTRY( "FillExtraLymanLine()" );
1109 
1110  (*(*t).Hi()).status() = LEVEL_INACTIVE;
1111 
1112  /* atomic number or charge and stage: */
1113  (*(*t).Hi()).nelem() = (int)(nelem+1);
1114  (*(*t).Hi()).IonStg() = (int)(nelem+1-ipISO);
1115 
1116  (*(*t).Hi()).n() = nHi;
1117 
1118  /* statistical weight is same as statistical weight of corresponding LyA. */
1119  (*(*t).Hi()).g() = iso_sp[ipISO][nelem].st[iso_ctrl.nLyaLevel[ipISO]].g();
1120 
1121  /* \todo add correct configuration, or better still link to standard level */
1122  (*(*t).Hi()).chConfig() = "ExtraLyman level (probably duplicate)";
1123 
1124  /* energies */
1125  Enerwn = iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd * RYD_INF * ( 1. - 1./POW2((double)nHi) );
1126 
1127  /* transition energy in various units:*/
1128  (*t).EnergyWN() = (realnum)(Enerwn);
1129  (*t).WLAng() = (realnum)(1.0e8/ Enerwn/ RefIndex(Enerwn));
1130  (*(*t).Hi()).energy().set( Enerwn, "cm^-1" );
1131 
1132  if( ipISO == ipH_LIKE )
1133  {
1134  Aul = H_Einstein_A( nHi, 1, 1, 0, nelem+1 );
1135  }
1136  else
1137  {
1138  if( nelem == ipHELIUM )
1139  {
1140  /* A simple fit for the calculation of Helium lyman Aul's. */
1141  Aul = (1.508e10) / pow((double)nHi,2.975);
1142  }
1143  else
1144  {
1145  /* Fit to values given in
1146  * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., &
1147  * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */
1148  /* originally astro.ph. 0201454 */
1149  Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
1150  }
1151  }
1152 
1153  (*t).Emis().Aul() = (realnum)Aul;
1154 
1155  (*(*t).Hi()).lifetime() = iso_state_lifetime( ipISO, nelem, nHi, 1 );
1156 
1157  (*t).Emis().dampXvel() = (realnum)( 1.f / (*(*t).Hi()).lifetime() / PI4 / (*t).EnergyWN() );
1158 
1159  (*t).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
1160 
1161  (*t).Emis().gf() = (realnum)(GetGF((*t).Emis().Aul(), (*t).EnergyWN(), (*(*t).Hi()).g()));
1162 
1163  /* derive the abs coef, call to function is Emis().gf(), wl (A), g_low */
1164  (*t).Emis().opacity() = (realnum)(abscf((*t).Emis().gf(), (*t).EnergyWN(), (*(*t).Lo()).g()));
1165 
1166  /* create array indices that will blow up */
1167  (*t).ipCont() = INT_MIN;
1168  (*t).Emis().ipFine() = INT_MIN;
1169 
1170  {
1171  /* option to print particulars of some line when called
1172  * a prettier print statement is near where chSpin is defined below
1173  * search for "pretty print" */
1174  enum {DEBUG_LOC=false};
1175  if( DEBUG_LOC )
1176  {
1177  fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
1178  nelem+1,
1179  nHi,
1180  (*t).Emis().Aul() ,
1181  (*t).Emis().opacity()
1182  );
1183  }
1184  }
1185  return;
1186 }
1187 
1188 /* calculate radiative lifetime of an individual iso state */
1189 double iso_state_lifetime( long ipISO, long nelem, long n, long l )
1190 {
1191  /* >>refer hydro lifetimes Horbatsch, M. W., Horbatsch, M. and Hessels, E. A. 2005, JPhysB, 38, 1765 */
1192 
1193  double tau, t0, eps2;
1194  /* mass of electron */
1195  double m = ELECTRON_MASS;
1196  /* nuclear mass */
1197  double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT;
1198  double mu = (m*M)/(M+m);
1199  long z = 1;
1200  long Z = nelem + 1 - ipISO;
1201 
1202  DEBUG_ENTRY( "iso_state_lifetime()" );
1203 
1204  /* this should not be used for l=0 per the Horbatsch et al. paper */
1205  ASSERT( l > 0 );
1206 
1207  eps2 = 1. - ( l*l + l + 8./47. - (l+1.)/69./n ) / POW2( (double)n );
1208 
1209  t0 = 3. * H_BAR * powi( (double)n, 5 ) /
1210  ( 2. * pow4( (double)( z * Z ) ) * powi( FINE_STRUCTURE, 5 ) * mu * pow2( SPEEDLIGHT ) ) *
1211  pow2( (m + M)/(Z*m + z*M) );
1212 
1213  tau = t0 * ( 1. - eps2 ) /
1214  ( 1. + 19./88.*( (1./eps2 - 1.) * log( 1. - eps2 ) + 1. -
1215  0.5 * eps2 - 0.025 * eps2 * eps2 ) );
1216 
1217  if( ipISO == ipHE_LIKE )
1218  {
1219  /* iso_state_lifetime is not spin specific, must exclude helike triplet here. */
1220  tau /= 3.;
1221  /* this is also necessary to correct the helike lifetimes */
1222  tau *= 1.1722 * pow( (double)nelem, 0.1 );
1223  }
1224 
1225  /* would probably need a new lifetime algorithm for any other iso sequences. */
1226  ASSERT( ipISO <= ipHE_LIKE );
1227  ASSERT( tau > 0. );
1228 
1229  return tau;
1230 }
1231 
1232 /* calculate cascade probabilities, branching ratios, and associated errors. */
1233 void iso_cascade( long ipISO, long nelem )
1234 {
1235  /* The sum of all A's coming out of a given n,
1236  * Below we assert a monotonic trend. */
1237  double *SumAPerN;
1238 
1239  long int i, j, ipLo, ipHi;
1240 
1241  DEBUG_ENTRY( "iso_cascade()" );
1242 
1243  SumAPerN = ((double*)MALLOC( (size_t)(iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double )));
1244  memset( SumAPerN, 0, (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double ) );
1245 
1246  /* Initialize some ground state stuff, easier here than in loops. */
1247  iso_sp[ipISO][nelem].CascadeProb[0][0] = 1.;
1248  if( iso_ctrl.lgRandErrGen[ipISO] )
1249  {
1250  iso_sp[ipISO][nelem].fb[0].SigmaAtot = 0.;
1251  iso_sp[ipISO][nelem].ex[0][0].SigmaCascadeProb = 0.;
1252  }
1253 
1254  /***************************************************************************/
1255  /****** Cascade probabilities, Branching ratios, and associated errors *****/
1256  /***************************************************************************/
1257  for( ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1258  {
1259  double SumAs = 0.;
1260 
1266  /* initialize variables. */
1267  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipHi] = 1.;
1268  iso_sp[ipISO][nelem].CascadeProb[ipHi][0] = 0.;
1269  iso_sp[ipISO][nelem].BranchRatio[ipHi][0] = 0.;
1270 
1271  if( iso_ctrl.lgRandErrGen[ipISO] )
1272  {
1273  iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = 0.;
1274  iso_sp[ipISO][nelem].ex[ipHi][ipHi].SigmaCascadeProb = 0.;
1275  }
1276 
1277  long ipLoStart = 0;
1278  if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) )
1279  ipLoStart = 1;
1280 
1281  for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1282  {
1283  SumAs += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1284  }
1285 
1286  for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1287  {
1288  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1289  {
1290  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1291  iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] = 0.;
1292  continue;
1293  }
1294 
1295  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1296  iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] =
1297  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() / SumAs;
1298 
1299  ASSERT( iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] <= 1.0000001 );
1300 
1301  SumAPerN[N_(ipHi)] += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1302 
1303  /* there are some negative energy transitions, where the order
1304  * has changed, but these are not optically allowed, these are
1305  * same n, different L, forbidden transitions */
1306  ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() > 0. ||
1307  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA );
1308 
1309  if( iso_ctrl.lgRandErrGen[ipISO] )
1310  {
1311  ASSERT( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] >= 0. );
1312  /* Uncertainties in A's are added in quadrature, square root is taken below. */
1313  iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot +=
1314  pow2( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] *
1315  (double)iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() );
1316  }
1317  }
1318 
1319  if( iso_ctrl.lgRandErrGen[ipISO] )
1320  {
1321  /* Uncertainties in A's are added in quadrature above, square root taken here. */
1322  iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = sqrt( iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot );
1323  }
1324 
1325  /* cascade probabilities */
1326  for( i=0; i<ipHi; i++ )
1327  {
1328  for( ipLo=0; ipLo<=i; ipLo++ )
1329  {
1330  iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] += iso_sp[ipISO][nelem].BranchRatio[ipHi][i] * iso_sp[ipISO][nelem].CascadeProb[i][ipLo];
1331  }
1332  }
1333 
1334  if( iso_ctrl.lgRandErrGen[ipISO] )
1335  {
1336  for( ipLo=0; ipLo<ipHi; ipLo++ )
1337  {
1338  double SigmaCul = 0.;
1339  for( i=ipLo; i<ipHi; i++ )
1340  {
1341  if( iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul() > iso_ctrl.SmallA )
1342  {
1343  /* Uncertainties in A's and cascade probabilities */
1344  double SigmaA = iso_sp[ipISO][nelem].ex[ipHi][i].Error[IPRAD] *
1345  iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul();
1346  SigmaCul +=
1347  pow2(SigmaA*iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime()) +
1348  pow2(iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*iso_sp[ipISO][nelem].BranchRatio[ipHi][i]*
1349  iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime()) +
1350  pow2(iso_sp[ipISO][nelem].ex[i][ipLo].SigmaCascadeProb*iso_sp[ipISO][nelem].BranchRatio[ipHi][i]);
1351  }
1352  }
1353  SigmaCul = sqrt(SigmaCul);
1354  iso_sp[ipISO][nelem].ex[ipHi][ipLo].SigmaCascadeProb = SigmaCul;
1355  }
1356  }
1357  }
1358 
1359  /************************************************************************/
1360  /*** Allowed decay conversion probabilities. See Robbins68b, Table 1. ***/
1361  /************************************************************************/
1362  {
1363  enum {DEBUG_LOC=false};
1364 
1365  if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) )
1366  {
1367  /* To output Bm(n,l; ipLo), set ipLo, hi_l, and hi_s accordingly. */
1368  long int hi_l,hi_s;
1369  double Bm;
1370 
1371  /* these must be set for following output to make sense
1372  * as is, a dangerous bit of code - set NaN for safety */
1373  hi_s = -100000;
1374  hi_l = -100000;
1375  ipLo = -100000;
1376  /* tripS to 2^3P */
1377  //hi_l = 0, hi_s = 3, ipLo = ipHe2p3P0;
1378 
1379  /* tripD to 2^3P */
1380  //hi_l = 2, hi_s = 3, ipLo = ipHe2p3P0;
1381 
1382  /* tripP to 2^3S */
1383  //hi_l = 1, hi_s = 3, ipLo = ipHe2s3S;
1384 
1385  ASSERT( hi_l != iso_sp[ipISO][nelem].st[ipLo].l() );
1386 
1387  fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
1388  fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");
1389 
1390  for( ipHi=ipHe2p3P2; ipHi<iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipHi++ )
1391  {
1392  /* Pick out excitations from metastable 2tripS to ntripP. */
1393  if( (iso_sp[ipISO][nelem].st[ipHi].l() == 1) && (iso_sp[ipISO][nelem].st[ipHi].S() == 3) )
1394  {
1395  fprintf(ioQQQ,"\n%ld\t",iso_sp[ipISO][nelem].st[ipHi].n());
1396  j = 0;
1397  Bm = 0;
1398  for( i = ipLo; i<=ipHi; i++)
1399  {
1400  if( (iso_sp[ipISO][nelem].st[i].l() == hi_l) && (iso_sp[ipISO][nelem].st[i].S() == hi_s) )
1401  {
1402  if( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
1403  {
1404  Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * ( iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P0] +
1405  iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P1] + iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P2] );
1406  }
1407  else
1408  Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * iso_sp[ipISO][nelem].BranchRatio[i][ipLo];
1409 
1410  if( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
1411  {
1412  j++;
1413  if(j == 3)
1414  {
1415  fprintf(ioQQQ,"%2.4e\t",Bm);
1416  Bm = 0;
1417  }
1418  }
1419  else
1420  {
1421  fprintf(ioQQQ,"%2.4e\t",Bm);
1422  Bm = 0;
1423  }
1424  }
1425  }
1426  }
1427  }
1428  fprintf(ioQQQ,"\n\n");
1429  }
1430  }
1431 
1432  /******************************************************/
1433  /*** Lifetimes should increase monotonically with ***/
1434  /*** increasing n...Make sure the A's decrease. ***/
1435  /******************************************************/
1436  for( i=2; i < iso_sp[ipISO][nelem].n_HighestResolved_max; ++i)
1437  {
1438  ASSERT( (SumAPerN[i] > SumAPerN[i+1]) || opac.lgCaseB );
1439  }
1440 
1441  {
1442  enum {DEBUG_LOC=false};
1443  if( DEBUG_LOC /* && (ipISO == ipH_LIKE) && (nelem == ipHYDROGEN) */)
1444  {
1445  for( i = 2; i<= (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max); ++i)
1446  {
1447  fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]);
1448  }
1449  }
1450  }
1451 
1452  free( SumAPerN );
1453 
1454  return;
1455 }
1456 
1458 /* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002 */
1459 /* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001. */
1460 STATIC void iso_satellite( void )
1461 {
1462  DEBUG_ENTRY( "iso_satellite()" );
1463 
1464  for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
1465  {
1466  for( long nelem = ipISO; nelem < LIMELM; nelem++ )
1467  {
1468  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
1469  {
1470  for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
1471  {
1472  TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
1473  (*tr).Zero();
1474 
1475  /* Make approximation that all levels have energy of H-like 2s level */
1476  /* Lines to 1s2s have roughly energy of parent Ly-alpha. So lines to 1snL will have an energy
1477  * smaller by the difference between nL and 2s energies. Therefore, the following has
1478  * energy of parent Ly-alpha MINUS the difference between daughter level and daughter n=2 level. */
1479  (*tr).WLAng() = (realnum)(RYDLAM/
1480  ((iso_sp[ipISO-1][nelem].fb[0].xIsoLevNIonRyd - iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd) -
1481  (iso_sp[ipISO][nelem].fb[1].xIsoLevNIonRyd- iso_sp[ipISO][nelem].fb[i].xIsoLevNIonRyd)) );
1482 
1483  (*tr).EnergyWN() = 1.e8f / (*tr).WLAng();
1484 
1485  (*tr).Emis().iRedisFun() = ipCRDW;
1486  /* this is not the usual nelem, is it atomic not C scale. */
1487  (*(*tr).Hi()).nelem() = nelem + 1;
1488  (*(*tr).Hi()).IonStg() = nelem + 1 - ipISO;
1489  fixit("what should the stat weight of the upper level be? For now say 2.");
1490  (*(*tr).Hi()).g() = 2.f;
1491 
1492  /* \todo add correct configuration, or better still link to standard level */
1493  (*(*tr).Hi()).chConfig() = "Satellite level (probably duplicate)";
1494 
1495  // The lower level must already be initialized.
1496  ASSERT( (*(*tr).Lo()).g() == iso_sp[ipISO][nelem].st[i].g() );
1497  //(*(*tr).Lo()).g() = iso_sp[ipISO][nelem].st[i].g();
1498  (*tr).Emis().PopOpc() =
1499  (*(*tr).Lo()).Pop();
1500 
1501  (*tr).Emis().pump() = 0.;
1502 
1503  }
1504  }
1505  }
1506  }
1507 
1508  return;
1509 }
1510 
1511 void iso_satellite_update( long nelem )
1512 {
1513  double ConBoltz, LTE_pop=SMALLFLOAT+FLT_EPSILON, factor1, ConvLTEPOP;
1514 
1515  DEBUG_ENTRY( "iso_satellite_update()" );
1516 
1517  for( long ipISO = ipHE_LIKE; ipISO < MIN2(NISO,nelem+1); ipISO++ )
1518  {
1519  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
1520  {
1521  /* This Boltzmann factor is exp( +ioniz energy / Te ). For simplicity, we make
1522  * the fair approximation that all of the autoionizing levels have an energy
1523  * equal to the parents n=2. */
1524  ConBoltz = dsexp(iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd/phycon.te_ryd);
1525 
1526  for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
1527  {
1528  double dr_rate = iso_sp[ipISO][nelem].fb[i].DielecRecomb * iso_ctrl.lgDielRecom[ipISO];
1529 
1530  TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
1531  (*tr).Emis().xObsIntensity() =
1532  (*tr).Emis().xIntensity() =
1533  dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO] *
1534  ERG1CM * (*tr).EnergyWN();
1535 
1536  /* We set line intensity above using a rate, but here we need a transition probability.
1537  * We can obtain this by dividing dr_rate by the population of the autoionizing level.
1538  * We assume this level is in statistical equilibrium. */
1539  factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
1540  (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
1541 
1542  /* term in () is stat weight of electron * ion */
1543  ConvLTEPOP = powpq(factor1,3,2)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;
1544 
1545  if( ConBoltz >= SMALLDOUBLE )
1546  {
1547  /* The energy used to calculate ConBoltz above
1548  * should be negative since this is above the continuum, but
1549  * to be safe we calculate ConBoltz with a positive energy above
1550  * and multiply by it here instead of dividing. */
1551  LTE_pop = (*(*tr).Hi()).g() * ConBoltz * ConvLTEPOP;
1552  }
1553 
1554  LTE_pop = max( LTE_pop, 1e-30f );
1555 
1556  /* Now the transition probability is simply dr_rate/LTE_pop. */
1557  (*tr).Emis().Aul() = (realnum)(dr_rate/LTE_pop);
1558  (*tr).Emis().Aul() =
1559  max( iso_ctrl.SmallA, (*tr).Emis().Aul() );
1560 
1561  (*tr).Emis().gf() = (realnum)GetGF(
1562  (*tr).Emis().Aul(),
1563  (*tr).EnergyWN(),
1564  (*(*tr).Hi()).g());
1565 
1566  (*tr).Emis().gf() =
1567  max( 1e-20f, (*tr).Emis().gf() );
1568 
1569  (*(*tr).Hi()).Pop() = LTE_pop * dense.xIonDense[nelem][nelem+1-ipISO] * dense.eden;
1570  // In the approximation used here, DepartCoef is unity by definition.
1571  (*(*tr).Hi()).DepartCoef() = 1.;
1572 
1573  (*tr).Emis().PopOpc() =
1574  (*(*tr).Lo()).Pop() -
1575  (*(*tr).Hi()).Pop() *
1576  (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
1577 
1578  (*tr).Emis().opacity() =
1579  (realnum)(abscf((*tr).Emis().gf(),
1580  (*tr).EnergyWN(),
1581  (*(*tr).Lo()).g()));
1582 
1583  /* a typical transition probability is of order 1e10 s-1 */
1584  double lifetime = 1e-10;
1585 
1586  (*tr).Emis().dampXvel() = (realnum)(
1587  (1.f/lifetime)/PI4/(*tr).EnergyWN());
1588  }
1589  }
1590  }
1591 
1592  return;
1593 }
1594 
1595 long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed )
1596 {
1597  DEBUG_ENTRY( "iso_get_total_num_levels()" );
1598 
1599  long tot_num_levels;
1600 
1601  /* return the number of levels up to and including nmaxResolved PLUS
1602  * the number (numCollapsed) of collapsed n-levels */
1603 
1604  if( ipISO == ipH_LIKE )
1605  {
1606  tot_num_levels = (long)( nmaxResolved * 0.5 *( nmaxResolved + 1 ) ) + numCollapsed;
1607  }
1608  else if( ipISO == ipHE_LIKE )
1609  {
1610  tot_num_levels = nmaxResolved*nmaxResolved + nmaxResolved + 1 + numCollapsed;
1611  }
1612  else
1613  TotalInsanity();
1614 
1615  return tot_num_levels;
1616 }
1617 
1618 void iso_update_num_levels( long ipISO, long nelem )
1619 {
1620  DEBUG_ENTRY( "iso_update_num_levels()" );
1621 
1622  /* This is the minimum resolved nmax. */
1623  ASSERT( iso_sp[ipISO][nelem].n_HighestResolved_max >= 3 );
1624 
1625  iso_sp[ipISO][nelem].numLevels_max =
1626  iso_get_total_num_levels( ipISO, iso_sp[ipISO][nelem].n_HighestResolved_max, iso_sp[ipISO][nelem].nCollapsed_max );
1627 
1628  if( iso_sp[ipISO][nelem].numLevels_max > iso_sp[ipISO][nelem].numLevels_malloc )
1629  {
1630  fprintf( ioQQQ, "The number of levels for ipISO %li, nelem %li, has been increased since the initial coreload.\n",
1631  ipISO, nelem );
1632  fprintf( ioQQQ, "This cannot be done.\n" );
1633  cdEXIT(EXIT_FAILURE);
1634  }
1635 
1636  /* set local copies to the max values */
1637  iso_sp[ipISO][nelem].numLevels_local = iso_sp[ipISO][nelem].numLevels_max;
1638  iso_sp[ipISO][nelem].nCollapsed_local = iso_sp[ipISO][nelem].nCollapsed_max;
1639  iso_sp[ipISO][nelem].n_HighestResolved_local = iso_sp[ipISO][nelem].n_HighestResolved_max;
1640 
1641  /* find the largest number of levels in any element in all iso sequences
1642  * we will allocate one matrix for ionization solver, and just use a piece of that memory
1643  * for smaller models. */
1644  max_num_levels = MAX2( max_num_levels, iso_sp[ipISO][nelem].numLevels_max);
1645 
1646  return;
1647 }
1648 
1649 void iso_collapsed_Aul_update( long ipISO, long nelem )
1650 {
1651  DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
1652 
1653  long ipFirstCollapsed = iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max;
1654 
1655  for( long ipLo=0; ipLo<ipFirstCollapsed; ipLo++ )
1656  {
1657  long spin = iso_sp[ipISO][nelem].st[ipLo].S();
1658 
1659  /* calculate effective Aul's from collapsed levels */
1660  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++ )
1661  {
1662  // NB - 2nd dim of CachedAs is indexed properly by IndexIfAllResolved, but ipLo is always resolved here, so there is no difference.
1663  realnum Auls[2] = {
1664  iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][0],
1665  iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] };
1666 
1667  realnum EffectiveAul =
1668  Auls[0]*spin*(2.f*(L_(ipLo)+1.f)+1.f);
1669 
1670  /* this is for n,L-1 -> n',L
1671  * make sure L-1 exists. */
1672  if( L_(ipLo) > 0 )
1673  {
1674  EffectiveAul +=
1675  Auls[1]*spin*(2.f*(L_(ipLo)-1.f)+1.f);
1676  }
1677 
1678  if( ipISO==ipH_LIKE )
1679  EffectiveAul /= (2.f*nHi*nHi);
1680  else if( ipISO==ipHE_LIKE )
1681  EffectiveAul /= (4.f*nHi*nHi);
1682  else
1683  TotalInsanity();
1684 
1685  long ipHi = iso_sp[ipISO][nelem].QuantumNumbers2Index[nHi][ L_(ipLo)+1 ][spin];
1686 
1687  /* FINALLY, put the effective A in the proper Emis structure. */
1688  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = EffectiveAul;
1689 
1690  ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0. );
1691  }
1692  }
1693 
1694  return;
1695 }
1696 
1697 void iso_collapsed_lifetimes_update( long ipISO, long nelem )
1698 {
1699  DEBUG_ENTRY( "iso_collapsed_lifetimes_update()" );
1700 
1701  for( long ipHi=iso_sp[ipISO][nelem].numLevels_max- iso_sp[ipISO][nelem].nCollapsed_max; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1702  {
1703  double totalAul = SMALLFLOAT;
1704 
1705  for( long ipLo=0; ipLo < ipHi; ipLo++ )
1706  {
1707  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1708  continue;
1709 
1710  totalAul += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1711  }
1712 
1713  /* sum of A's was just stuffed, now invert for lifetime. */
1714  iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./totalAul;
1715 
1716  for( long ipLo=0; ipLo < ipHi; ipLo++ )
1717  {
1718  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
1719  continue;
1720 
1721  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
1722  continue;
1723 
1724  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
1725  totalAul/(PI4*iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()));
1726 
1727  ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
1728  }
1729  }
1730 
1731  return;
1732 }
1733 
1734 #if 0
1735 STATIC void Prt_AGN_table( void )
1736 {
1737  /* the designation of the levels, chLevel[n][string] */
1738  multi_arr<char,2> chLevel(max_num_levels,10);
1739 
1740  /* create spectroscopic designation of labels */
1741  for( long ipLo=0; ipLo < iso_sp[ipISO][ipISO].numLevels_max-iso_sp[ipISO][ipISO].nCollapsed_max; ++ipLo )
1742  {
1743  long nelem = ipISO;
1744  sprintf( &chLevel.front(ipLo) , "%li %li%c", N_(ipLo), S_(ipLo), chL[MIN2(20,L_(ipLo))] );
1745  }
1746 
1747  /* option to print cs data for AGN */
1748  /* create spectroscopic designation of labels */
1749  {
1750  /* option to print particulars of some line when called */
1751  enum {AGN=false};
1752  if( AGN )
1753  {
1754 # define NTEMP 6
1755  double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
1756  double telog[NTEMP] ,
1757  cs ,
1758  ratecoef;
1759  long nelem = ipHELIUM;
1760  fprintf( ioQQQ,"trans");
1761  for( long i=0; i < NTEMP; ++i )
1762  {
1763  telog[i] = log10( te[i] );
1764  fprintf( ioQQQ,"\t%.3e",te[i]);
1765  }
1766  for( long i=0; i < NTEMP; ++i )
1767  {
1768  fprintf( ioQQQ,"\t%.3e",te[i]);
1769  }
1770  fprintf(ioQQQ,"\n");
1771 
1772  for( long ipHi=ipHe2s3S; ipHi< iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max; ++ipHi )
1773  {
1774  for( long ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
1775  {
1776 
1777  /* deltaN = 0 transitions may be wrong because
1778  * COLL_CONST below is only correct for electron colliders */
1779  if( N_(ipHi) == N_(ipLo) )
1780  continue;
1781 
1782  /* print the designations of the lower and upper levels */
1783  fprintf( ioQQQ,"%s - %s",
1784  &chLevel.front(ipLo) , &chLevel.front(ipHi) );
1785 
1786  /* print the interpolated collision strengths */
1787  for( long i=0; i < NTEMP; ++i )
1788  {
1789  phycon.alogte = telog[i];
1790  /* print cs */
1791  cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
1792  fprintf(ioQQQ,"\t%.2e", cs );
1793  }
1794 
1795  /* print the rate coefficients */
1796  for( long i=0; i < NTEMP; ++i )
1797  {
1798  phycon.alogte = telog[i];
1799  phycon.te = exp10(telog[i] );
1800  tfidle(false);
1801  cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
1802  /* collisional deexcitation rate */
1803  ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso_sp[ipHE_LIKE][nelem].st[ipLo].g() *
1804  sexp( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyK() / phycon.te );
1805  fprintf(ioQQQ,"\t%.2e", ratecoef );
1806  }
1807  fprintf(ioQQQ,"\n");
1808  }
1809  }
1810  cdEXIT(EXIT_FAILURE);
1811  }
1812  }
1813 
1814  return;
1815 }
1816 #endif
pow4
T pow4(T a)
Definition: cddefines.h:995
t_iso_sp::numLevels_malloc
long int numLevels_malloc
Definition: iso.h:533
atmdat_adfa.h
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:450
EmissionProxy::opacity
realnum & opacity() const
Definition: emission.h:650
null_mole
molecule * null_mole
Definition: mole_species.cpp:50
t_iso_sp::st
qList st
Definition: iso.h:482
exp10
double exp10(double x)
Definition: cddefines.h:1368
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:65
iso_assign_quantum_numbers
STATIC void iso_assign_quantum_numbers(void)
Definition: iso_create.cpp:770
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
multi_arr::alloc
void alloc()
Definition: container_classes.h:1113
opac
t_opac opac
Definition: opacity.cpp:5
ipSatelliteLines
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:34
t_iso_sp::CachedAs
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:596
iso_satellite
STATIC void iso_satellite(void)
Definition: iso_create.cpp:1460
abscf
double abscf(double gf, double enercm, double gl)
Definition: lines_service.cpp:113
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat_adfa.h:61
lgHydroMalloc
bool lgHydroMalloc
Definition: cdinit.cpp:32
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:246
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
NISO
const int NISO
Definition: cddefines.h:311
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:49
t_isoCTRL::lgFSM
bool lgFSM[NISO]
Definition: iso.h:426
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:48
t_hydro::HLineWidth
realnum HLineWidth
Definition: hydrogenic.h:100
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:46
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:250
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:518
iso_get_total_num_levels
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
Definition: iso_create.cpp:1595
t_isoCTRL::ipResoRedist
int ipResoRedist[NISO]
Definition: iso.h:394
phycon
t_phycon phycon
Definition: phycon.cpp:6
helike_einsta.h
HeCollidSetup
void HeCollidSetup(void)
Definition: helike_cs.cpp:241
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:999
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:167
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_iso_sp::ipTrans
multi_arr< long, 2 > ipTrans
Definition: iso.h:477
t_isoCTRL::lgRandErrGen
bool lgRandErrGen[NISO]
Definition: iso.h:430
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:43
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:481
Transitions
vector< vector< TransitionList > > Transitions
Definition: taulines.cpp:30
molecule
Definition: mole.h:142
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
dsexp
double dsexp(double x)
Definition: service.cpp:1038
HeCSInterp
realnum HeCSInterp(long nelem, long ipHi, long ipLo, long Collider)
Definition: helike_cs.cpp:390
IPRAD
#define IPRAD
Definition: iso.h:88
t_iso_sp::IndexIfAllResolved
multi_arr< long, 3 > IndexIfAllResolved
Definition: iso.h:492
tfidle
STATIC void tfidle(bool lgForceUpdate)
Definition: temp_change.cpp:151
dense
t_dense dense
Definition: global.cpp:15
ipTRIPLET
Definition: iso.h:84
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:209
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:391
iso_collapsed_lifetimes_update
void iso_collapsed_lifetimes_update(long ipISO, long nelem)
Definition: iso_create.cpp:1697
max_num_levels
long int max_num_levels
Definition: iso.cpp:13
FillExtraLymanLine
STATIC void FillExtraLymanLine(const TransitionList::iterator &t, long ipISO, long nelem, long nHi)
Definition: iso_create.cpp:1104
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_iso_sp::n_HighestResolved_local
long int n_HighestResolved_local
Definition: iso.h:538
nLine
static long int nLine
Definition: save_line.cpp:271
multi_arr::clone
const multi_geom< d, ALLOC > & clone() const
Definition: container_classes.h:1779
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:554
lines_service.h
t_iso_sp::BranchRatio
multi_arr< double, 2 > BranchRatio
Definition: iso.h:480
helike_energy
double helike_energy(long nelem, long ipLev)
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:570
iso_update_num_levels
void iso_update_num_levels(long ipISO, long nelem)
Definition: iso_create.cpp:1618
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:536
L_
#define L_(A_)
Definition: iso.h:23
chL
char chL[21]
Definition: iso_create.cpp:37
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:477
POW2
#define POW2
Definition: cddefines.h:979
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:352
energy
double energy(const genericState &gs)
Definition: generic_state.cpp:49
ipH1s
const int ipH1s
Definition: iso.h:29
multi_arr
Definition: container_classes.h:938
STATIC
#define STATIC
Definition: cddefines.h:118
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:610
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:447
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
iso_satellite_update
void iso_satellite_update(long nelem)
Definition: iso_create.cpp:1511
N_
#define N_(A_)
Definition: iso.h:22
t_isoCTRL::chISO
const char * chISO[NISO]
Definition: iso.h:348
TransitionProxy::ipEmis
int & ipEmis() const
Definition: transition.h:455
M
static const int M
Definition: thirdparty.cpp:3955
mole
t_mole_local mole
Definition: mole.cpp:8
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:789
multi_arr::invalidate
void invalidate()
Definition: container_classes.h:1054
t_iso_sp::tr
TransitionList * tr
Definition: iso.h:483
iso_zero
STATIC void iso_zero(void)
Definition: iso_create.cpp:481
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:174
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:385
realnum
float realnum
Definition: cddefines.h:124
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:468
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
ipLY_A
const int ipLY_A
Definition: cddefines.h:346
t_prt_matrix::species
char species[CHARS_SPECIES]
Definition: prt.h:116
max
long max(int a, long b)
Definition: cddefines.h:817
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
molecule::index
int index
Definition: mole.h:194
S_
#define S_(A_)
Definition: iso.h:24
powi
double powi(double, long int)
Definition: service.cpp:594
helike_transprob
realnum helike_transprob(long nelem, long ipHi, long ipLo)
Definition: helike_einsta.cpp:1283
iso_collapsed_Aul_update
void iso_collapsed_Aul_update(long ipISO, long nelem)
Definition: iso_create.cpp:1649
SatelliteLines
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:35
ipDOUBLET
Definition: iso.h:84
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
prt
t_prt prt
Definition: prt.cpp:14
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:160
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
t_iso_sp
Definition: iso.h:470
iso_allocate
STATIC void iso_allocate(void)
Definition: iso_create.cpp:523
t_isoCTRL::nLyaLevel
int nLyaLevel[NISO]
Definition: iso.h:397
hydro_bauman.h
ipH2p
const int ipH2p
Definition: iso.h:31
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
makeChemical
string makeChemical(long nelem, long ion)
Definition: species.cpp:929
t_iso_sp::CascadeProb
multi_arr< double, 2 > CascadeProb
Definition: iso.h:479
t_iso_sp::ex
multi_arr< extra_tr, 2 > ex
Definition: iso.h:478
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
DoFSMixing
void DoFSMixing(long nelem, long ipLoSing, long ipHiSing)
Definition: helike_einsta.cpp:1097
t_isoCTRL::ipLyaRedist
int ipLyaRedist[NISO]
Definition: iso.h:394
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1077
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:92
iso_recomb_setup
void iso_recomb_setup(long ipISO)
Definition: iso_radiative_recomb.cpp:842
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
iso_recomb_malloc
void iso_recomb_malloc(void)
Definition: iso_radiative_recomb.cpp:769
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
LIMELM
const int LIMELM
Definition: cddefines.h:308
pow2
T pow2(T a)
Definition: cddefines.h:981
t_prt::matrix
t_prt_matrix matrix
Definition: prt.h:238
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:50
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:630
iso_state_lifetime
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
Definition: iso_create.cpp:1189
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:350
iso_cascade
void iso_cascade(long ipISO, long nelem)
Definition: iso_create.cpp:1233
t_phycon::te_ryd
double te_ryd
Definition: phycon.h:27
iso_create
void iso_create(void)
Definition: iso_create.cpp:111
iso_setRedisFun
void iso_setRedisFun(long ipISO, long nelem, long ipLo, long ipHi)
Definition: iso_create.cpp:48
t_dense::eden
double eden
Definition: dense.h:201
ipSINGLET
Definition: iso.h:84
t_isoCTRL::ipSubRedist
int ipSubRedist[NISO]
Definition: iso.h:394
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
t_iso_sp::lgPrtMatrix
bool lgPrtMatrix
Definition: iso.h:603
t_isoCTRL::stat_ion
realnum stat_ion[NISO]
Definition: iso.h:382
t_phycon::alogte
double alogte
Definition: phycon.h:92
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1826
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:9
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:524
fixit
#define fixit(a)
Definition: cddefines.h:417
HelikeTransProbSetup
void HelikeTransProbSetup(void)
Definition: helike_einsta.cpp:1333
t_phycon::te
double te
Definition: phycon.h:21
iso_setOpacity
void iso_setOpacity(long ipISO, long nelem, long ipLo, long ipHi)
Definition: iso_create.cpp:82
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:349
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:690
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:2356
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:529
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:519
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:468
t_isoCTRL::nLyman_malloc
long int nLyman_malloc[NISO]
Definition: iso.h:352
AnonStates
qList AnonStates("AnonStates", 1)
t_phycon::te32
double te32
Definition: phycon.h:58
ipCRDW
const int ipCRDW
Definition: cddefines.h:344
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:631
isoSpecies
static vector< species > isoSpecies(NISO)
iso_recomb_auxiliary_free
void iso_recomb_auxiliary_free(void)
Definition: iso_radiative_recomb.cpp:829
hydro_transprob
realnum hydro_transprob(long nelem, long ipHi, long ipLo)
Definition: hydroeinsta.cpp:46
Generated by   doxygen 1.8.5