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

Generated for cloudy by doxygen 1.8.1.1