cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_createpointers.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 /*ContCreatePointers set up pointers for lines and continua called by cloudy after input read in
4  * and continuum mesh has been set */
5 /*ipShells assign continuum energy pointers to shells for all atoms,
6  * called by ContCreatePointers */
7 /*LimitSh sets upper energy limit to subshell integrations */
8 /*ContBandsCreate - read set of continuum bands to enter total emission into line stack*/
9 #include "cddefines.h"
10 #include "iso.h"
11 #include "secondaries.h"
12 #include "taulines.h"
13 #include "elementnames.h"
14 #include "ionbal.h"
15 #include "rt.h"
16 #include "opacity.h"
17 #include "yield.h"
18 #include "he.h"
19 #include "fe.h"
20 #include "rfield.h"
21 #include "oxy.h"
22 #include "trace.h"
23 #include "hmi.h"
24 #include "heavy.h"
25 #include "atmdat_adfa.h"
26 #include "ipoint.h"
27 #include "h2.h"
28 #include "continuum.h"
29 #include "freebound.h"
30 #include "two_photon.h"
31 #include "dense.h"
32 #include "lines_service.h"
33 
34 /*LimitSh sets upper energy limit to subshell integrations */
35 STATIC long LimitSh(long int ion,
36  long int nshell,
37  long int nelem);
38 
39 STATIC void ipShells(
40  /* nelem is the atomic number on the C scale, Li is 2 */
41  long int nelem);
42 
43 /*ContBandsCreate - read set of continuum bands to enter total emission into line*/
45  /* chFile is optional filename, if void then use default bands,
46  * if not void then use file specified,
47  * return value is 0 for success, 1 for failure */
48  const char chFile[] );
49 
50 STATIC inline void print_emline_fine( const char *LineGroup, const TransitionProxy &tr )
51 {
52  fprintf( ioQQQ, "%-10s -> '%s'\t Energy Ang= %.4e\t Ryd= %.4e\t Fine= %ld\t fine_ener= %.4e\n",
53  LineGroup,
54  tr.chLabel().c_str(),
55  tr.EnergyAng(),
56  tr.EnergyRyd(),
57  tr.Emis().ipFine(),
58  rfield.fine_anu[ tr.Emis().ipFine() ] );
59 }
60 
62 {
63  /* counter to say whether pointers have ever been evaluated */
64  static int nCalled = 0;
65 
66  DEBUG_ENTRY( "ContCreatePointers()" );
67 
68  /* create the hydrogen atom for this core load, routine creates space then zeros it out
69  * on first call, on second and later calls it only zeros things out */
70  iso_create();
71 
72  /* create internal static variables needed to do the H2 molecule */
73  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
74  (*diatom)->init();
75 
76  /* nCalled is local static variable defined 0 when defined.
77  * it is incremented below, so that space only allocated one time per coreload. */
78  if( nCalled > 0 )
79  {
80  if( trace.lgTrace )
81  fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
82  return;
83  }
84  else
85  {
86  if( trace.lgTrace )
87  fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
88  ++nCalled;
89  }
90 
91  for( long i=0; i < rfield.nflux_with_check; i++ )
92  {
93  /* this is array of labels for lines and continua, set to blanks at first */
94  rfield.chContLabel[i] = " ";
95  rfield.chLineLabel[i] = " ";
96  }
97 
98  /* we will generate a set of array indices to ionization edges for
99  * the first thirty elements. First set all array indices to
100  * totally bogus values so we will crash if misused */
101  for( long nelem=0; nelem<LIMELM; ++nelem )
102  {
103  if( dense.lgElmtOn[nelem] )
104  {
105  for( long ion=0; ion<LIMELM; ++ion )
106  {
107  for( long nshells=0; nshells<7; ++nshells )
108  {
109  for( long j=0; j<3; ++j )
110  {
111  opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
112  }
113  }
114  }
115  }
116  }
117 
118  /* pointer to excited state of O+2 */
119  opac.o3exc = 3.855;
120  opac.ipo3exc[0] = ipContEnergy(opac.o3exc,"O3ex");
121 
122  /* main hydrogenic arrays - THIS OCCURS TWICE!! H and He here, then the
123  * remaining hydrogenic species near the bottom. This is so that H and HeII get
124  * their labels stuffed into the arrays, and the rest of the hydrogenic series
125  * get whatever is left over after the level 1 lines.
126  * to find second block, search on "ipZ=2" */
127  /* NB note that no test for H or He being on exists here - we will always
128  * define the H and He arrays even when He is off, since we need to
129  * know where the he edges are for the bookkeeping that occurs in continuum
130  * binning routines */
131  /* this loop is over H, He-like only - DO NOT CHANGE TO NISO */
132  for( long ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
133  {
134  /* this will be over HI, HeII, then HeI only */
135  for( long nelem=ipISO; nelem < 2; nelem++ )
136  {
137  /* generate label for this ion */
138  string chLab = chIonLbl( nelem+1, nelem+1-ipISO );
139 
140  /* array index for continuum edges for ground */
141  iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd, chLab.c_str());
142  for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
143  {
144  /* array index for continuum edges for excited levels */
145  iso_sp[ipISO][nelem].fb[ipHi].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd, chLab.c_str());
146 
147  /* define all line array indices */
148  for( long ipLo=0; ipLo < ipHi; ipLo++ )
149  {
150  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
151  continue;
152 
153  /* some lines have negative or zero energy */
154  /* >>chng 03 apr 22, this check was if less than or equal to zero,
155  * changed to lowest energy point so that ultra low energy transitions are
156  * not considered. */
157  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() < rfield.emm() )
158  continue;
159 
160  /* some energies are negative for inverted levels */
161  iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() =
162  ipLineEnergy(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd(), chLab.c_str(),
163  iso_sp[ipISO][nelem].fb[ipLo].ipIsoLevNIonCon);
164  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine() =
165  ipFineCont(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() );
166  /* check that energy scales are the same, to within energy resolution of arrays */
167 # ifndef NDEBUG
168  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine() > 0 )
169  {
170  realnum anuCoarse = rfield.anu(iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1);
171  realnum anuFine = rfield.fine_anu[iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine()];
172  realnum widCoarse = rfield.widflx(iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1);
173  realnum widFine = anuFine - rfield.fine_anu[iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine()-1];
174  realnum width = MAX2( widFine , widCoarse );
175  /* NB - can only assert more than width of coarse cell
176  * since close to ionization limit, coarse index is
177  * kept below next ionization edge
178  * >>chng 05 mar 02, pre factor below had been 1.5, chng to 2
179  * tripped near H grnd photo threshold */
180  ASSERT( fabs(anuCoarse - anuFine) / anuCoarse <
181  2.*width/anuCoarse);
182  }
183 # endif
184  }
185  }/* ipHi loop */
186  }/* nelem loop */
187  }/* ipISO */
188 
189  /* make sure the HeII Balmer pointers are above the HI Lyman pointer
190  *
191  * H and He are special because of their large abundances -
192  * the OTS fields they produce are very important in establishing the ionization in the gas.
193  * The OTS fields are placed in a single cell at an element's threshold.
194  * When that shell's photoionization rate is evaluated the first cell,
195  * with the OTS field, is not included - that would be double counting.
196  * The self-OTS is taken into account with changes to the recombination coefficient.
197  *
198  * Note that this code makes implicit assumptions about the abundance pattern, which is a bug
199  *
200  * The edges in the frequency mesh (in mesh.cpp) are chosen such that the test below should
201  * always pass. If these ASSERTs trip, update the energies of the H end He edges in mesh.cpp
202  * to be the same as the ones used here (i.e. iso_sp[][].fb[].xIsoLevNIonRyd). */
203  ASSERT( iso_sp[ipH_LIKE][ipHELIUM].fb[ipH2s].ipIsoLevNIonCon > iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon );
204  ASSERT( iso_sp[ipH_LIKE][ipHELIUM].fb[ipH2p].ipIsoLevNIonCon > iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon );
205 
206  /* we will define an array of either 1 or 0 to show whether photooelectron
207  * energy is large enough to produce secondary ionizations
208  * below 100eV no secondary ionization - this is the threshold*/
209  secondaries.ipSecIon = ipoint(7.353);
210 
211  /* this is highest energy where k-shell opacities are counted
212  * can be adjusted with "set kshell" command */
214 
215  /* pointers for molecules
216  * H2+ dissociation energy 2.647 eV but cs small below 0.638 Ryd */
217  opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
218  opac.ih2pnt[1] = ipoint(1.03);
219  // Excited state H2+
220  opac.ih2pnt_ex[0] = ipContEnergy(0.052,"H2+*");
221  opac.ih2pnt_ex[1] = ipoint(1.03);
222 
223  //pointers for most prominent PAH features
224  {
225  /* energies given to ipContEnergy are only to put lave in the right place
226  * wavelengths are rough observed values of blends
227  * 7.6 microns */
228  (void) ipContEnergy(0.0117, "PAH " );
229 
230  /* feature near 6.2 microns */
231  (void) ipContEnergy(0.0147, "PAH " );
232 
233  /* 3.3 microns */
234  (void) ipContEnergy(0.028, "PAH " );
235 
236  /* 11.2 microns */
237  (void) ipContEnergy(0.0080, "PAH " );
238 
239  /* 12.3 microns */
240  (void) ipContEnergy(0.0077, "PAH " );
241 
242  /* 13.5 microns */
243  (void) ipContEnergy(0.0069, "PAH " );
244  }
245 
246  /* fix pointers for hydrogen and helium */
247 
248  /* pointer to Bowen 374A resonance line */
249  he.ip660 = ipLineEnergy(1.38,"He 2",0);
250 
251  /* pointer to energy defining effect x-ray column */
252  rt.ipxry = ipoint(73.5);
253 
254  /* pointer to Hminus edge at 0.754eV */
255  hmi.iphmin = ipContEnergy(HMINUSIONPOT,"H- ");
256 
257  /* pointer to threshold for H2 photoionization at 15.4 eV */
258  fixit("need to generalize ionization energy and label!");
259  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
260  (*diatom)->ip_photo_opac_thresh = ipContEnergy( 15.4/EVRYD , "H2 ");
261 
262  hmi.iheh1 = ipoint(1.6);
263  hmi.iheh2 = ipoint(2.3);
264 
265  /* pointer to carbon k-shell ionization */
266  opac.ipCKshell = ipoint(20.6);
267 
268  /* pointer to threshold for pair production */
269  opac.ippr = ipoint(7.51155e4) + 1;
270 
271  /* pointer to x-ray - gamma ray bound; 100 kev */
273 
274  /* confirm that labels are in correct location */
275  ASSERT( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1].c_str(), "H 1" ) ==0 );
276  ASSERT( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon-1].c_str(), "H 1" ) ==0 );
277 
278  ASSERT( strcmp( rfield.chContLabel[iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1].c_str(), "He 2" ) ==0 );
279 
280  /* these are indices for centers of B and V filters,
281  * taken from table on page 202 of Allen, AQ, 3rd ed */
282  /* the B filter array offset */
283  rfield.ipB_filter = ipoint( RYDLAM / WL_B_FILT );
284  /* the V filter array offset */
285  rfield.ipV_filter = ipoint( RYDLAM / WL_V_FILT );
286 
287  /* these are the lower and upper bounds for the G0 radiation field
288  * used by Tielens & Hollenbach in their PDR work */
289  rfield.ipG0_TH85_lo = ipoint( 6.0 / EVRYD );
290  rfield.ipG0_TH85_hi = ipoint( 13.6 / EVRYD );
291 
292  /* this is the limits for Draine & Bertoldi Habing field */
293  rfield.ipG0_DB96_lo = ipoint( RYDLAM / 1110. );
294  rfield.ipG0_DB96_hi = ipoint( RYDLAM / 912. );
295 
296  /* this is special form of G0 that could be used in some future version, for now,
297  * use default, TH85 */
298  rfield.ipG0_spec_lo = ipoint( 6.0 / EVRYD );
299  rfield.ipG0_spec_hi = ipoint( RYDLAM / 912. );
300 
301  /* this index is to 1000A to obtain the extinction at 1000A */
302  rfield.ip1000A = ipoint( RYDLAM / 1000. );
303 
304  /* following order of elements is in roughly decreasing abundance
305  * when ipShells gets a cell for the valence IP it does so through
306  * ipContEnergy, which makes sure that no two ions get the same cell
307  * earliest elements have most precise ip mapping */
308 
309  /* set up shell pointers for hydrogen */
310  {
311  long nelem = ipHYDROGEN;
312  long ion = 0;
313 
314  /* the number of shells */
315  Heavy.nsShells[nelem][0] = 1;
316 
317  /*pointer to ionization threshold in energy array*/
318  Heavy.ipHeavy[nelem][ion] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
319  opac.ipElement[nelem][ion][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
320 
321  /* upper limit to energy integration */
322  opac.ipElement[nelem][ion][0][1] = rfield.nflux_with_check;
323  }
324 
325  if( dense.lgElmtOn[ipHELIUM] )
326  {
327  /* set up shell pointers for helium */
328  long nelem = ipHELIUM;
329  long ion = 0;
330 
331  /* the number of shells */
332  Heavy.nsShells[nelem][0] = 1;
333 
334  /*pointer to ionization threshold in energy array*/
335  Heavy.ipHeavy[nelem][ion] = iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon;
336  opac.ipElement[nelem][ion][0][0] = iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon;
337 
338  /* upper limit to energy integration */
339  opac.ipElement[nelem][ion][0][1] = rfield.nflux_with_check;
340 
341  /* (hydrogenic) helium ion */
342  ion = 1;
343  /* the number of shells */
344  Heavy.nsShells[nelem][1] = 1;
345 
346  /*pointer to ionization threshold in energy array*/
347  Heavy.ipHeavy[nelem][ion] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
348  opac.ipElement[nelem][ion][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;
349 
350  /* upper limit to energy integration */
351  opac.ipElement[nelem][ion][0][1] = rfield.nflux_with_check;
352  }
353 
354  /* check that ionization potential of neutral carbon valence shell is
355  * positive */
356  ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
357 
358  /* now fill in all sub-shell ionization array indices for elements heavier than He,
359  * this must be done after previous loop on iso.ipIsoLevNIonCon[ipH_LIKE] since hydrogenic species use
360  * iso.ipIsoLevNIonCon[ipH_LIKE] rather than ipoint in getting array index within continuum array */
361  for( long i=NISO; i<LIMELM; ++i )
362  {
363  /* i is the atomic number on the c scale, 2 for Li */
364  ipShells(i);
365  }
366 
367  /* most of these are set in ipShells, but not h-like or he-like, so do these here */
368  Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD* 0.9998787;
369  Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD* 0.9998787;
370  Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD* 0.9998787;
371  for( long nelem=2; nelem<LIMELM; ++nelem )
372  {
373  Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD* 0.9998787;
374  Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
375  if( dense.lgElmtOn[nelem])
376  {
377  /* now confirm that all are properly set */
378  for( long j=0; j<=nelem; ++j )
379  {
380  ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
381  }
382  for( long j=0; j<nelem; ++j )
383  {
384  ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
385  }
386  }
387  }
388 
389  /* array indices to bound Compton electron recoil ionization of all elements */
390  for( long nelem=0; nelem<LIMELM; ++nelem )
391  {
392  if( dense.lgElmtOn[nelem])
393  {
394  for( long ion=0; ion<nelem+1; ++ion )
395  {
396  /* this is the threshold energy to Compton ionize valence shell electrons */
397  double energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
398  /* the array index for this energy */
399  ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
400  }
401  }
402  }
403 
404  /* oxygen pointers for excited states
405  * IO3 is pointer to O++ exc state, is set above */
406  oxy.i2d = ipoint(1.242);
407  oxy.i2p = ipoint(1.367);
408  opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");//2s^2 2p^4, ^1D level, J=2, energy relative to ground level is ~0.1446 Ry
409  opac.ipo1exc[1] = ipoint(2.0);// upper limit to range of energies where opacity for 1D absorption is defined
410 
411  /* upper limit for excited state photoionization
412  * do not use ipContEnergy since only upper limit */
413  opac.ipo3exc[1] = ipoint(5.0);
414 
415  /* upper level of 4363 */
416  opac.o3exc3 = 3.646;
417  opac.ipo3exc3[0] = ipContEnergy(opac.o3exc3,"O3ex");
418  opac.ipo3exc3[1] = ipoint(5.0);
419 
420  /* >>chng 97 jan 27, move nitrogen after oxygen so that O gets the
421  * most accurate pointers
422  * Nitrogen
423  * in1(1) is thresh for photo from excited state */
424  opac.in1[0] = ipContEnergy(0.893,"N1ex");
425 
426  /* upper limit */
427  opac.in1[1] = ipoint(2.);
428 
430  {
431  fprintf( ioQQQ, " ContCreatePointers:%ld energy cells used. N(1R):%4ld N(1.8):%4ld N(4Ryd):%4ld N(O3)%4ld N(x-ray):%5ld N(rcoil)%5ld\n",
433  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon,
434  iso_sp[ipHE_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon,
435  iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon,
436  opac.ipo3exc[0],
437  opac.ipCKshell,
439 
440 
441  fprintf( ioQQQ, " ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n",
443 
444  fprintf( ioQQQ, " ContCreatePointers: H pointers;" );
445  for( long i=0; i <= 6; i++ )
446  {
447  fprintf( ioQQQ, "%4ld%4ld", i, iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].ipIsoLevNIonCon );
448  }
449  fprintf( ioQQQ, "\n" );
450 
451  fprintf( ioQQQ, " ContCreatePointers: Oxy pnters;" );
452 
453  for( long i=1; i <= 8; i++ )
454  {
455  fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
456  }
457  fprintf( ioQQQ, "\n" );
458 
459  }
460 
461  /* Magnesium
462  * following is energy for phot of MG+ from exc state producing 2798 */
463  opac.ipmgex = ipoint(0.779);
464 
465  /* lower, upper edges of Ca+ excited term photoionizaiton */
466  opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
467  opac.ica2ex[1] = ipoint(1.);
468 
469  /* set up factors and pointers for Fe continuum fluorescence */
470  fe.ipfe10 = ipoint(2.605);
471 
472  /* following is WL(CM)**2/(8PI) * conv fac for RYD to NU *A21 */
473  fe.pfe10 = (realnum)(2.00e-18/rfield.widflx(fe.ipfe10-1));
474 
475  /* this is 353 pump, f=0.032 */
476  fe.pfe11a = (realnum)(4.54e-19/rfield.widflx(fe.ipfe10-1));
477 
478  /* this is 341.1 f=0.012 */
479  fe.pfe11b = (realnum)(2.75e-19/rfield.widflx(fe.ipfe10-1));
480  fe.pfe14 = (realnum)(1.15e-18/rfield.widflx(fe.ipfe10-1));
481 
482  /* set up energy pointers for line optical depth arrays
483  * this also increments flux, sets other parameters for lines */
484 
485  /*Beginning of the dBaseLines*/
486  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
487  {
488  /*Put null terminated line label into chLab*/
489  char chLab[NCHLAB];
490  strncpy(chLab,dBaseSpecies[ipSpecies].chLabel,NCHLAB-1);
491  chLab[NCHLAB-1]='\0';
492 
493  for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
494  em != dBaseTrans[ipSpecies].Emis().end(); ++em)
495  {
496  /* upper level lifetime to calculate the damping parameter.*/
497  (*em).dampXvel() = (realnum)(1./
498  dBaseStates[ipSpecies][em->Tran().ipHi()].lifetime()/em->Tran().EnergyWN()/PI4);
499  (*em).damp() = -1000.0;
500 
501  // Data intended for ADAS have a large number of transitions with Aul = 1e-30 - these are
502  // not real, but are there to keep the ADAS codes happy. We do not want to transfer them.
503  static const double minAul = 1e-29;
504  if( (*em).Aul() > minAul )
505  {
506  (*em).Tran().ipCont() = ipLineEnergy((*em).Tran().EnergyRyd(), chLab ,0);
507  (*em).ipFine() = ipFineCont((*em).Tran().EnergyRyd() );
508  }
509  else
510  {
511  (*em).Tran().ipCont() = -1;
512  (*em).ipFine() = -1;
513  }
514 
515  // these are edge cases which did pass the minAul above, but still underflowed
516  // when converted to (real)gf - conversion depends on energy of transition so
517  // Aul - gf mapping is not linear
518  if ((*em).gf() <= 0.0)
519  {
520  (*em).Tran().ipCont() = -1;
521  (*em).ipFine() = -1;
522  }
523 
524  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
525  (*em).opacity() =
526  (realnum)(abscf((*em).gf(),(*em).Tran().EnergyWN(),
527  (*(*em).Tran().Lo()).g()));
528  }
529  }
530  /*end of the dBaseLines*/
531 
532  /* set the ipCont struc element for the H2 molecule */
533  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
534  (*diatom)->H2_ContPoint();
535 
536  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
537  {
538  /* do remaining part of the iso sequences */
539  for( long nelem=2; nelem < LIMELM; nelem++ )
540  {
541  if( dense.lgElmtOn[nelem])
542  {
543  string chLab = chIonLbl( nelem+1, nelem+1-ipISO );
544 
545  /* array index for continuum edges */
546  iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon =
547  ipContEnergy(iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd, chLab.c_str());
548 
549  for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
550  {
551  /* array index for continuum edges */
552  iso_sp[ipISO][nelem].fb[ipHi].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd, chLab.c_str());
553 
554  /* define all line pointers */
555  for( long ipLo=0; ipLo < ipHi; ipLo++ )
556  {
557 
558  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
559  continue;
560 
561  /* some lines have negative or zero energy */
562  /* >>chng 03 apr 22, this check was if less than or equal to zero,
563  * changed to lowest energy point so that ultra low energy transitions are
564  * not considered. */
565  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() < rfield.emm() )
566  continue;
567 
568  iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() =
569  ipLineEnergy(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd(), chLab.c_str(),
570  iso_sp[ipISO][nelem].fb[ipLo].ipIsoLevNIonCon);
571  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ipFine() =
572  ipFineCont(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyRyd() );
573  }
574  }
575  iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd, chLab.c_str());
576  }
577  }
578  }
579  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
580  {
581  /* this will be over HI, HeII, then HeI only */
582  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
583  {
584  if( dense.lgElmtOn[nelem])
585  {
586  /* these are the extra Lyman lines */
587  for( long ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
588  {
589  long ipLo = 0;
590  /* some energies are negative for inverted levels */
591  char chLab[NCHLAB];
592  strncpy(chLab,"LyEx",NCHLAB-1);
593  chLab[NCHLAB-1]='\0';
594  TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
595  (*tr).ipCont() =
596  ipLineEnergy((*tr).EnergyRyd() , chLab ,
597  iso_sp[ipISO][nelem].fb[ipLo].ipIsoLevNIonCon);
598 
599  (*tr).Emis().ipFine() =
600  ipFineCont((*tr).EnergyRyd() );
601  }
602 
603  if( iso_ctrl.lgDielRecom[ipISO] )
604  {
605  ASSERT( ipISO>ipH_LIKE );
606  for( long ipLo=0; ipLo<iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
607  {
608  char chLab[NCHLAB];
609  strncpy(chLab,"SatL",NCHLAB-1);
610  chLab[NCHLAB-1]='\0';
611  SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].ipCont() = ipLineEnergy(
612  SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].EnergyRyd() , chLab ,
613  0);
614 
615  SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].Emis().ipFine() =
616  ipFineCont(SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]].EnergyRyd() );
617  }
618  }
619  }
620  }
621  }
622 
623  fixit("is this redundant?");
624  /* for He-like sequence the majority of the transitions are bogus - A set to special value in this case */
625  {
626  long ipISO = ipHE_LIKE;
627  /* do remaining part of the he iso sequence */
628  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
629  {
630  if( dense.lgElmtOn[nelem])
631  {
632  for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
633  {
634  for( long ipLo=0; ipLo < ipHi; ipLo++ )
635  {
636  TransitionProxy tr = iso_sp[ipISO][nelem].trans(ipHi,ipLo);
637  if( tr.ipCont() <= 0 )
638  continue;
639 
640  if( fabs(tr.Emis().Aul() - iso_ctrl.SmallA) < SMALLFLOAT )
641  {
642  /* iso_ctrl.SmallA is value assigned to bogus transitions */
643  tr.ipCont() = -1;
644  tr.Emis().ipFine() = -1;
645  }
646  }
647  }
648  }
649  }
650  }
651 
652  /* inner shell transitions */
653  for( size_t i=0; i<UTALines.size(); ++i )
654  {
655  if( UTALines[i].Emis().Aul() > 0. )
656  {
657 
658  // dampXvel is derived in atmdat_readin because autoionization rates
659  // must be included in the total rate for the damping constant
660  ASSERT( UTALines[i].Emis().dampXvel() >0. );
661 
662  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
663  UTALines[i].Emis().opacity() =
664  (realnum)(abscf( UTALines[i].Emis().gf(), UTALines[i].EnergyWN(), (*UTALines[i].Lo()).g()));
665 
666  /* get pointer to energy in continuum mesh */
667  UTALines[i].ipCont() = ipLineEnergy(UTALines[i].EnergyRyd(), chIonLbl(UTALines[i]).c_str(),0 );
668  UTALines[i].Emis().ipFine() = ipFineCont(UTALines[i].EnergyRyd() );
669  {
670  enum{ DEBUG_LOC = false };
671  if( DEBUG_LOC && UTALines[i].chLabel() == "Ar 7 43.5239A" )
672  {
673  print_emline_fine( "UTA", UTALines[i] );
674  }
675  }
676 
677  /* find heating per absorption,
678  * first find threshold for this shell in ergs */
679  /* ionization threshold in erg */
680  double thresh = Heavy.Valence_IP_Ryd[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1] *EN1RYD;
681  UTALines[i].Coll().heat() = (UTALines[i].EnergyErg()-thresh);
682  ASSERT( UTALines[i].Coll().heat()> 0. );
683  }
684  }
685 
686  /* level 2 heavy element lines */
687  for( long i=0; i < nWindLine; i++ )
688  {
689  /* derive the A, call to function is gf, wl (A), g_up */
690  TauLine2[i].Emis().Aul() =
691  (realnum)(eina(TauLine2[i].Emis().gf(),
692  TauLine2[i].EnergyWN(),(*TauLine2[i].Hi()).g()));
693 
694  /* coefficient needed for damping constant - units cm s-1 */
695  TauLine2[i].Emis().dampXvel() =
696  (realnum)(TauLine2[i].Emis().Aul()/
697  TauLine2[i].EnergyWN()/PI4);
698 
699  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
700  TauLine2[i].Emis().opacity() =
701  (realnum)(abscf(TauLine2[i].Emis().gf(),
702  TauLine2[i].EnergyWN(),(*TauLine2[i].Lo()).g()));
703 
704  /* get pointer to energy in continuum mesh */
705  TauLine2[i].ipCont() = ipLineEnergy(TauLine2[i].EnergyRyd(), chIonLbl(TauLine2[i]).c_str(),0 );
706  TauLine2[i].Emis().ipFine() = ipFineCont(TauLine2[i].EnergyRyd() );
707  /*if( TauLine2[i].ipCont()==751 )
708  fprintf(ioQQQ,"( atom_level2 %s\n", chLab);*/
709  }
710 
711  /* hyperfine structure lines */
712  for( size_t i=0; i < HFLines.size(); i++ )
713  {
714  ASSERT( HFLines[i].Emis().Aul() > 0. );
715  /* coefficient needed for damping constant */
716  HFLines[i].Emis().dampXvel() =
717  (realnum)(HFLines[i].Emis().Aul()/
718  HFLines[i].EnergyWN()/PI4);
719  HFLines[i].Emis().damp() = 1e-20f;
720 
721  /* derive the abs coefficient, call to function is gf, wl (A), g_low */
722  HFLines[i].Emis().opacity() =
723  (realnum)(abscf(HFLines[i].Emis().gf(),
724  HFLines[i].EnergyWN(),
725  (*HFLines[i].Lo()).g()));
726  /* gf from this and 21 cm do not agree, A for HFS is 10x larger than level1 dat */
727  /*fprintf(ioQQQ,"HFLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
728  i,HFLines[i].Emis().opacity() , HFLines[i].Emis().gf() , HFLines[i].Emis().Aul() , HFLines[i].EnergyWN(),(*HFLines[i].Lo()).g());*/
729 
730  /* get pointer to energy in continuum mesh */
731  HFLines[i].ipCont() = ipLineEnergy(HFLines[i].EnergyRyd() , chIonLbl(HFLines[i]).c_str(),0 );
732  HFLines[i].Emis().ipFine() = ipFineCont(HFLines[i].EnergyRyd() );
733  }
734 
735  /* the group of inner shell fluorescent lines */
736  for( long i=0; i < t_yield::Inst().nlines(); ++i )
737  {
738  string chLab = chIonLbl( t_yield::Inst().nelem(i)+1, t_yield::Inst().ion_emit(i)+1 );
739  t_yield::Inst().set_ipoint( i, ipLineEnergy( t_yield::Inst().energy(i) , chLab.c_str() , 0 ) );
740  }
741 
742  /* ================================================================================== */
743  /* two photon two-photon 2-nu 2 nu 2 photon 2-photon */
744 
745  /* now loop over the two iso-sequences with two photon continua */
746  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
747  {
748  /* set up two photon emission */
749  for( long nelem=ipISO; nelem<LIMELM; ++nelem )
750  {
751  if( dense.lgElmtOn[nelem] )
752  {
753  // upper level for two-photon emission in H and He iso sequences
754  // The 1+ipISO is not rigorous, it just works for H- and He-like
755  const int TwoS = (1+ipISO);
756  double Aul;
757  /* 2s two photon */
758  if( ipISO==ipH_LIKE )
759  Aul = 8.226*powi((double)(nelem+1.),6);
760  else
761  {
762  ASSERT( ipISO==ipHE_LIKE );
763  /* >>refer Helike 2pho Derevianko, A., & Johnson, W. R. 1997, Phys. Rev. A, 56, 1288
764  * numbers are not explicitly given in this paper for Z=21-24,26,27,and 29.
765  * So numbers given here are interpolated. */
766  fixit("where is 51.02 from? Value is 51.3 from the Derevianko & Johnson paper cited above.");
767  const double As2nuFrom1S[29] = {51.02,1940.,1.82E+04,9.21E+04,3.30E+05,9.44E+05,2.31E+06,5.03E+06,1.00E+07,
768  1.86E+07,3.25E+07,5.42E+07,8.69E+07,1.34E+08,2.02E+08,2.96E+08,4.23E+08,5.93E+08,8.16E+08,
769  1.08E+09,1.43E+09,1.88E+09,2.43E+09,3.25E+09,3.95E+09,4.96E+09,6.52E+09,7.62E+09,9.94E+09};
770  Aul = As2nuFrom1S[nelem-1];
771  }
772 
773  TwoPhotonSetup( iso_sp[ipISO][nelem].TwoNu, TwoS, 0,
774  Aul,
775  iso_sp[ipISO][nelem].trans(TwoS,0),
776  ipISO, nelem );
777  }
778  }
779  }
780 
781  // add He-like 2nu 2^3S - 1^1S
782  {
783  const long ipISO = ipHE_LIKE;
784  for( long nelem=ipISO; nelem<LIMELM; ++nelem )
785  {
786  if( dense.lgElmtOn[nelem] )
787  {
788  /* Important clarification, according to Derevianko & Johnson (see ref above), 2^3S can decay
789  * to ground in one of two ways: through a two-photon process, or through a single-photon M1 decay,
790  * but the M1 rates are about 10^4 greater that the two-photon decays throughout the entire
791  * sequence. Thus these numbers, are much weaker than the effective decay rate, but should probably
792  * be treated in as a two-photon decay at some point */
793  // >> refer He As Drake, G. W. F., Victor, G. A., & Dalgarno, A. 1969, Physical Review, 180, 25
794  const double As2nuFrom3S[29] = {4.09e-9,1.25E-06,5.53E-05,8.93E-04,8.05E-03,4.95E-02,2.33E-01,8.94E-01,2.95E+00,
795  8.59E+00,2.26E+01,5.49E+01,1.24E+02,2.64E+02,5.33E+02,1.03E+03,1.91E+03,3.41E+03,5.91E+03,
796  9.20E+03,1.50E+04,2.39E+04,3.72E+04,6.27E+04,8.57E+04,1.27E+05,2.04E+05,2.66E+05,4.17E+05};
797 
798  TwoPhotonSetup( iso_sp[ipISO][nelem].TwoNu, ipHe2s3S, ipHe1s1S,
799  As2nuFrom3S[nelem-1],
800  iso_sp[ipISO][nelem].trans(ipHe2s3S,ipHe1s1S),
801  ipISO, nelem );
802  }
803  }
804  }
805 
806  {
807  /* this is an option to print out one of the two photon continua */
808  enum {DEBUG_LOC=false};
809  if( DEBUG_LOC )
810  {
811  const int nCRS = 21;
812  double ener[nCRS]={
813  0., 0.03738, 0.07506, 0.1124, 0.1498, 0.1875,
814  0.225, 0.263, 0.300, 0.3373, 0.375, 0.4127,
815  0.4500, 0.487, 0.525, 0.5625, 0.6002, 0.6376,
816  0.6749, 0.7126, 0.75};
817 
818  long nelem = ipHYDROGEN;
819  long ipISO = ipHYDROGEN;
820  two_photon& tnu = iso_sp[ipISO][nelem].TwoNu[0];
821 
822  double limit = tnu.ipTwoPhoE;
823 
824  for( long i=0; i < nCRS; i++ )
825  {
826  fprintf(ioQQQ,"%.3e\t%.3e\n", ener[i] ,
827  atmdat_2phot_shapefunction( ener[i]/0.75, ipISO, nelem ) );
828  }
829 
830  double xnew = 0.;
832  for( long i=0; i < limit; i++ )
833  {
834  double fach = tnu.As2nu[i]/2.*rfield.anu2(i)/rfield.widflx(i)*EN1RYD;
835  fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n",
836  rfield.anu(i) ,
837  tnu.As2nu[i] / rfield.widflx(i) ,
838  fach );
839  xnew += tnu.As2nu[i];
840  }
841  fprintf(ioQQQ," sum is %.3e\n", xnew );
843  }
844  }
845 
846  {
847  enum {DEBUG_LOC=false};
848  if( DEBUG_LOC )
849  {
850  for( long i=0; i<11; ++i )
851  {
852  (*TauDummy).WLAng() = (realnum)(PI * exp10((double)i));
853  fprintf(ioQQQ,"%.2f\t%s\n", (*TauDummy).WLAng() , chLineLbl(*TauDummy).c_str());
854  }
856  }
857  }
858 
859  /* option to print out whole thing with "trace lines" command */
860  if( trace.lgTrLine )
861  {
862  fprintf( ioQQQ, " WL(Ang) E(RYD) IP gl gu gf A damp abs K\n" );
863 
864  /*Atomic Or Molecular lines*/
865  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
866  {
867  for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
868  em != dBaseTrans[ipSpecies].Emis().end(); ++em)
869  {
870  long iWL_Ang = (long)(*em).Tran().WLAng();
871 
872  if( iWL_Ang > 1000000 )
873  {
874  iWL_Ang /= 10000;
875  }
876  else if( iWL_Ang > 10000 )
877  {
878  iWL_Ang /= 1000;
879  }
880  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
881  chLineLbl((*em).Tran()).c_str(), iWL_Ang, RYDLAM/(*em).Tran().WLAng(),
882  (*em).Tran().ipCont(), (long)((*(*em).Tran().Lo()).g()),
883  (long)((*(*em).Tran().Hi()).g()),(*em).gf(),
884  (*em).Aul(),(*em).dampXvel(),
885  (*em).opacity());
886  }
887  }
888 
889  for( long i=0; i < nWindLine; i++ )
890  {
891  long iWL_Ang = (long)TauLine2[i].WLAng();
892 
893  if( iWL_Ang > 1000000 )
894  {
895  iWL_Ang /= 10000;
896  }
897  else if( iWL_Ang > 10000 )
898  {
899  iWL_Ang /= 1000;
900  }
901  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
902  chLineLbl(TauLine2[i]).c_str(), iWL_Ang, RYDLAM/TauLine2[i].WLAng(),
903  TauLine2[i].ipCont(), (long)((*TauLine2[i].Lo()).g()),
904  (long)((*TauLine2[i].Hi()).g()), TauLine2[i].Emis().gf(),
905  TauLine2[i].Emis().Aul(), TauLine2[i].Emis().dampXvel(),
906  TauLine2[i].Emis().opacity() );
907  }
908  for( size_t i=0; i < HFLines.size(); i++ )
909  {
910  long iWL_Ang = (long)HFLines[i].WLAng();
911 
912  if( iWL_Ang > 1000000 )
913  {
914  iWL_Ang /= 10000;
915  }
916  else if( iWL_Ang > 10000 )
917  {
918  iWL_Ang /= 1000;
919  }
920  fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n",
921  chLineLbl(HFLines[i]).c_str(), iWL_Ang, RYDLAM/HFLines[i].WLAng(),
922  HFLines[i].ipCont(), (long)((*HFLines[i].Lo()).g()),
923  (long)((*HFLines[i].Hi()).g()), HFLines[i].Emis().gf(),
924  HFLines[i].Emis().Aul(), HFLines[i].Emis().dampXvel(),
925  HFLines[i].Emis().opacity() );
926  }
927  }
928 
929  /* this is an option to kill fine structure line optical depths */
930  if( !rt.lgFstOn )
931  {
932  /*Atomic or Molecular Lines-Humeshkar Nemala*/
933  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
934  {
935  for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
936  em != dBaseTrans[ipSpecies].Emis().end(); ++em)
937  {
938  if((*em).Tran().EnergyWN() < 10000. )
939  {
940  (*em).opacity() = 0.;
941  }
942  }
943  }
944  }
945 
946  /* read in continuum bands data set */
947  ContBandsCreate( "" );
948 
949  /* we're done adding lines and states to the stacks.
950  * This flag is used to make sure we never add them again in this coreload. */
951  lgLinesAdded = true;
952  lgStatesAdded = true;
953 
955 
956  return;
957 }
958 
959 /*ipShells assign continuum energy pointers to shells for all atoms,
960  * called by ContCreatePointers */
962  /* nelem is the atomic number on the C scale, Li is 2 */
963  long int nelem)
964 {
965  long int
966  imax,
967  ion,
968  nelec,
969  ns,
970  nshell;
971  /* following value cannot be used - will be set to proper threshold */
972  double thresh=-DBL_MAX;
973 
974  DEBUG_ENTRY( "ipShells()" );
975 
976  ASSERT( nelem >= NISO);
977  ASSERT( nelem < LIMELM );
978 
979  /* fills in pointers to valence shell ionization threshold
980  * PH1(a,b,c,d)
981  * a=1 => thresh, others fitting parameters
982  * b atomic number
983  * c number of electrons
984  * d shell number 7-1 */
985 
986  /* threshold in Ryd
987  * ion=0 for atom, up to nelem-1 for helium like, hydrogenic is elsewhere */
988  for( ion=0; ion < nelem; ion++ )
989  {
990  string chLab = chIonLbl( nelem+1, ion+1 );
991 
992  /* this is the iso sequence - must not redo sequence if done as iso */
993  long int ipISO = nelem-ion;
994 
995  /* number of bound electrons */
996  nelec = ipISO+1;
997 
998  /* nsShells(nelem,ion) is the number of shells for ion with nelec electrons,
999  * physical not c scale */
1000  imax = Heavy.nsShells[nelem][ion];
1001 
1002  /* loop on all inner shells, valence shell */
1003  for( nshell=0; nshell < imax; nshell++ )
1004  {
1005  /* ionization potential of this shell in rydbergs */
1006  thresh = (double)(t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
1007  if( thresh <= 0.1 )
1008  {
1009  /* negative ip shell does not exist, set upper limit
1010  * to less than lower limit so this never looped upon
1011  * these are used as flags by LimitSh to check whether
1012  * this is a real shell - if 1 or 2 is changed - change LimitSh!! */
1013  opac.ipElement[nelem][ion][nshell][0] = 2;
1014  opac.ipElement[nelem][ion][nshell][1] = 1;
1015  }
1016  else
1017  {
1018  /* this is lower bound to energy range for this shell */
1019  /* >>chng 02 may 27, change to version of ip with label, so that
1020  * inner shell edges will appear */
1021  /*opac.ipElement[nelem][ion][nshell][0] = ipoint(thresh);*/
1022  opac.ipElement[nelem][ion][nshell][0] =
1023  ipContEnergy( thresh , chLab.c_str() );
1024 
1025  /* this is upper bound to energy range for this shell
1026  * LimitSh is an integer function, returns pointer
1027  * to threshold of next major shell. For k-shell it
1028  * returns the values KshellLimit, default=7.35e4
1029  * >>chng 96 sep 26, had been below, result zero cross sec at
1030  * many energies where opacity project did not produce state specific
1031  * cross section */
1032  opac.ipElement[nelem][ion][nshell][1] =
1033  LimitSh(ion+1, nshell+1,nelem+1);
1034  ASSERT( opac.ipElement[nelem][ion][nshell][1] > 0);
1035  }
1036  }
1037 
1038  ASSERT( imax > 0 && imax <= 7 );
1039 
1040  /* this will be index pointing to valence edge */
1041  /* [0] is pointer to threshold in energy array */
1042  opac.ipElement[nelem][ion][imax-1][0] =
1043  ipContEnergy(thresh, chLab.c_str());
1044 
1045  /* pointer to valence electron ionization potential */
1046  Heavy.ipHeavy[nelem][ion] = opac.ipElement[nelem][ion][imax-1][0];
1047  ASSERT( Heavy.ipHeavy[nelem][ion]>0 );
1048 
1049  /* ionization potential of valence shell in Ryd
1050  * thresh was evaluated above, now has last value, the valence shell */
1051  Heavy.Valence_IP_Ryd[nelem][ion] = thresh;
1052 
1053  Heavy.xLyaHeavy[nelem][ion] = 0.;
1054  if( ipISO >= NISO )
1055  {
1056  /* this is set of 3/4 of valence shell IP, this is important
1057  * source of ots deep in cloud */
1058  Heavy.ipLyHeavy[nelem][ion] =
1059  ipLineEnergy(thresh*0.75, chLab.c_str(), 0);
1060 
1061  Heavy.ipBalHeavy[nelem][ion] =
1062  ipLineEnergy(thresh*0.25, chLab.c_str(), 0);
1063  }
1064  else
1065  {
1066  /* do not treat this simple way since done exactly with iso
1067  * sequences */
1068  Heavy.ipLyHeavy[nelem][ion] = -1;
1069  Heavy.ipBalHeavy[nelem][ion] = -1;
1070  }
1071  }
1072 
1073  /* above loop did up to hydrogenic, now do hydrogenic -
1074  * hydrogenic is special since arrays already set up */
1075  Heavy.nsShells[nelem][nelem] = 1;
1076 
1077  /* this is lower limit to range */
1078  /* hydrogenic photoionization set to special hydro array
1079  * this is pointer to threshold energy */
1080  /* this statement is in ContCreatePointers but has not been done when this routine called */
1081  /*iso_sp[ipH_LIKE][ipZ].fb[ipLo].ipIsoLevNIonCon = ipContEnergy(iso_sp[ipH_LIKE][ipZ].fb[ipLo].xIsoLevNIonRyd,chLab);*/
1082  /*opac.ipElement[nelem][nelem][0][0] = iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon;*/
1083  opac.ipElement[nelem][nelem][0][0] = ipoint( t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787 );
1084  ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );
1085 
1086  /* this is the high-energy limit */
1087  opac.ipElement[nelem][nelem][0][1] = continuum.KshellLimit;
1088 
1089  Heavy.ipHeavy[nelem][nelem] = opac.ipElement[nelem][nelem][0][0];
1090 
1091  /* this is for backwards computability with Cambridge code */
1092  if( trace.lgTrace && trace.lgPointBug )
1093  {
1094  for( ion=0; ion < (nelem+1); ion++ )
1095  {
1096  fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n",
1097  nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
1098  , Heavy.nsShells[nelem][ion] );
1099  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1100  {
1101  fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n",
1102  ns+1, Heavy.chShell[ns], rfield.anu(opac.ipElement[nelem][ion][ns][0]-1)*
1103  EVRYD, rfield.anu(opac.ipElement[nelem][ion][ns][1]-1)*EVRYD );
1104  }
1105  }
1106  }
1107  return;
1108 }
1109 
1110 /*LimitSh sets upper energy limit to subshell integrations */
1111 STATIC long LimitSh(long int ion,
1112  long int nshell,
1113  long int nelem)
1114 {
1115  long int LimitSh_v;
1116 
1117  DEBUG_ENTRY( "LimitSh()" );
1118 
1119  /* this routine returns the high-energy limit to the energy range
1120  * for photoionization of a given shell
1121  * */
1122  if( nshell == 1 )
1123  {
1124  /* this limit is high-energy limit to code unless changed with set kshell */
1125  LimitSh_v = continuum.KshellLimit;
1126 
1127  }
1128  else if( nshell == 2 )
1129  {
1130  /* this is 2s shell, upper limit is 1s
1131  * >>chng 96 oct 08, up to high-energy limit
1132  * LimitSh = ipElement(nelem,ion , 1,1)-1 */
1133  LimitSh_v = continuum.KshellLimit;
1134 
1135  }
1136  else if( nshell == 3 )
1137  {
1138  /* this is 2p shell, upper limit is 1s
1139  * >>chng 96 oct 08, up to high-energy limit
1140  * LimitSh = ipElement(nelem,ion , 1,1)-1 */
1141  LimitSh_v = continuum.KshellLimit;
1142 
1143  }
1144  else if( nshell == 4 )
1145  {
1146  /* this is 3s shell, upper limit is 2p
1147  * >>chng 96 oct 08, up to K-shell edge
1148  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1149  LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1150 
1151  }
1152  else if( nshell == 5 )
1153  {
1154  /* this is 3p shell, upper limit is 2p
1155  * >>chng 96 oct 08, up to K-shell edge
1156  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1157  LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1158 
1159  }
1160  else if( nshell == 6 )
1161  {
1162  /* this is 3d shell, upper limit is 2p
1163  * >>chng 96 oct 08, up to K-shell edge
1164  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
1165  LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
1166 
1167  }
1168  else if( nshell == 7 )
1169  {
1170  /* this is 4s shell, upper limit is 3d */
1171  if( opac.ipElement[nelem-1][ion-1][5][0] < 3 )
1172  {
1173  /* this is check for empty shell 6, 3d
1174  * if so then set to 3p instead */
1175  LimitSh_v = opac.ipElement[nelem-1][ion-1][4][0] -
1176  1;
1177  }
1178  else
1179  {
1180  LimitSh_v = opac.ipElement[nelem-1][ion-1][5][0] -
1181  1;
1182  }
1183  /* >>chng 96 sep 26, set upper limit down to 2s */
1184  LimitSh_v = opac.ipElement[nelem-1][ion-1][2][0] - 1;
1185 
1186  }
1187  else
1188  {
1189  fprintf( ioQQQ, " LimitSh cannot handle nshell as large as%4ld\n",
1190  nshell );
1192  }
1193  return( LimitSh_v );
1194 }
1195 
1196 /*ContBandsCreate - read set of continuum bands to enter total emission into line*/
1198  /* chFile is optional filename, if void then use default bands,
1199  * if not void then use file specified,
1200  * return value is 0 for success, 1 for failure */
1201  const char chFile[] )
1202 {
1203  char chLine[FILENAME_PATH_LENGTH_2];
1204  const char* chFilename;
1205  FILE *ioDATA;
1206 
1207  bool lgEOL;
1208  long int i,k;
1209 
1210  /* keep track of whether we have been called - want to be
1211  * called a total of one time */
1212  static bool lgCalled=false;
1213 
1214  DEBUG_ENTRY( "ContBandsCreate()" );
1215 
1216  /* do nothing if second or later call*/
1217  if( lgCalled )
1218  {
1219  /* success */
1220  return;
1221  }
1222  lgCalled = true;
1223 
1224  /* use default filename if void string, else use file specified */
1225  chFilename = ( strlen(chFile) == 0 ) ? "continuum_bands.ini" : chFile;
1226 
1227  /* get continuum band data */
1228  if( trace.lgTrace )
1229  {
1230  fprintf( ioQQQ, " ContBandsCreate opening %s:", chFilename );
1231  }
1232 
1233  ioDATA = open_data( chFilename, "r" );
1234 
1235  /* now count how many bands are in the file */
1236  continuum.nContBand = 0;
1237 
1238  /* first line is a magic number and does not count as a band*/
1239  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1240  {
1241  fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
1243  }
1244  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1245  {
1246  /* we want to count the lines that do not start with #
1247  * since these contain data */
1248  if( chLine[0] != '#')
1249  ++continuum.nContBand;
1250  }
1251 
1252  /* now rewind the file so we can read it a second time*/
1253  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
1254  {
1255  fprintf( ioQQQ, " ContBandsCreate could not rewind %s.\n", chFilename );
1257  }
1258 
1260  continuum.chContBandLabels = (char **)MALLOC(sizeof(char *)*(unsigned)(continuum.nContBand) );
1261  continuum.ipContBandLow = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
1262  continuum.ipContBandHi = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
1263  continuum.BandEdgeCorrLow = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
1264  continuum.BandEdgeCorrHi = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
1265 
1266  /* now make second dim, id wavelength, and lower and upper bounds */
1267  for( i=0; i<continuum.nContBand; ++i )
1268  {
1269  /* array of labels, each 4 long plus 0 at [4] */
1270  continuum.chContBandLabels[i] = (char *)MALLOC(sizeof(char)*(unsigned)(5) );
1271  }
1272 
1273  /* first line is a versioning magic number - now confirm that it is valid */
1274  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1275  {
1276  fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
1278  }
1279  /* bands_continuum magic number here <- this string is in band_continuum.dat
1280  * with comment to search for this to find magic number */
1281  {
1282  long int m1 , m2 , m3,
1283  // the magic number at the start of the data file
1284  myr = 11, mmo = 9, mdy = 10;
1285 
1286  i = 1;
1287  m1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1288  m2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1289  m3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1290  if( ( m1 != myr ) ||
1291  ( m2 != mmo ) ||
1292  ( m3 != mdy ) )
1293  {
1294  fprintf( ioQQQ,
1295  " ContBandsCreate: the version of the data file %s I found (%li %li %li)is not the current version (%li %li %li).\n",
1296  chFilename ,
1297  m1 , m2 , m3 ,
1298  myr , mmo , mdy );
1299  fprintf( ioQQQ,
1300  " ContBandsCreate: you need to update this file.\n");
1302  }
1303  }
1304 
1305  /* now read in data again, but save it this time */
1306  k = 0;
1307  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1308  {
1309  /* we want to count the lines that do not start with #
1310  * since these contain data */
1311  if( chLine[0] != '#')
1312  {
1313  /* copy 4 char label plus termination */
1314  strncpy( continuum.chContBandLabels[k] , chLine , 4 );
1315  continuum.chContBandLabels[k][4] = 0;
1316 
1317  /* now get central band wavelength
1318  * >>chng 06 aug 11 from 4 to 6, the first 4 char are labels and
1319  * these can contain numbers, next comes a space, then the number */
1320  i = 6;
1321  continuum.ContBandWavelength[k] = (realnum)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1322  /* >>chng 06 feb 21, multiply by 1e4 to convert micron wavelength into Angstroms,
1323  * which is assumed by the code. before this correction the band centroid
1324  * wavelength was given in the output incorrectly listed as Angstroms.
1325  * results were correct just label was wrong */
1326  continuum.ContBandWavelength[k] *= 1e4f;
1327 
1328  /* these are short and long wave limits, which are high and
1329  * low energy limits - these are now wl in microns but are
1330  * converted to Angstroms */
1331  double xHi = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL)*1e4;
1332  double xLow = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL)*1e4;
1333  if( lgEOL )
1334  {
1335  fprintf( ioQQQ, " There should have been 3 numbers on this band line. Sorry.\n" );
1336  fprintf( ioQQQ, " string==%s==\n" ,chLine );
1338  }
1339 
1340  {
1341  enum {DEBUG_LOC=false};
1342  if( DEBUG_LOC )
1343  {
1344  fprintf(ioQQQ, "READ:%s\n", chLine );
1345  fprintf(ioQQQ, "GOT: %s %g %g %g\n",continuum.chContBandLabels[k],
1346  continuum.ContBandWavelength[k] , xHi , xLow );
1347  }
1348  }
1349 
1350  /* make sure bands bounds are in correct order, shorter - longer wavelength*/
1351  if( xHi >= xLow )
1352  {
1353  fprintf( ioQQQ, " ContBandWavelength band %li "
1354  "edges are in improper order.\n" ,k);
1355  fprintf(ioQQQ,"band: %s %.3e %.3e %.3e \n",
1358  xHi ,
1359  xLow);
1361  }
1362 
1363  // check that central wavelength is indeed between the limits
1364  // xHi & xLow are hi and low energy limits to band so logic reversed
1365  if( continuum.ContBandWavelength[k] < xHi ||
1366  continuum.ContBandWavelength[k] > xLow )
1367  {
1368  fprintf( ioQQQ, " ContBandWavelength band %li central "
1369  "wavelength not within band.\n" ,k);
1370  fprintf(ioQQQ,"band: %s %.3e %li %li \n",
1373  continuum.ipContBandHi[k] ,
1376  }
1377 
1378  /* get continuum index - RYDLAM is 911.6A = 1 Ryd so 1e4 converts
1379  * micron to Angstrom - xHi is high energy (not wavelength)
1380  * edge of the band */
1381  continuum.ipContBandHi[k] = ipoint( RYDLAM / xHi );
1382  continuum.ipContBandLow[k] = ipoint( RYDLAM / xLow );
1383 
1384  /* fraction of first and last bin to include */
1387  (realnum)(RYDLAM/xLow)) /
1390  continuum.BandEdgeCorrHi[k] = ((realnum)(RYDLAM/xHi) -
1394  /*fprintf(ioQQQ,"DEBUG bands_continuum %s %.3e %li %li \n",
1395  continuum.chContBandLabels[k],
1396  continuum.ContBandWavelength[k],
1397  continuum.ipContBandHi[k] ,
1398  continuum.ipContBandLow[k]);*/
1399 
1400  if( trace.lgTrace && trace.lgConBug )
1401  {
1402  if( k==0 )
1403  fprintf( ioQQQ, " ContCreatePointer trace bands\n");
1404  fprintf( ioQQQ,
1405  " band %ld label %s low wl= %.3e low ipnt= %li "
1406  " hi wl= %.3e hi ipnt= %li \n",
1407  k,
1409  xLow,
1411  xHi,
1412  continuum.ipContBandHi[k] );
1413  }
1414 # if 0
1415  // hazy table giving band properties
1416 # include "prt.h"
1417  fprintf(ioQQQ,
1418  "DEBUG %s & ",
1421  fprintf(ioQQQ," & ");
1422  prt_wl( ioQQQ , xHi );
1423  fprintf(ioQQQ," -- ");
1424  prt_wl( ioQQQ , xLow );
1425  fprintf(ioQQQ,"\\\\ \n");
1426 # endif
1427  ++k;
1428  }
1429  }
1430  /* now validate this incoming data */
1431  for( i=0; i<continuum.nContBand; ++i )
1432  {
1433  /* make sure all are positive */
1434  if( continuum.ContBandWavelength[i] <=0. )
1435  {
1436  fprintf( ioQQQ, " ContBandWavelength band %li has non-positive entry.\n",i );
1438  }
1439  }
1440 
1441  fclose(ioDATA);
1442  return;
1443 }
long int iphmin
Definition: hmi.h:128
double emm() const
Definition: mesh.h:84
long int & ipFine() const
Definition: emission.h:448
static bool lgCalled
Definition: cddrive.cpp:429
realnum * fine_anu
Definition: rfield.h:393
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:751
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:222
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:230
size_t size(void) const
Definition: transition.h:331
long int ipG0_spec_hi
Definition: rfield.h:253
TransitionList UTALines("UTALines",&AnonStates)
double exp10(double x)
Definition: cddefines.h:1383
long int iheh1
Definition: hmi.h:69
const int ipHE_LIKE
Definition: iso.h:65
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:296
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
realnum EnerGammaRay
Definition: rfield.h:446
long int ipG0_DB96_hi
Definition: rfield.h:250
char chIonStage[LIMELM+1][CHARS_ION_STAGE]
Definition: elementnames.h:29
double o3exc
Definition: opacity.h:295
double widflx(size_t i) const
Definition: mesh.h:147
t_opac opac
Definition: opacity.cpp:5
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:34
t_Heavy Heavy
Definition: heavy.cpp:5
t_fe fe
Definition: fe.cpp:5
const double WL_V_FILT
Definition: rfield.h:13
double abscf(double gf, double enercm, double gl)
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat_adfa.h:61
vector< string > chContLabel
Definition: rfield.h:215
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:310
void set_ipoint(long n, long val)
Definition: yield.h:73
long int ipo1exc[3]
Definition: opacity.h:222
long ipmgex
Definition: opacity.h:299
long int * ipContBandLow
Definition: continuum.h:88
double eina(double gf, double enercm, double gup)
const int ipHe2s3S
Definition: iso.h:46
const int ipOXYGEN
Definition: cddefines.h:355
realnum * BandEdgeCorrHi
Definition: continuum.h:91
long int KshellLimit
Definition: continuum.h:98
long int ipG0_spec_lo
Definition: rfield.h:253
realnum xLyaHeavy[LIMELM][LIMELM]
Definition: heavy.h:21
long ipTwoPhoE
Definition: two_photon.h:41
TransitionList HFLines("HFLines",&AnonStates)
long int ipG0_TH85_hi
Definition: rfield.h:246
long int i2d
Definition: oxy.h:38
long int ipo3exc[3]
Definition: opacity.h:222
long ipFineCont(double energy_ryd)
long int ipEnerGammaRay
Definition: rfield.h:447
FILE * ioQQQ
Definition: cddefines.cpp:7
long int ipG0_TH85_lo
Definition: rfield.h:246
realnum pfe11a
Definition: fe.h:19
long int ipfe10
Definition: fe.h:18
const int ipHe1s1S
Definition: iso.h:43
vector< freeBound > fb
Definition: iso.h:481
TransitionList TauLine2("TauLine2",&AnonStates)
double anu(size_t i) const
Definition: mesh.h:111
long int nSpecies
Definition: taulines.cpp:22
double o3exc3
Definition: opacity.h:296
bool lgPointBug
Definition: trace.h:34
t_dense dense
Definition: global.cpp:15
static t_ADfA & Inst()
Definition: cddefines.h:209
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int nflux_with_check
Definition: rfield.h:51
STATIC long LimitSh(long int ion, long int nshell, long int nelem)
realnum SmallA
Definition: iso.h:391
long int ipSecIon
Definition: secondaries.h:55
t_trace trace
Definition: trace.cpp:5
#define MALLOC(exp)
Definition: cddefines.h:556
bool lgLinesAdded
Definition: taulines.cpp:12
t_ionbal ionbal
Definition: ionbal.cpp:8
realnum EnergyAng() const
Definition: transition.h:100
long int ipG0_DB96_lo
Definition: rfield.h:250
vector< two_photon > TwoNu
Definition: iso.h:598
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
bool lgConBug
Definition: trace.h:97
void ContCreatePointers(void)
realnum pfe10
Definition: fe.h:19
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
long ipContEnergy(double energy, const char *chLabel)
Definition: cont_ipoint.cpp:31
realnum pfe14
Definition: fe.h:19
double energy(const genericState &gs)
const int ipH1s
Definition: iso.h:29
void checkTransitionListOfLists(vector< TransitionList > &list)
Definition: taulines.cpp:39
double & heat() const
Definition: collision.h:224
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
EmissionList::reference Emis() const
Definition: transition.h:447
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
long ipLineEnergy(double energy, const char *chLabel, long ipIonEnergy)
Definition: cont_ipoint.cpp:68
t_continuum continuum
Definition: continuum.cpp:6
STATIC void print_emline_fine(const char *LineGroup, const TransitionProxy &tr)
long int ipCKshell
Definition: opacity.h:310
t_rfield rfield
Definition: rfield.cpp:9
long & ipCont() const
Definition: transition.h:489
bool lgDielRecom[NISO]
Definition: iso.h:385
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
STATIC void ipShells(long int nelem)
vector< diatomics * > diatoms
Definition: h2.cpp:8
bool lgFstOn
Definition: rt.h:187
#define cdEXIT(FAIL)
Definition: cddefines.h:484
realnum * ContBandWavelength
Definition: continuum.h:87
double powi(double, long int)
Definition: service.cpp:690
int nlines() const
Definition: yield.h:76
double anu2(size_t i) const
Definition: mesh.h:115
bool lgTrLine
Definition: trace.h:43
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:35
long nWindLine
Definition: cdinit.cpp:19
realnum EnergyKshell
Definition: continuum.h:99
vector< string > chLineLabel
Definition: rfield.h:212
long int nContBand
Definition: continuum.h:85
double anumin(size_t i) const
Definition: mesh.h:139
bool lgElmtOn[LIMELM]
Definition: dense.h:160
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
string chLabel() const
Definition: transition.cpp:271
long int ipV_filter
Definition: rfield.h:242
const int ipH2p
Definition: iso.h:31
long int ipLyHeavy[LIMELM][LIMELM-1]
Definition: heavy.h:11
long int ipo3exc3[3]
Definition: opacity.h:222
#define ASSERT(exp)
Definition: cddefines.h:617
long int ipBalHeavy[LIMELM][LIMELM-1]
Definition: heavy.h:11
long int ippr
Definition: opacity.h:222
long int ip660
Definition: he.h:21
const int ipH2s
Definition: iso.h:30
const double WL_B_FILT
Definition: rfield.h:16
void TwoPhotonSetup(vector< two_photon > &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem)
Definition: two_photon.cpp:10
t_he he
Definition: he.cpp:5
const int ipH_LIKE
Definition: iso.h:64
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
const int LIMELM
Definition: cddefines.h:307
realnum * BandEdgeCorrLow
Definition: continuum.h:91
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
const int ipHELIUM
Definition: cddefines.h:349
long int in1[3]
Definition: opacity.h:222
void iso_create(void)
Definition: iso_create.cpp:111
vector< qList > dBaseStates
Definition: taulines.cpp:16
long int ip1000A
Definition: rfield.h:256
STATIC void ContBandsCreate(const char chFile[])
vector< species > dBaseSpecies
Definition: taulines.cpp:15
realnum pfe11b
Definition: fe.h:19
double EnergyRyd() const
Definition: transition.h:95
t_oxy oxy
Definition: oxy.cpp:5
const int NCHLAB
Definition: cddefines.h:303
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
vector< realnum > As2nu
Definition: two_photon.h:46
char chShell[7][3]
Definition: heavy.h:31
long int ipB_filter
Definition: rfield.h:242
bool lgStatesAdded
Definition: taulines.cpp:11
vector< TransitionList > AllTransitions
Definition: taulines.cpp:9
long int numLevels_max
Definition: iso.h:524
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:72
long int iheh2
Definition: hmi.h:69
long ica2ex[2]
Definition: opacity.h:299
double anumax(size_t i) const
Definition: mesh.h:143
#define fixit(a)
Definition: cddefines.h:416
long int ih2pnt[2]
Definition: opacity.h:222
t_hmi hmi
Definition: hmi.cpp:5
t_secondaries secondaries
Definition: secondaries.cpp:5
long int ipxry
Definition: rt.h:181
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
long int * ipContBandHi
Definition: continuum.h:88
long int i2p
Definition: oxy.h:38
const int ipHYDROGEN
Definition: cddefines.h:348
realnum & Aul() const
Definition: emission.h:668
long int ih2pnt_ex[2]
Definition: opacity.h:222
double atmdat_2phot_shapefunction(double EbyE2nu, long ipISO, long nelem)
long int nLyman_malloc[NISO]
Definition: iso.h:352
EmissionList & Emis()
Definition: transition.h:363
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
char ** chContBandLabels
Definition: continuum.h:86
long int ** ipCompRecoil
Definition: ionbal.h:156
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:393
t_rt rt
Definition: rt.cpp:5