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