cloudy  trunk
atmdat_readin.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*atmdat_readin read in some data files, but only if this is very first call,
4  * called by Cloudy */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "taulines.h"
8 #include "mewecoef.h"
9 #include "iterations.h"
10 #include "heavy.h"
11 #include "mole.h"
12 #include "input.h"
13 #include "h2.h"
14 #include "yield.h"
15 #include "trace.h"
16 #include "lines.h"
17 #include "lines_service.h"
18 #include "ionbal.h"
19 #include "struc.h"
20 #include "geometry.h"
21 #include "dense.h"
22 #include "dynamics.h"
23 #include "elementnames.h"
24 #include "hyperfine.h"
25 #include "atmdat.h"
26 #include "iso.h"
27 /* */
28 /* this was needed to get array to crash out of bounds if not set.
29  * std limits on limits.h did not work with visual studio! */
30 const long INTBIG = 2000000000;
31 
32 /* these are the individual pointers to the level 1 lines, they are set to
33  * very large negative.
34  * NB NB NB!!
35  * these occur two times in the code!!
36  * They are declared in taulines.h, and defined here */
37 long ipT1656=INTBIG, ipT9830=INTBIG,
38  ipT8727=INTBIG, ipT1335=INTBIG, ipT1909=INTBIG,
39  ipT977=INTBIG, ipT1550=INTBIG, ipT1548=INTBIG, ipT386=INTBIG,
40  ipT310=INTBIG, ipc31175=INTBIG, ipT291=INTBIG, ipT280=INTBIG,
41  ipT274=INTBIG, ipT270=INTBIG, ipT312=INTBIG, ipT610=INTBIG,
42  ipT370=INTBIG, ipT157=INTBIG, ipT1085=INTBIG,
43  ipT990=INTBIG, ipT1486=INTBIG, ipT765=INTBIG, ipT1243=INTBIG,
44  ipT1239=INTBIG, ipT374g=INTBIG, ipT374x=INTBIG, ipT1200=INTBIG,
45  ipT2140=INTBIG, ipT671=INTBIG, ipT315=INTBIG, ipT324=INTBIG,
46  ipT333=INTBIG, ipT209=INTBIG, ipT122=INTBIG, ipT205=INTBIG,
47  ipT57=INTBIG, ipT6300=INTBIG, ipT6363=INTBIG, ipT5577=INTBIG,
48  ipT834=INTBIG, ipT1661=INTBIG, ipT1666=INTBIG, ipT835=INTBIG,
49  ipT789=INTBIG, ipT630=INTBIG, ipT1304=INTBIG, ipSi10_606=INTBIG,
50  ipT1039=INTBIG, ipT8446=INTBIG, ipT4368=INTBIG, ipTOI13=INTBIG,
51  ipTOI11=INTBIG, ipTOI29=INTBIG, ipTOI46=INTBIG, ipTO1025=INTBIG,
52  ipT304=INTBIG, ipT1214=INTBIG, ipT150=INTBIG, ipT146=INTBIG,
53  ipT63=INTBIG, ipTO88=INTBIG, ipT52=INTBIG, ipT26=INTBIG,
54  ipT1032=INTBIG, ipT1037=INTBIG, ipT770=INTBIG, ipT780=INTBIG,
55  ipxNe0676=INTBIG, ipT895=INTBIG, ipT88=INTBIG, ipTNe13=INTBIG,
56  ipTNe36=INTBIG, ipTNe16=INTBIG, ipTNe14=INTBIG, ipTNe24=INTBIG,
57  ipT5895=INTBIG, ipfsNa373=INTBIG,ipfsNa490=INTBIG, ipfsNa421=INTBIG,
58  ipxNa6143=INTBIG, ipxNa6862=INTBIG,ipxNa0746=INTBIG, ipMgI2853=INTBIG,
59  ipMgI2026=INTBIG, ipT2796=INTBIG, ipT2804=INTBIG,
60  ipT705=INTBIG, ipT4561=INTBIG, ipxMg51325=INTBIG, ipxMg52417=INTBIG,
61  ipxMg52855=INTBIG,ipxMg71190=INTBIG, ipxMg72261=INTBIG,
62  ipxMg72569=INTBIG,ipxMg08303=INTBIG, ipTMg610=INTBIG, ipTMg625=INTBIG,
63  ipT58=INTBIG, ipTMg4=INTBIG, ipTMg14=INTBIG, ipTMg6=INTBIG,
64  ipfsMg790=INTBIG, ipfsMg755=INTBIG, ipAlI3957=INTBIG, ipAlI3090=INTBIG,
65  ipT1855=INTBIG, ipT1863=INTBIG, ipT2670=INTBIG,
66  ipAl529=INTBIG, ipAl6366=INTBIG, ipAl6912=INTBIG, ipAl8575=INTBIG,
67  ipAl8370=INTBIG, ipAl09204=INTBIG, ipT639=INTBIG,
68  ipTAl550=INTBIG, ipTAl568=INTBIG, ipTAl48=INTBIG, ipSii2518=INTBIG,
69  ipSii2215=INTBIG, ipT1808=INTBIG,
70  ipT1207=INTBIG, ipT1895=INTBIG, ipT1394=INTBIG, ipT1403=INTBIG,
71  ipT1527=INTBIG, ipT1305=INTBIG, ipT1260=INTBIG, ipSi619=INTBIG,
72  ipSi10143=INTBIG, ipTSi499=INTBIG, ipTSi521=INTBIG, ipTSi41=INTBIG,
73  ipTSi35=INTBIG, ipTSi25=INTBIG, ipTSi65=INTBIG,
74  ipTSi3=INTBIG, ipTSi4=INTBIG, ipP0260=INTBIG, ipP0233=INTBIG,
75  ipP0318=INTBIG, ipP713=INTBIG, ipT1256=INTBIG, ipT1194=INTBIG,
76  ipTS1720=INTBIG, ipT1198=INTBIG, ipT786=INTBIG,
77  ipTS19=INTBIG, ipTS34=INTBIG, ipTS11=INTBIG, ipfsCl214=INTBIG, ipfsCl233=INTBIG,
78  ipCl04203=INTBIG, ipCl04117=INTBIG, ipCl973=INTBIG,
79  ipTAr7=INTBIG, ipTAr9=INTBIG, ipTAr22=INTBIG, ipTAr13=INTBIG,
80  ipTAr8=INTBIG, ipAr06453=INTBIG, ipKI7745=INTBIG, ipxK03462=INTBIG, ipxK04598=INTBIG,
81  ipxK04154=INTBIG, ipxK07319=INTBIG, ipCaI4228=INTBIG, ipT3934=INTBIG,
82  ipT3969=INTBIG, ipT8498=INTBIG, ipT8542=INTBIG,
83  ipT8662=INTBIG, ipT7291=INTBIG, ipT7324=INTBIG, ipTCa3=INTBIG,
84  ipSc05231=INTBIG, ipSc13264=INTBIG,
85  ipFeI3884=INTBIG, ipFeI3729=INTBIG, ipFeI3457=INTBIG,
86  ipFeI3021=INTBIG, ipFeI2966=INTBIG, ipTuv3=INTBIG,
87  ipTr48=INTBIG, ipTFe16=INTBIG, ipTFe26=INTBIG, ipTFe34=INTBIG,
88  ipTFe35=INTBIG, ipTFe46=INTBIG, ipTFe56=INTBIG, ipT1122=INTBIG,
89  ipT191=INTBIG, ipCo11527=INTBIG;
90 long ipS4_1405=INTBIG,ipS4_1398=INTBIG,ipS4_1424=INTBIG,ipS4_1417=INTBIG,ipS4_1407=INTBIG,
91  ipO4_1400=INTBIG,ipO4_1397=INTBIG,ipO4_1407=INTBIG,ipO4_1405=INTBIG,ipO4_1401=INTBIG,
92  ipN3_1749=INTBIG,ipN3_1747=INTBIG,ipN3_1754=INTBIG,ipN3_1752=INTBIG,ipN3_1751=INTBIG,
93  ipC2_2325=INTBIG,ipC2_2324=INTBIG,ipC2_2329=INTBIG,ipC2_2328=INTBIG,ipC2_2327=INTBIG,
94  ipSi2_2334=INTBIG,ipSi2_2329=INTBIG,ipSi2_2350=INTBIG,ipSi2_2344=INTBIG,ipSi2_2336=INTBIG,
95  ipS1_25m=INTBIG,ipS1_56m=INTBIG,ipCl1_11m=INTBIG,
96  ipFe1_24m=INTBIG,ipFe1_35m=INTBIG,ipFe1_54m=INTBIG,ipFe1_111m=INTBIG,
97  ipNi1_7m=INTBIG ,ipNi1_11m=INTBIG,ipSi1_130m=INTBIG,ipSi1_68m=INTBIG,
98  ipNI_pumpDirect[NI_NDP]={INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG},
99  ipNI_pumpIndirect=INTBIG;
100 
101 /* above are the individual pointers to the level 1 lines, they are set to
102  * very large negative.
103  * NB NB NB!!
104  * these occur two times in the code!!
105  * They are declared in taulines.h, and defined here */
106 
107 /* definition for whether level 2 lines are enabled, will be set to -1
108  * with no level2 command */
109 /*long nWindLine = NWINDDIM;*/
110 /*realnum TauLine2[NWINDDIM][NTA];*/
111 /*realnum **TauLine2;*/
112 #include "atomfeii.h"
113 
114 // NB NB - IS_TOP should always be the last entry of this enum!
115 typedef enum { IS_NONE, IS_K_SHELL, IS_L1_SHELL, IS_L2_SHELL, IS_TOP } exc_type;
116 
117 struct t_BadnellLevel
118 {
119  string config;
120  int irsl;
121  int S;
122  int L;
123  int g; // 2*J+1
124  realnum energy; // in cm^-1
125  bool lgAutoIonizing;
126  exc_type WhichShell;
127  t_BadnellLevel() : irsl(0), S(0), L(0), g(0), energy(0.f), lgAutoIonizing(false), WhichShell(IS_NONE) {}
128 };
129 
130 // read Hummer and Storey 98 He1 photoionization cross-section data
131 STATIC void read_SH98_He1_cross_sections(void);
132 
133 STATIC void read_Helike_cross_sections(void);
134 
135 // read autoionization data from Badnell data file
136 STATIC void ReadBadnellAIData(const string& fnam, // filename containing the Badnell data
137  long nelem, // nelem is on C scale
138  long ion, // ion is on C scale
139  TransitionList& UTA, // UTA lines will be pushed on this stack
140  bitset<IS_TOP> Skip); // option to skip transitions from a particular shell
141 // simple helper functions for ReadBadnellAIData
142 inline void InitTransition(const TransitionProxy& t);
143 inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl);
144 
145 // UTA lines below this absorption oscillator strength value will be ignored -
146 // F+13 paper plotted f not gf
147 const realnum f_cutoff = 1.e-4f;
148 
149 void atmdat_readin(void)
150 {
151  long int i,
152  iCase ,
153  ion,
154  ipDens ,
155  ipISO ,
156  ipTemp ,
157  j,
158  ig0,
159  ig1,
160  imax,
161  nelem,
162  nelec,
163  magic1,
164  magic2,
165  mol;
166 
167  char cha;
168  char chS2[3];
169 
170  long ipZ;
171  bool lgEOL;
172 
173  FILE *ioDATA;
174  char chLine[FILENAME_PATH_LENGTH_2] ,
175  chFilename[FILENAME_PATH_LENGTH_2];
176 
177  static bool lgFirstCall = true;
178 
179  DEBUG_ENTRY( "atmdat_readin()" );
180 
181  /* do nothing if not first call */
182  if( !lgFirstCall )
183  {
184  /* do not do anything, but make sure that number of zones has not increased */
185  bool lgTooBig = false;
186  for( j=0; j < iterations.iter_malloc; j++ )
187  {
188  if( geometry.nend[j]>=struc.nzlim )
189  lgTooBig = true;
190  }
191  if( lgTooBig )
192  {
193  fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
194  fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
195  fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
196  fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
197  fprintf(ioQQQ," Sorry.\n");
198  cdEXIT(EXIT_FAILURE);
199  }
200  return;
201  }
202 
203  lgFirstCall = false; /* do not reevaluate again */
204 
205  /* make sure that molecules have been initialized - this will fail
206  * if this routine is called before size of molecular network is known */
207  if( !mole_global.num_total )
208  {
209  /* mole_global.num_comole_calc can't be zero */
210  TotalInsanity();
211  }
212 
213  /* create space for the structure variables */
214  /* nzlim will be limit, and number allocated */
215  /* >>chng 01 jul 28, define this var, do all following mallocs */
216  for( j=0; j < iterations.iter_malloc; j++ )
217  {
218  struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
219  }
220 
221  /* sloppy, but add one extra for safely */
222  ++struc.nzlim;
223 
224  struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
225 
226  struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
227 
228  struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
229 
230  struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
231 
232  struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
233 
234  struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
235 
236  struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
237 
238  struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
239 
240  struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
241 
242  struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
243 
244  struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
245 
246  struc.AccelTotalOutward = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
247 
248  struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
249 
250  struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
251 
252  struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
253 
254  struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
255 
256  struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
257 
258  struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
259 
260  struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
261 
262  struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
263 
264  struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
265 
266  struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
267 
268  struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
269 
270  if (mole_global.num_calc != 0)
271  {
272  struc.molecules = ((realnum**)MALLOC( (size_t)(mole_global.num_calc)*sizeof(realnum* )));
273  }
274  else
275  {
276  struc.molecules = NULL;
277  }
278 
279  for(mol=0;mol<mole_global.num_calc;mol++)
280  {
281  struc.molecules[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
282  }
283 
284  struc.H2_abund = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
285 
286  /* create space for gas phase abundances array, first create space for the elements */
287  struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)LIMELM );
288  /* last the zones */
289  for( ipZ=0; ipZ< LIMELM;++ipZ )
290  {
291  struc.gas_phase[ipZ] =
292  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
293  }
294 
295  /* create space for struc.xIonDense array, first create space for the zones */
296  struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM) );
297 
298  for( ipZ=0; ipZ<LIMELM;++ipZ )
299  {
300  /* space for the ionization stage */
301  struc.xIonDense[ipZ] =
302  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+1) );
303  /* now create diagonal of space for zones */
304  for( ion=0; ion < (LIMELM+1); ++ion )
305  {
306  struc.xIonDense[ipZ][ion] =
307  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
308  }
309  }
310 
311  struc.StatesElem = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(LIMELM) );
312  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
313  {
314  if( dense.lgElmtOn[nelem] )
315  {
316  struc.StatesElem[nelem] = (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(nelem+1) );
317  for( long ion=0; ion<nelem+1; ion++ )
318  {
319  long ipISO = nelem-ion;
320  if( ipISO < NISO )
321  {
322  struc.StatesElem[nelem][ion] = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
323  for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
324  {
325  struc.StatesElem[nelem][ion][level] =
326  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
327  }
328  }
329  else
330  {
331  fixit(); // for now, point non-iso ions to NULL
332  struc.StatesElem[nelem][ion] = NULL;
333  }
334  }
335  }
336  }
337 
338  /* some structure variables */
339  for( i=0; i < struc.nzlim; i++ )
340  {
341  struc.testr[i] = 0.;
342  struc.volstr[i] = 0.;
343  struc.drad_x_fillfac[i] = 0.;
344  struc.histr[i] = 0.;
345  struc.hiistr[i] = 0.;
346  struc.ednstr[i] = 0.;
347  struc.o3str[i] = 0.;
348  struc.heatstr[i] = 0.;
349  struc.coolstr[i] = 0.;
350  struc.pressure[i] = 0.;
351  struc.pres_radiation_lines_curr[i] = 0.;
352  struc.GasPressure[i] = 0.;
353  struc.DenParticles[i] = 0.;
354  struc.depth[i] = 0.;
355  for(mol=0;mol<mole_global.num_calc;mol++)
356  {
357  struc.molecules[mol][i] = 0.;
358  }
359  struc.H2_abund[i] = 0.;
360  }
361 
362  /* allocate space for some arrays used by dynamics routines, and zero out vars */
363  DynaCreateArrays( );
364 
365  /*************************************************************
366  * *
367  * get the level 1 lines, with real collision data set *
368  * *
369  *************************************************************/
370 
371  if( trace.lgTrace )
372  fprintf( ioQQQ," atmdat_readin reading level1.dat\n");
373 
374  ioDATA = open_data( "level1.dat", "r" );
375 
376  /* first line is a version number and does not count */
377  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
378  {
379  fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
380  cdEXIT(EXIT_FAILURE);
381  }
382 
383  /* count how many lines are in the file, ignoring all lines
384  * starting with '#' */
385  nLevel1 = 0;
386  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
387  {
388  /* we want to count the lines that do not start with #
389  * since these contain data */
390  if( chLine[0] != '#')
391  ++nLevel1;
392  }
393 
394  /* now malloc the TauLines array */
395  TauLines.resize(nLevel1+1);
396  AllTransitions.push_back(TauLines);
397 
398  for( i=0; i<nLevel1+1; i++ )
399  {
400  TauLines[i].Junk();
401  TauLines[i].AddLoState();
402  TauLines[i].AddHiState();
403  TauLines[i].AddLine2Stack();
404  }
405 
406  /* now rewind the file so we can read it a second time*/
407  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
408  {
409  fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
410  cdEXIT(EXIT_FAILURE);
411  }
412 
413  /* check that magic number is ok */
414  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
415  {
416  fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
417  cdEXIT(EXIT_FAILURE);
418  }
419  i = 1;
420  /* level 1 magic number */
421  nelem = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
422  nelec = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
423  ion = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
424 
425  /* magic number
426  * the following is the set of numbers that appear at the start of level1.dat */
427  int nYr=12, nMonth=1, nDay=13;
428  if( ( nelem != nYr ) || ( nelec != nMonth ) || ( ion != nDay ) )
429  {
430  fprintf( ioQQQ,
431  " atmdat_readin: the version of level1.dat is not the current version.\n" );
432  fprintf( ioQQQ,
433  " Please obtain the current version from the Cloudy web site.\n" );
434  fprintf( ioQQQ,
435  " I expected to find the number %i %i %i and got %2.2li %2.2li %2.2li instead.\n" ,
436  nYr, nMonth, nDay,
437  nelem , nelec , ion );
438  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
439  cdEXIT(EXIT_FAILURE);
440  }
441 
442  /* this starts at 1 not 0 since zero is reserved for the
443  * dummy line - we counted number of lines above */
444  i = 1;
445  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
446  {
447  if( i >= (nLevel1+1) )
448  TotalInsanity();
449 
450  /* only look at lines without '#' in first col */
451  if( chLine[0] != '#')
452  {
453  /* first two cols of line has chem element symbol,
454  * try to match it with the list of names
455  * there are no H lines in the level 1 set so zero
456  * is an invalid result */
457 
458  /* copy first two char into chS2 and null terminate */
459  strncpy( chS2 , chLine , 2);
460  chS2[2] = 0;
461 
462  ipZ = 0;
463  for( j=0; j<LIMELM; ++j)
464  {
465  if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
466  {
467  /* ipZ is the actual atomic number starting from 1 */
468  ipZ = j + 1;
469  break;
470  }
471  }
472 
473  /* this happens if no valid chemical element in first two cols on this line */
474  if( ipZ == 0 )
475  {
476  fprintf( ioQQQ,
477  " atmdat_readin could not identify chem symbol on this level 1line:\n");
478  fprintf( ioQQQ, "%s\n", chLine );
479  fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
480  cdEXIT(EXIT_FAILURE);
481  }
482 
483  /* now stuff them into the array, will convert over to proper units later */
484  (*TauLines[i].Hi()).nelem() = (int)ipZ;
485  /* start the scan early on the line -- the element label will not be
486  * picked up, first number will be the ion stage after the label, as in c 4 */
487  j = 1;
488  (*TauLines[i].Hi()).IonStg() = (int)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
489  if( lgEOL )
490  {
491  fprintf( ioQQQ, " There should have been a number on this level1 line 1. Sorry.\n" );
492  fprintf( ioQQQ, "string==%s==\n" ,chLine );
493  cdEXIT(EXIT_FAILURE);
494  }
495 
496  (*TauLines[i].Lo()).nelem() = (*TauLines[i].Hi()).nelem();
497  (*TauLines[i].Lo()).IonStg() = (*TauLines[i].Hi()).IonStg();
498 
499  TauLines[i].WLAng() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
500  if( lgEOL )
501  {
502  fprintf( ioQQQ, " There should have been a number on this level1 line 2. Sorry.\n" );
503  fprintf( ioQQQ, "string==%s==\n" ,chLine );
504  cdEXIT(EXIT_FAILURE);
505  }
506 
507  TauLines[i].EnergyWN() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
508  if( lgEOL )
509  {
510  fprintf( ioQQQ, " There should have been a number on this level1 line 3. Sorry.\n" );
511  fprintf( ioQQQ, "string==%s==\n" ,chLine );
512  cdEXIT(EXIT_FAILURE);
513  }
514 
515  (*TauLines[i].Lo()).g() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
516  if( lgEOL )
517  {
518  fprintf( ioQQQ, " There should have been a number on this level1 line 4. Sorry.\n" );
519  fprintf( ioQQQ, "string==%s==\n" ,chLine );
520  cdEXIT(EXIT_FAILURE);
521  }
522 
523  (*TauLines[i].Hi()).g() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
524  if( lgEOL )
525  {
526  fprintf( ioQQQ, " There should have been a number on this level1 line 5. Sorry.\n" );
527  fprintf( ioQQQ, "string==%s==\n" ,chLine );
528  cdEXIT(EXIT_FAILURE);
529  }
530 
531  /* gf is log if negative */
532  TauLines[i].Emis().gf() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
533  if( lgEOL )
534  {
535  fprintf( ioQQQ, " There should have been a number on this level1 line 6. Sorry.\n" );
536  fprintf( ioQQQ, "string==%s==\n" ,chLine );
537  cdEXIT(EXIT_FAILURE);
538  }
539 
540  if( TauLines[i].Emis().gf() < 0. )
541  TauLines[i].Emis().gf() = (realnum)pow((realnum)10.f,TauLines[i].Emis().gf());
542 
543  /* Emis().Aul() is log if negative */
544  TauLines[i].Emis().Aul() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
545  if( lgEOL )
546  {
547  fprintf( ioQQQ, " There should have been a number on this level1 line 7. Sorry.\n" );
548  fprintf( ioQQQ, "string==%s==\n" ,chLine );
549  cdEXIT(EXIT_FAILURE);
550  }
551  if( TauLines[i].Emis().Aul() < 0. )
552  TauLines[i].Emis().Aul() = (realnum)pow((realnum)10.f,TauLines[i].Emis().Aul());
553 
554  TauLines[i].Emis().iRedisFun() = (int)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
555  if( lgEOL )
556  {
557  fprintf( ioQQQ, " There should have been a number on this level1 line 8. Sorry.\n" );
558  fprintf( ioQQQ, "string==%s==\n" ,chLine );
559  cdEXIT(EXIT_FAILURE);
560  }
561 
562  /* this is new in c94.01 and returns nothing (0.) for most lines */
563  TauLines[i].Coll().col_str() = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
564 
565  /* finally increment i */
566  ++i;
567  }
568  }
569 
570  /* now close the file */
571  fclose(ioDATA);
572  if( trace.lgTrace )
573  fprintf( ioQQQ, " reading level1.dat OK\n" );
574 
575 
576  /*************************************************************
577  * *
578  * get the level 2 line, opacity project, data set *
579  * *
580  *************************************************************/
581 
582  /* nWindLine is initialized to the dimension of the vector when it is
583  * initialized in the definition at the start of this file.
584  * it is set to -1 with the "no level2" command, which
585  * stops us from trying to establish this vector */
586  if( nWindLine > 0 )
587  {
588  /* begin level2 level 2 wind block */
589  /* open file with level 2 line data */
590 
591  /* create the TauLine2 emline array */
592  TauLine2.resize(nWindLine);
593  AllTransitions.push_back(TauLine2);
594  cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum )));
595 
596  /* first initialize entire array to dangerously large negative numbers */
597  for( i=0; i< nWindLine; ++i )
598  {
599  /* >>chng 99 jul 16, from setting all t[] to flt_max, to call for
600  * following, each member of structure set to own type of impossible value */
601  TauLine2[i].Junk();
602 
603  TauLine2[i].AddHiState();
604  TauLine2[i].AddLoState();
605  TauLine2[i].AddLine2Stack();
606  }
607 
608  if( trace.lgTrace )
609  fprintf( ioQQQ," atmdat_readin reading level2.dat\n");
610 
611  ioDATA = open_data( "level2.dat", "r" );
612 
613  /* get magic number off first line */
614  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
615  {
616  fprintf( ioQQQ, " level2.dat error getting magic number\n" );
617  cdEXIT(EXIT_FAILURE);
618  }
619 
620  /* scan off three numbers, should be the yr mn dy of creation date */
621  i = 1;
622  nelem = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
623  nelec = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
624  ion = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
625 
626  /* level2.dat magic number
627  * the following is the set of numbers that appear at the start of level2.dat */
628  if( lgEOL || ( nelem != 9 ) || ( nelec != 11 ) || ( ion != 18 ) )
629  {
630  fprintf( ioQQQ,
631  " atmdat_readin: the version of level2.dat is not the current version.\n" );
632  fprintf( ioQQQ,
633  " I expected to find the number 09 11 18 and got %2.2li %2.2li %2.2li instead.\n" ,
634  nelem , nelec , ion );
635  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
636  fprintf( ioQQQ, "Please obtain the correct version.\n");
637  cdEXIT(EXIT_FAILURE);
638  }
639 
640  /* now get the actual data */
641  i = 0;
642  while( i < nWindLine )
643  {
644  /* this must be double for sscanf to work below */
645  double tt[7];
646  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
647  {
648  fprintf( ioQQQ, " level2.dat error getting line %li\n", i );
649  cdEXIT(EXIT_FAILURE);
650  }
651  /* option to skip any line starting with # */
652  if( chLine[0]!='#' )
653  {
654  /*printf("string is %s\n",chLine );*/
655  sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " ,
656  &tt[0] ,
657  &tt[1] ,
658  &tt[2] ,
659  &tt[3] ,
660  &tt[4] ,
661  &tt[5] ,
662  &tt[6] );
663  /* these are readjusted into their final form in the structure
664  * in routine lines_setup*/
665  (*TauLine2[i].Hi()).nelem() = (int)tt[0];
666  (*TauLine2[i].Hi()).IonStg() = (int)tt[1];
667  (*TauLine2[i].Lo()).g() = (realnum)tt[2];
668  (*TauLine2[i].Hi()).g() = (realnum)tt[3];
669  TauLine2[i].Emis().gf() = (realnum)tt[4];
670  TauLine2[i].EnergyWN() = (realnum)tt[5];
671  cs1_flag_lev2[i] = (realnum)tt[6];
672  ++i;
673  }
674  }
675  /* get magic number off last line */
676  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
677  {
678  fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
679  cdEXIT(EXIT_FAILURE);
680  }
681  sscanf( chLine , "%ld" , &magic2 );
682  if( 999 != magic2 )
683  {
684  fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n",
685  magic2 );
686  cdEXIT(EXIT_FAILURE);
687  }
688  fclose( ioDATA );
689  if( trace.lgTrace )
690  fprintf( ioQQQ," reading level2.dat OK\n");
691 
692  /* end of level2 block*/
693  }
694 
695  /* the UTA line sets - trace will print summary of various data sources */
696 
697  /* reserve space for all data sets, no worries if too small though... */
698  UTALines.reserve( 4000 );
699  AllTransitions.push_back(UTALines);
700 
701  // version of element symbols in lower case and without spaces....
702  const char* chElmSymLC[] =
703  { "h", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al", "si", "p",
704  "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn", "fe", "co", "ni", "cu", "zn" };
705 
706  // save cite for UTA sources, insure no double counting
707  char chUTA_ref[LIMELM][LIMELM][5];
708  for( long nelem=0; nelem < LIMELM; ++nelem )
709  for( long ion=0; ion <= nelem; ++ion )
710  strcpy( chUTA_ref[nelem][ion] , "" );
711 
712  /* first read in the Badnell data */
713  for( ipISO=ipLI_LIKE; ipISO < ipAL_LIKE; ++ipISO )
714  {
715  for( nelem=ipISO; nelem < LIMELM; ++nelem )
716  {
717  // ion = 0 for neutral atom
718  long ion = nelem - ipISO;
719  strcat( chUTA_ref[nelem][ion] , "B" );
720 
721  bitset<IS_TOP> Skip;
722  if( ipISO < ipNA_LIKE )
723  {
724  // construct file name
725  ostringstream oss;
726  // Badnell calls Li-like series He-like, etc...
727  oss << "UTA/nrb00_" << chElmSymLC[ipISO-1] << "_";
728  // Badnell uses ion = 1 for neutral atom
729  oss << chElmSymLC[nelem] << ion+1 << "ic1-2.dat";
730  // now read the data...
731  ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
732  }
733  else
734  {
735  // from Na-like onwards both K-shell (ic1-3) and L-shell (ic2-3)
736  // excitation are treated in two separate files
737  ostringstream oss;
738  oss << "UTA/nrb00_" << chElmSymLC[ipISO-1] << "_";
739  oss << chElmSymLC[nelem] << ion+1 << "ic1-3.dat";
740  ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
741 
742  // Kisielius L1, L2-shell data sets take precedence for Na, Mg, and Al-like iron
743  if( ionbal.lgInnerShell_Kisielius && nelem == ipIRON &&
744  ipISO >= ipNA_LIKE && ipISO <= ipAL_LIKE )
745  {
746  Skip[IS_L1_SHELL] = 1;
747  Skip[IS_L2_SHELL] = 1;
748  }
749 
750  ostringstream oss2;
751  oss2 << "UTA/nrb00_" << chElmSymLC[ipISO-1] << "_";
752  oss2 << chElmSymLC[nelem] << ion+1 << "ic2-3.dat";
753  ReadBadnellAIData( oss2.str(), nelem, ion, UTALines, Skip );
754  }
755  }
756  }
757 
758  /* these are the statistical weights for the ground levels of all ions of iron */
759  const realnum StatWeightGroundLevelIron[] =
760  { 9.f, 10.f, 9.f, 6.f, 1.f, 4.f, 5.f, 4.f, 1.f, 4.f, 5.f, 4.f, 1.f,
761  2.f, 1.f, 2.f, 1.f, 4.f, 5.f, 4.f, 1.f, 2.f, 1.f, 2.f, 1.f, 2.f };
762 
763  // blank line that will be pushed on the UTA line stack
764  qList BlankStates(1);
765  TransitionList BlankList("BlankList",&BlankStates,1);
766  TransitionList::iterator BlankLine = BlankList.begin();
767  (*BlankLine).Junk();
768 
769  /* next read in the Gu file */
770  if( ionbal.lgInnerShell_Gu06 )
771  {
772  /* read the Gu et al. (2006) data
773  * >>refer Fe UTA Gu, M. F., Holczer T., Behar E., & Kahn S. M. 2006, ApJ 641, 1227-1232 */
774  if( trace.lgTrace )
775  fprintf( ioQQQ," atmdat_readin reading UTA_Gu06.dat\n");
776 
777  FILE *ioGU06 = open_data( "UTA/UTA_Gu06.dat", "r" );
778 
779  nelem = 0;
780  /* get up to magic number in Gu 06 file - there is a large header of
781  * comments with the first data the magic number */
782  while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
783  {
784  /* we want to break on first line that does not start with #
785  * since that contains data */
786  if( chLine[0] != '#')
787  break;
788  }
789  /* now get Gu magic number */
790  sscanf( chLine, "%li %li %li", &nelem, &nelec, &ion );
791 
792  /* is magic number correct?
793  * the following is the set of numbers that appear at the start of UTA_Gu06.dat 2006 11 17 */
794  if( nelem != 2007 || nelec != 1 || ion != 23 )
795  {
796  fprintf( ioQQQ,
797  " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
798  fprintf( ioQQQ,
799  " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
800  nelem , nelec , ion );
801  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
802  cdEXIT(EXIT_FAILURE);
803  }
804 
805  int nelemGu =-1, ionGu=-1;
806  while( read_whole_line( chLine, (int)sizeof(chLine), ioGU06 ) != NULL )
807  {
808  if( chLine[0] != '#' )
809  {
810  long int i2;
811  double WLAng, Aul, oscill, Aauto;
812 
813  sscanf( chLine, "%4li%5li%8lf%13lf%12lf",
814  &ion, &i2, &WLAng, &Aul, &Aauto );
815  sscanf( &chLine[54], "%13lf", &oscill );
816 
817  // avoid duplication of ions: anything upto and including the Mg-like
818  // series is covered by the Badnell data set, Al-like is covered by
819  // the Kisielius data set if it is turned on.
820  ipISO = ipIRON - ion + 1;
821  int ipThres = ionbal.lgInnerShell_Kisielius ? ipAL_LIKE : ipMG_LIKE;
822  if( ipISO <= ipThres )
823  continue;
824 
825  UTALines.push_back( *BlankLine );
826  InitTransition( UTALines.back() );
827 
828  /* all these are iron, first number was ion stage with 0 the atom */
829  (*UTALines.back().Hi()).nelem() = ipIRON+1;
830 
831  /* now do stage of ionization */
832  ASSERT( ion > 0 && ion <= ipIRON );
833  /* the ion stage - 1 is atom - this data file has different
834  * format from other two - others are 0 for atom */
835  (*UTALines.back().Hi()).IonStg() = ion;
836  if( ipIRON!=nelemGu || ion!=ionGu )
837  {
838  // one label per ion
839  nelemGu = ipIRON;
840  ionGu = ion;
841  strcat( chUTA_ref[ipIRON][ion-1] , "G" );
842  }
843 
844  /* these are the statistical weights
845  * lower levels are not included in the original data file */
846  if( strstr_s( chLine, "(J=1/2)" ) != NULL )
847  (*UTALines.back().Hi()).g() = 2.f;
848  else if( strstr_s( chLine, "(J=1)" ) != NULL )
849  (*UTALines.back().Hi()).g() = 3.f;
850  else if( strstr_s( chLine, "(J=3/2)" ) != NULL )
851  (*UTALines.back().Hi()).g() = 4.f;
852  else if( strstr_s( chLine, "(J=2)" ) != NULL )
853  (*UTALines.back().Hi()).g() = 5.f;
854  else if( strstr_s( chLine, "(J=5/2)" ) != NULL )
855  (*UTALines.back().Hi()).g() = 6.f;
856  else if( strstr_s( chLine, "(J=3)" ) != NULL )
857  (*UTALines.back().Hi()).g() = 7.f;
858  else if( strstr_s( chLine, "(J=7/2)" ) != NULL )
859  (*UTALines.back().Hi()).g() = 8.f;
860  else if( strstr_s( chLine, "(J=4)" ) != NULL )
861  (*UTALines.back().Hi()).g() = 9.f;
862  else if( strstr_s( chLine, "(J=9/2)" ) != NULL )
863  (*UTALines.back().Hi()).g() = 10.f;
864  else if( strstr_s( chLine, "(J=5)" ) != NULL )
865  (*UTALines.back().Hi()).g() = 11.f;
866  else if( strstr_s( chLine, "(J=11/2)" ) != NULL )
867  (*UTALines.back().Hi()).g() = 12.f;
868  else
869  TotalInsanity();
870  (*UTALines.back().Lo()).g() = StatWeightGroundLevelIron[ion-1];
871 
872  /* wavelength in Angstroms */
873  UTALines.back().WLAng() = (realnum)WLAng;
874  UTALines.back().EnergyWN() = (realnum)(1e8/WLAng);
875 
876  /* store branching ratio for autoionization */
877  double frac_ioniz = Aauto/(Aul + Aauto);
878  ASSERT( frac_ioniz >= 0. && frac_ioniz <= 1. );
879  UTALines.back().Emis().AutoIonizFrac() = (realnum)frac_ioniz;
880 
881  /* save gf scanned from line */
882  UTALines.back().Emis().gf() = (*UTALines.back().Lo()).g() * (realnum)oscill;
883 
884  /* this is true spontaneous rate for doubly excited state to inner
885  * shell UTA, and is used for pumping, and also relaxing back to inner shell */
886  UTALines.back().Emis().Aul() =
887  (realnum)eina( UTALines.back().Emis().gf(),
888  UTALines.back().EnergyWN(), (*UTALines.back().Hi()).g() );
889 
890  ASSERT( fp_equal_tol( (realnum)Aul, UTALines.back().Emis().Aul(), 1.e-3f*(realnum)Aul ) );
891 
892  UTALines.back().Emis().iRedisFun() = ipPRD;
893 
894  UTALines.back().Emis().dampXvel() = (realnum)(
895  (UTALines.back().Emis().Aul()+Aauto) /
896  UTALines.back().EnergyWN()/PI4);
897  ASSERT( UTALines.back().Emis().dampXvel()>0. );
898 
899  // ignore line if too weak
900  if( UTALines.back().Emis().gf() < StatWeightGroundLevelIron[ion-1] * f_cutoff )
901  UTALines.pop_back();
902  }
903  }
904 
905  fclose( ioGU06 );
906 
907  if( trace.lgTrace )
908  fprintf( ioQQQ, " reading UTA_Gu06.dat OK\n" );
909  }
910  else
911  {
912  /* this branch, get the Behar 01 data */
913  /* >>refer Fe UTA Behar, E., Sako, M., & Kahn, S. M. 2001, ApJ, 563, 497-504 */
914  if( trace.lgTrace )
915  fprintf( ioQQQ," atmdat_readin reading UTA_Behar.dat\n");
916 
917  FILE *ioBEHAR = open_data( "UTA/UTA_Behar.dat", "r" );
918 
919  /* UTA_Behar.dat magic number */
920  nelem = 0;
921  if( read_whole_line( chLine, (int)sizeof(chLine), ioBEHAR ) != NULL )
922  sscanf( chLine, "%li %li %li", &nelem, &nelec, &ion );
923 
924  /* magic number
925  * the following is the set of numbers that appear at the start of UTA_Behar.dat 2002 8 19 */
926  if( nelem != 2002 || nelec != 10 || ion != 22 )
927  {
928  fprintf( ioQQQ,
929  " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" );
930  fprintf( ioQQQ,
931  " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
932  nelem , nelec , ion );
933  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
934  cdEXIT(EXIT_FAILURE);
935  }
936 
937  /* read in the Behar 01 data */
938  int nelemBehar=-1, ionBehar=-1;
939  while( read_whole_line( chLine, (int)sizeof(chLine), ioBEHAR ) != NULL )
940  {
941  if( chLine[0] != '#' )
942  {
943  long int i1, i2, i3;
944  double f1, f2, oscill;
945  double frac_relax;
946 
947  sscanf( chLine ,"%li\t%li\t%li\t%lf\t%lf\t%lf\t%lf",
948  &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
949 
950  // skip line if too weak
951  if( oscill < f_cutoff )
952  continue;
953 
954  // avoid duplication of ions: anything upto and including the Mg-like
955  // series is covered by the Badnell data set, Al-like is covered by
956  // the Kisielius data set if it is turned on.
957  ipISO = ipIRON - i1;
958  int ipThres = ionbal.lgInnerShell_Kisielius ? ipAL_LIKE : ipMG_LIKE;
959  if( ipISO <= ipThres )
960  continue;
961 
962  UTALines.push_back( *BlankLine );
963  InitTransition( UTALines.back());
964 
965  /* all these are iron, first number was ion stage with 0 the atom */
966  (*UTALines.back().Hi()).nelem() = ipIRON+1;
967 
968  /* now do stage of ionization */
969  ASSERT( i1 >= 0 && i1 < ipIRON );
970  /* NB - the plus one is because 1 will be subtracted when used,
971  * in the original data file i1 is 0 for the atom */
972  (*UTALines.back().Hi()).IonStg() = i1 + 1;
973  if( ipIRON!=nelemBehar || i1!=ionBehar )
974  {
975  // one label per ion
976  nelemBehar = ipIRON;
977  ionBehar = i1;
978  strcat( chUTA_ref[ipIRON][i1] , "b" );
979  }
980 
981  /* wavelength in Angstroms */
982  UTALines.back().WLAng() = (realnum)f1;
983  UTALines.back().EnergyWN() = 1.e8f/(realnum)f1;
984 
985  (*UTALines.back().Lo()).g() = StatWeightGroundLevelIron[i1];
986  /* now derive g_up by comparing Aul (is f2) and the oscillator strength
987  * this procedure will cause flu_test to be a factor g_up too small */
988  double flu_test = GetGF( f2, UTALines.back().EnergyWN(), 1. )/(*UTALines.back().Lo()).g();
989  /* so dividing the real oscillator strength by flu_test yields g_up... */
990  (*UTALines.back().Hi()).g() = (realnum)nint( oscill/flu_test );
991 
992  ASSERT( i2 == 0 || fp_equal( (*UTALines.back().Hi()).g(), (realnum)i2 ) );
993  ASSERT( i3 == 0 || fp_equal( (*UTALines.back().Lo()).g(), (realnum)i3 ) );
994 
995  /* this is true spontaneous rate for doubly excited state to inner shell UTA,
996  * and is used for pumping, and also relaxing back to inner shell */
997  UTALines.back().Emis().Aul() = (realnum)f2;
998 
999  UTALines.back().Emis().gf() = (*UTALines.back().Lo()).g() * (realnum)oscill;
1000 
1001  UTALines.back().Emis().iRedisFun() = ipPRD;
1002 
1003  UTALines.back().Emis().dampXvel() = (realnum)(
1004  (UTALines.back().Emis().Aul()/frac_relax) /
1005  UTALines.back().EnergyWN()/PI4);
1006  ASSERT( UTALines.back().Emis().dampXvel()>0. );
1007 
1008  /* store branching ratio for autoionization */
1009  ASSERT( frac_relax >= 0.f && frac_relax <= 1.f );
1010  UTALines.back().Emis().AutoIonizFrac() = 1.f-(realnum)frac_relax;
1011  }
1012  }
1013 
1014  fclose( ioBEHAR );
1015 
1016  if( trace.lgTrace )
1017  fprintf( ioQQQ, " reading UTA_Behar.dat OK\n" );
1018  }
1019 
1020  if( ionbal.lgInnerShell_Kisielius )
1021  {
1022  /* last read in the Romas Kisielius data
1023  *>>refer Fe UTA Kisielius, R., Hibbert, A.. Ferland, G. J., et al. 2003, MNRAS, 344, 696 */
1024  if( trace.lgTrace )
1025  fprintf( ioQQQ," atmdat_readin reading UTA_Kisielius.dat\n");
1026 
1027  FILE *ioROMAS = open_data( "UTA/UTA_Kisielius.dat", "r" );
1028 
1029  // magic number
1030  while( read_whole_line( chLine , (int)sizeof(chLine) , ioROMAS ) != NULL )
1031  {
1032  // skip any #
1033  if( chLine[0] != '#')
1034  break;
1035  }
1036  sscanf( chLine, "%li %li %li", &nelem, &nelec, &ion );
1037  if( nelem != 11 || nelec != 8 || ion != 25 )
1038  {
1039  fprintf( ioQQQ,
1040  " atmdat_readin: the version of UTA_Kisielius.dat is not the current version.\n" );
1041  fprintf( ioQQQ,
1042  " I expected to find the number 11 8 25 and got %li %li %li instead.\n" ,
1043  nelem , nelec , ion );
1044  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1045  cdEXIT(EXIT_FAILURE);
1046  }
1047 
1048  long int nRomasUsed = 0 , nRomasTotal = 0;
1049  int nelemRomas=-1 , ionRomas=-1;
1050  FILE *ioROMASused;
1051  bool lgSaveRomasUsed = false;
1052  if( lgSaveRomasUsed )
1053  {
1054  if( (ioROMASused=fopen("RomasUsed.txt","w"))==NULL)
1055  {
1056  fprintf(ioQQQ,"could not open RomasUsed.txt\n");
1057  cdEXIT(EXIT_FAILURE);
1058  }
1059  }
1060  while( read_whole_line( chLine, (int)sizeof(chLine), ioROMAS ) != NULL )
1061  {
1062  /* only look at lines without '#' in first col */
1063  if( chLine[0] != '#' )
1064  {
1065  long int i1, i2, i3;
1066  double f1, f2, oscill;
1067  double frac_relax;
1068 
1069  ++nRomasTotal;
1070  sscanf( chLine, "%li\t%li\t%li\t%lf\t%lf\t%lf\t%lf",
1071  &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
1072 
1073  // skip line if too weak
1074  // Fe+13 has 25 000 lines, most are vastly weak and do not add to integrated rate
1075  // the cutoff of 1e-4 reduces the number of Romas UTA lines from 84224 to 573
1076  if( oscill < f_cutoff )
1077  continue;
1078 
1079  if( lgSaveRomasUsed )
1080  fprintf(ioROMASused , "%s" , chLine);
1081 
1082  /* For Fe XIV both levels of the 2P* ground term are present in the data.
1083  * following true, ignore the transitions from the excited level
1084  * false, assume ground term populated by statistical weight if
1085  * "fudge 0" also appears in input stream */
1086  const bool lgAllowSplitFe14 = false;
1087  if( lgAllowSplitFe14 || i2 == StatWeightGroundLevelIron[i1] )
1088  {
1089  ++nRomasUsed;
1090  UTALines.push_back( *BlankLine );
1091  InitTransition( UTALines.back());
1092 
1093  /* all these are iron, first number was ion stage with 0 the atom */
1094  (*UTALines.back().Hi()).nelem() = ipIRON+1;
1095 
1096  /* now do stage of ionization */
1097  ASSERT( i1 >= 0 && i1 < ipIRON );
1098  /* NB - the plus one is because 1 will be subtracted when used,
1099  * in the original data file i1 is 0 for the atom */
1100  (*UTALines.back().Hi()).IonStg() = i1 + 1;
1101 
1102  realnum facpop = 1.;
1103  // allow low or high density limit for level populations in ground term of Fe XIV
1104  if( i1 == 13 )
1105  {
1106  // Fe XIV - fudge -1 returns number of fudge parameters, 0 if not specified
1107  if( lgAllowSplitFe14 && fudge(-1) )
1108  {
1109  if( i2 == StatWeightGroundLevelIron[i1] )
1110  facpop= 0.333;
1111  else
1112  facpop = 0.667;
1113  }
1114  else
1115  {
1116  if( i2 == StatWeightGroundLevelIron[i1] )
1117  facpop= 1.;
1118  else
1119  facpop = 1e-10;
1120  }
1121  }
1122  if( ipIRON!=nelemRomas || i1!=ionRomas )
1123  {
1124  // one label per ion
1125  nelemRomas = ipIRON;
1126  ionRomas = i1;
1127  strcat( chUTA_ref[ipIRON][i1] , "K" );
1128  }
1129 
1130  /* these were the statistical weights */
1131  (*UTALines.back().Hi()).g() = (realnum)i3;
1132  (*UTALines.back().Lo()).g() = (realnum)i2;
1133 
1134  UTALines.back().WLAng() = (realnum)f1;
1135  UTALines.back().EnergyWN() = 1.e8f/(realnum)f1;
1136  // print transitions contributing to 15.5A UTA feature
1137 # if 0
1138  if( i1==13 && f1>15.35 && f1<15.55)
1139  {
1140  fprintf(ioQQQ,"DEBUG %li\t%.5f\t%.3e\n",i2, f1 , oscill * facpop);
1141  }
1142 # endif
1143 
1144  /* this is true spontaneous rate for doubly excited state to inner shell,
1145  * and is used for pumping, and also relaxing back to inner shell */
1146  UTALines.back().Emis().gf() = (*UTALines.back().Lo()).g() * (realnum)oscill*facpop;
1147  UTALines.back().Emis().Aul() =
1148  (realnum)eina( UTALines.back().Emis().gf(),
1149  UTALines.back().EnergyWN(),
1150  (*UTALines.back().Hi()).g() );
1151 
1152  UTALines.back().Emis().iRedisFun() = ipPRD;
1153 
1154  /* store branching ratio for autoionization */
1155  ASSERT( frac_relax >= 0.f && frac_relax <= 1.f );
1156  UTALines.back().Emis().AutoIonizFrac() = 1.f-(realnum)frac_relax;
1157 
1158  // Romas data do not have autoionization rates so take typical
1159  // value of 1e15 s-1, suggested by Badnell OP data
1160  UTALines.back().Emis().dampXvel() = (realnum)(
1161  (UTALines.back().Emis().Aul()+1e15) /
1162  UTALines.back().EnergyWN()/PI4);
1163  ASSERT( UTALines.back().Emis().dampXvel()>0. );
1164  //fprintf(ioQQQ,"DEBUGGG %li %.3e\n", nRomasUsed, UTALines.back().Emis().dampXvel() );
1165  }
1166  }
1167  }
1168 
1169  fclose( ioROMAS );
1170  if( lgSaveRomasUsed )
1171  fclose( ioROMASused );
1172 
1173  if( trace.lgTrace )
1174  fprintf( ioQQQ, " reading UTA_Kisielius.dat OK,used %li lines from a total of %li\n" , nRomasUsed , nRomasTotal );
1175  }
1176 
1177  nUTA = (long)UTALines.size();
1178 
1179  /* option to dump UTA lines */
1180  if( trace.lgTrace )
1181  {
1182  fprintf(ioQQQ,"\nUTA data sources; B=Badnell 05; G==Gu 06, b=Behar, K=2011 paper\n");
1183  fprintf(ioQQQ," ion ");
1184  for( long ion=0; ion<=LIMELM; ++ion )
1185  fprintf(ioQQQ,"%4li",ion);
1186  fprintf(ioQQQ,"\n");
1187  for( long nelem=0; nelem<LIMELM; ++nelem )
1188  {
1189  fprintf(ioQQQ,"%4s ", elementnames.chElementSym[nelem] );
1190  for( long ion=0; ion<=nelem; ++ion )
1191  {
1192  fprintf(ioQQQ,"%4s",chUTA_ref[nelem][ion] );
1193  }
1194  fprintf(ioQQQ,"\n");
1195  }
1196  fprintf(ioQQQ," ion ");
1197  for( long ion=0; ion<=LIMELM; ++ion )
1198  fprintf(ioQQQ,"%4li",ion);
1199  fprintf(ioQQQ,"\n\n");
1200  }
1201 
1202  if( false )
1203  {
1204  for( i=0; i < nUTA; ++i )
1205  dprintf( ioQQQ, "%5ld %s %2ld wavl %7.3f glo %2g gup %2g Aul %.2e gf %.2e ai branch %.3f\n",
1206  i,
1207  elementnames.chElementSym[(*UTALines[i].Hi()).nelem()-1],
1208  (*UTALines[i].Hi()).IonStg(),
1209  UTALines[i].WLAng(),
1210  (*UTALines[i].Lo()).g(),
1211  (*UTALines[i].Hi()).g(),
1212  UTALines[i].Emis().Aul(),
1213  UTALines[i].Emis().gf(),
1214  UTALines[i].Emis().AutoIonizFrac() );
1215  cdEXIT(EXIT_SUCCESS);
1216  }
1217 
1218  /* end UTA */
1219 
1220  /* read in data for the set of hyperfine structure lines, and allocate
1221  * space for the transition HFLines[nHFLines] structure */
1222  HyperfineCreate();
1223 
1224  /* Make sure that if hybrid is on, then Stout/Chianti are on */
1225  if( atmdat.lgChiantiHybrid && !atmdat.lgChiantiOn)
1226  {
1227  TotalInsanity();
1228  }
1229  if( atmdat.lgStoutHybrid && !atmdat.lgStoutOn )
1230  {
1231  TotalInsanity();
1232  }
1233 
1234  /* read in atomic and molecular models from third-party databases */
1235  if( atmdat.lgLamdaOn || atmdat.lgChiantiOn || atmdat.lgStoutOn)
1236  database_readin();
1237  else
1238  nSpecies = 0;
1239 
1240  /* initialize the large block of level 1 real lines, and OP level 2 lines */
1241  lines_setup();
1242 
1243  /* mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat ========*/
1244  /* read in g-bar data taken from
1245  *>>refer all gbar Mewe, R., Gronenschild, E. H. B. M., van den Oord, G. H. J. 1985, A&AS, 62, 197 */
1246  /* open file with Mewe coefficients */
1247 
1248  if( trace.lgTrace )
1249  fprintf( ioQQQ," atmdat_readin reading mewe_gbar.dat\n");
1250 
1251  ioDATA = open_data( "mewe_gbar.dat", "r" );
1252 
1253  /* get magic number off first line */
1254  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1255  {
1256  fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" );
1257  cdEXIT(EXIT_FAILURE);
1258  }
1259  /* check the number */
1260  sscanf( chLine , "%ld" , &magic1 );
1261  if( magic1 != 9101 )
1262  {
1263  fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n",
1264  magic1 );
1265  cdEXIT(EXIT_FAILURE);
1266  }
1267 
1268  /* now get the actual data, indices are correct for c, in Fort went from 2 to 210 */
1269  for( i=1; i < 210; i++ )
1270  {
1271  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1272  {
1273  fprintf( ioQQQ, " mewe_gbar.dat error getting line %li\n", i );
1274  cdEXIT(EXIT_FAILURE);
1275  }
1276  /*printf("%s\n",chLine);*/
1277  double help[4];
1278  sscanf( chLine, "%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] );
1279  for( int l=0; l < 4; ++l )
1280  MeweCoef.g[i][l] = (realnum)help[l];
1281  }
1282 
1283  /* get magic number off last line */
1284  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1285  {
1286  fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
1287  cdEXIT(EXIT_FAILURE);
1288  }
1289 
1290  sscanf( chLine , "%ld" , &magic2 );
1291 
1292  if( magic1 != magic2 )
1293  {
1294  fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n",
1295  magic2 );
1296  cdEXIT(EXIT_FAILURE);
1297  }
1298 
1299  fclose( ioDATA );
1300 
1301  if( trace.lgTrace )
1302  fprintf( ioQQQ," reading mewe_gbar.dat OK \n");
1303 
1304  /* This is what remains of the t_yield initialization
1305  * this should not be in the constructor of t_yield
1306  * since it initializes a different struct */
1307  for( nelem=0; nelem < LIMELM; nelem++ )
1308  for( ion=0; ion < LIMELM; ion++ )
1309  Heavy.nsShells[nelem][ion] = LONG_MAX;
1310 
1311  /* now read in all auger yields
1312  * will do elements from li on up,
1313  * skip any line starting with #
1314  * this loop goes from lithium to Zn */
1315  for( nelem=2; nelem < LIMELM; nelem++ )
1316  {
1317  /* nelem is on the shifted C scale, so 2 is Li */
1318  for( ion=0; ion <= nelem; ion++ )
1319  {
1320  /* number of bound electrons, = atomic number for neutral */
1321  nelec = nelem - ion + 1;
1322  /* one of dima's routines to determine the number of electrons
1323  * for this species, nelem +1 to shift to physical number */
1324  /* subroutine atmdat_outer_shell(iz,in,imax,ig0,ig1)
1325  * iz - atomic number from 1 to 30 (integer)
1326  * in - number of electrons from 1 to iz (integer)
1327  * Output: imax - number of the outer shell
1328  */
1329  atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);
1330 
1331  ASSERT( imax > 0 && imax <= 10 );
1332 
1333  /* nsShells[nelem][ion] is outer shell number for ion with nelec electrons
1334  * on physics scale, with K shell being 1 */
1335  Heavy.nsShells[nelem][ion] = imax;
1336  }
1337  }
1338 
1339  /*************************************************************
1340  * *
1341  * get the Auger electron yield data set *
1342  * *
1343  *************************************************************/
1344 
1345  t_yield::Inst().init_yield();
1346 
1347  /****************************************************************
1348  * *
1349  * get the Hummer and Storey model case A, B results, these are *
1350  * the two data files e1b.dat and e2b.dat, for H and He *
1351  * *
1352  ****************************************************************/
1353 
1354  /* now get Hummer and Storey case b data, loop over H, He */
1355  /* >>chng 01 aug 08, generalized this to both case A and B, and all elements
1356  * up to HS_NZ, now 8 */
1357  for( ipZ=0; ipZ<HS_NZ; ++ipZ )
1358  {
1359  /* don't do the minor elements, Li-B */
1360  if( ipZ>1 && ipZ<5 ) continue;
1361 
1362  for( iCase=0; iCase<2; ++iCase )
1363  {
1364  /* open file with case B data */
1365  /* >>chng 01 aug 08, add HS_ to start of file names to indicate Hummer Storey origin */
1366  /* create file name for this charge
1367  * first character after e is charge number for this element,
1368  * then follows whether this is case A or case B data */
1369  sprintf( chFilename, "HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ? 'a' : 'b' );
1370 
1371  if( trace.lgTrace )
1372  fprintf( ioQQQ," atmdat_readin reading Hummer Storey emission file %s\n",chFilename );
1373 
1374  /* open the file */
1375  ioDATA = open_data( chFilename, "r" );
1376 
1377  /* read in the number of temperatures and densities*/
1378  i = fscanf( ioDATA, "%li %li ",
1379  &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );
1380 
1381  /* check that ntemp and nDensity are below NHSDIM,
1382  * set to 15 in atmdat_HS_caseB.h */
1383  assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
1384  assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM);
1385 
1386  /* loop reading in line emissivities for all temperatures*/
1387  for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
1388  {
1389  for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
1390  {
1391  long int junk, junk2 , ne;
1392  i = fscanf( ioDATA, " %lf %li %lf %c %li %ld ",
1393  &atmdat.Density[iCase][ipZ][ipDens], &junk ,
1394  &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 ,
1395  &atmdat.ncut[iCase][ipZ] );
1396 
1397  ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
1398  ASSERT( ne<=NLINEHS );
1399  for( j=0; j < ne; j++ )
1400  {
1401  i = fscanf( ioDATA, "%lf ",
1402  &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
1403  }
1404  }
1405  }
1406 
1407  /*this is end of read-in loop */
1408  fclose(ioDATA);
1409  if( trace.lgTrace )
1410  fprintf( ioQQQ," reading %s OK\n", chFilename );
1411 
1412 # if 0
1413  /* print list of densities and temperatures */
1414  for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
1415  {
1416  fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]);
1417  }
1418  fprintf(ioQQQ,"\n");
1419  for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
1420  {
1421  fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]);
1422  }
1423  fprintf(ioQQQ,"\n");
1424 # endif
1425  }
1426  }
1427 
1428  // read cross sections for neutral helium
1429  read_SH98_He1_cross_sections();
1430  // read cross sections for some he-like ions
1431  read_Helike_cross_sections();
1432 
1433  /* Verner's model FeII atom
1434  * do not call if Netzer model used, or it Fe(2) is zero
1435  * exception is when code is searching for first soln */
1436  /* read the atomic data from files */
1437  /* convert line arrays to internal form needed by code */
1438  FeIICreate();
1439 
1440  // set up spline coefficients for two-photon continuua
1441  atmdat_2phot_setSplineCoefs();
1442 
1443  return;
1444 }
1445 
1446 STATIC void read_SH98_He1_cross_sections(void)
1447 {
1448  DEBUG_ENTRY( "read_SH98_He1_cross_sections()" );
1449 
1450  FILE *ioDATA;
1451  bool lgEOL;
1452 
1453  char chPath[FILENAME_PATH_LENGTH_2],
1454  chDirectory[FILENAME_PATH_LENGTH_2],
1455  chLine[FILENAME_PATH_LENGTH_2];
1456 
1457  const int ipNUM_FILES = 10;
1458 
1459  char chFileNames[ipNUM_FILES][10] = {
1460  "p0202.3se",
1461  "p0202.3po",
1462  "p0202.3ge",
1463  "p0202.3fo",
1464  "p0202.3de",
1465  "p0202.1se",
1466  "p0202.1po",
1467  "p0202.1ge",
1468  "p0202.1fo",
1469  "p0202.1de" };
1470 
1471  HS_He1_Xsectn = ((double****)MALLOC( (size_t)(26)*sizeof(double***)));
1472  HS_He1_Energy = ((double****)MALLOC( (size_t)(26)*sizeof(double***)));
1473 
1474  // point these to NULL and use real quantum numbers rather than starting at 0
1475  HS_He1_Xsectn[0] = NULL;
1476  HS_He1_Energy[0] = NULL;
1477 
1478  for( long in = 1; in<=25; in++ )
1479  {
1480  // malloc n values of angular momentum, but not more than 5
1481  HS_He1_Xsectn[in] = ((double***)MALLOC( (size_t)(MIN2(5,in))*sizeof(double**)));
1482  HS_He1_Energy[in] = ((double***)MALLOC( (size_t)(MIN2(5,in))*sizeof(double**)));
1483  for( long il = 0; il<MIN2(5,in); il++ )
1484  {
1485  // malloc two values of spin
1486  HS_He1_Xsectn[in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
1487  HS_He1_Energy[in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
1488  HS_He1_Xsectn[in][il][0] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
1489  HS_He1_Energy[in][il][0] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
1490  HS_He1_Xsectn[in][il][1] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
1491  HS_He1_Energy[in][il][1] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
1492  }
1493  }
1494 
1495 # ifdef _WIN32
1496  strcpy( chDirectory, "sh98_he1\\pi\\" );
1497 # else
1498  strcpy( chDirectory, "sh98_he1/pi/" );
1499 # endif
1500 
1501  //HS_He1_data[25][<=4][2]
1502 
1503  for( long ipFile=0; ipFile<ipNUM_FILES; ipFile++ )
1504  {
1505  long S, L, index, N=0;
1506  long UNUSED P;
1507 
1508  strcpy( chPath, chDirectory );
1509  strcat( chPath, chFileNames[ipFile] );
1510  ioDATA = open_data( chPath, "r" );
1511 
1512  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1513  {
1514  long i=1, s;
1515  long i1, i2, i3;
1516  long numDataPoints;
1517 
1518  // first line (read above in while) is not needed except that "0 0 0" marks EOF
1519  i1 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1520  i2 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1521  i3 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1522  if( i1==0 && i2==0 && i3==0 )
1523  break;
1524 
1525  // don't need next two lines in each set
1526  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1527  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1528 
1529  i=1;
1530  // 4th line in each set identifies the quantum level
1531  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1532  S = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1533  L = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1534  P = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1535  index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1536 
1537  //indices start with unity
1538  ASSERT( index >= 1 );
1539 
1540  // index is energy order in that series and is therefore related
1541  // to principal quantum number. For triplet S must add one because
1542  // series starts at n=2.
1543  if( L==0 && S==3 )
1544  N = L + index + 1;
1545  else
1546  N = L + index;
1547 
1548  // data go up to n=25
1549  ASSERT( N<=25 );
1550 
1551  if( S==1 )
1552  s=0;
1553  else if( S==3 )
1554  s=1;
1555  else
1556  TotalInsanity();
1557 
1558  i=1;
1559  // 5th line in each set contains the number of energies
1560  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1561  //first value is not needed
1562  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1563  numDataPoints = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1564  // each set has exactly 811 data points, might as well assert this
1565  ASSERT( numDataPoints == NUM_HS98_DATA_POINTS );
1566 
1567  // don't need 6th line either
1568  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1569 
1570  // now begin reading in lines
1571  // there must be exactly numDataPoints ordered pairs,
1572  // throw an assert for any deviation
1573  for( long k=0; k<NUM_HS98_DATA_POINTS; k++ )
1574  {
1575  i=1;
1576  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1577  HS_He1_Energy[N][L][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1578  HS_He1_Xsectn[N][L][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1579  }
1580  }
1581 
1582  // we reached end of file, assert last quantum number was 25
1583  ASSERT( N==25 );
1584 
1585  fclose( ioDATA );
1586  }
1587 
1588  return;
1589 }
1590 
1591 STATIC void read_Helike_cross_sections(void)
1592 {
1593  DEBUG_ENTRY( "read_Helike_cross_sections()" );
1594 
1595  FILE *ioDATA;
1596  bool lgEOL;
1597 
1598  char chLine[FILENAME_PATH_LENGTH_2];
1599  char chFileName[23] = "helike_pcs_topbase.dat";
1600 
1601  // the data only go as high as n=10
1602  const int MaxN = 10;
1603  long last_i1=0;
1604 
1605  // data will be used up to calcium (nelem=19)
1606  // data exists up to iron but we will not use it because it has strong resonance
1607  // features, but is nonetheless very nearly hydrogenic
1608  OP_Helike_Xsectn = ((double*****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(double****)));
1609  OP_Helike_Energy = ((double*****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(double****)));
1610  OP_Helike_NumPts = ((long****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(long***)));
1611 
1612  // point these to NULL because we will never use them
1613  OP_Helike_Xsectn[ipHYDROGEN] = NULL;
1614  OP_Helike_Energy[ipHYDROGEN] = NULL;
1615  OP_Helike_NumPts[ipHYDROGEN] = 0;
1616  OP_Helike_Xsectn[ipHELIUM] = NULL;
1617  OP_Helike_Energy[ipHELIUM] = NULL;
1618  OP_Helike_NumPts[ipHELIUM] = 0;
1619 
1620  for( long nelem = ipLITHIUM; nelem<= ipCALCIUM; nelem++ )
1621  {
1622  // malloc principal quantum number
1623  OP_Helike_Xsectn[nelem] = ((double****)MALLOC( (size_t)(MaxN+1)*sizeof(double***)));
1624  OP_Helike_Energy[nelem] = ((double****)MALLOC( (size_t)(MaxN+1)*sizeof(double***)));
1625  OP_Helike_NumPts[nelem] = ((long***)MALLOC( (size_t)(MaxN+1)*sizeof(long**)));
1626 
1627  for( long in = 1; in<=MaxN; in++ )
1628  {
1629  // malloc angular momentum
1630  OP_Helike_Xsectn[nelem][in] = ((double***)MALLOC( (size_t)(in)*sizeof(double**)));
1631  OP_Helike_Energy[nelem][in] = ((double***)MALLOC( (size_t)(in)*sizeof(double**)));
1632  OP_Helike_NumPts[nelem][in] = ((long**)MALLOC( (size_t)(in)*sizeof(long*)));
1633  for( long il = 0; il<in; il++ )
1634  {
1635  // malloc two values of spin
1636  OP_Helike_Xsectn[nelem][in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
1637  OP_Helike_Energy[nelem][in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
1638  OP_Helike_NumPts[nelem][in][il] = ((long*)MALLOC( (size_t)(2)*sizeof(long)));
1639  // point these to NULL for now, won't know how many we need until we begin parsing
1640  OP_Helike_Xsectn[nelem][in][il][0] = NULL;
1641  OP_Helike_Energy[nelem][in][il][0] = NULL;
1642  OP_Helike_NumPts[nelem][in][il][0] = 0;
1643  OP_Helike_Xsectn[nelem][in][il][1] = NULL;
1644  OP_Helike_Energy[nelem][in][il][1] = NULL;
1645  OP_Helike_NumPts[nelem][in][il][1] = 0;
1646  }
1647  }
1648  }
1649 
1650  ioDATA = open_data( chFileName, "r" );
1651 
1652  // Header looks like this:
1653  // ================================================
1654  // I NZ NE ISLP ILV E(RYD) NP
1655  // ================================================
1656  // 1 3 2 100 1 -5.53159E+00 55
1657 
1658  // so skip the first three lines.
1659  for( long i=0; i<3; i++)
1660  {
1661  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1662  {
1663  fprintf( ioQQQ,"PROBLEM corruption in TOPbase Helike pcs datafile.\nSorry\n" );
1664  cdEXIT(EXIT_FAILURE);
1665  }
1666  }
1667 
1668  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1669  {
1670  long i=1;
1671  long n, l, s;
1672  long i1, i2, i3, i4, i5, i7;
1673  double i6;
1674 
1675  i1 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1676  i2 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1677  i3 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1678  i4 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1679  i5 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1680  i6 = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1681  i7 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1682 
1683  if( lgEOL )
1684  {
1685  fprintf( ioQQQ,"PROBLEM corruption in TOPbase Helike pcs datafile.\nSorry\n" );
1686  cdEXIT(EXIT_FAILURE);
1687  }
1688 
1689  // this marks end of data.
1690  if( i1==i2 && i1==i3 && i1==i4 && i1==i5 && i1==i7 && i1==-1 )
1691  break;
1692 
1693  // first parameter is level index (overall, not just for series or charge)
1694  // check that it is as expected
1695  ASSERT( i1>0 && i1==(last_i1 + 1) && i1<=795 );
1696  last_i1 = i1;
1697  // second parameter is nuclear charge
1698  ASSERT( (i2-1)<=ipCALCIUM && (i2-1)>=ipLITHIUM );
1699  // third parameter must be 2 for helike.
1700  ASSERT( i3==2 );
1701  // fourth parameter is (2S+1)*100+L*10+P
1702  ASSERT( i4>=100 && i4<400 );
1703  if( i4 >= 300 )
1704  s=1;
1705  else
1706  s=0;
1707  l = (i4 - (2*s+1)*100)/10;
1708  // data only goes up to l=2
1709  ASSERT( l<=2 );
1710  // fifth is index in the series, related to principal quantum number
1711  ASSERT( i5>=1 && i5<=10 );
1712  if( s==1 && l==0 )
1713  n = i5 + 1;
1714  else
1715  n = i5 + l;
1716  ASSERT( l<=MaxN );
1717  // sixth is threshhold energy, don't need but assert negative
1718  // \todo 3 save this and renorm cross-section with ratio of actual to recorded Eth?
1719  ASSERT( i6 < 0. );
1720  // seventh parameter is number of data points, can be zero, use to MALLOC otherwise
1721  OP_Helike_NumPts[i2-1][n][l][s] = i7;
1722  if( i7==0 )
1723  continue;
1724 
1725  ASSERT( i7 > 0 );
1726  ASSERT( OP_Helike_Xsectn[i2-1][n][l][s]==NULL );
1727  ASSERT( OP_Helike_Energy[i2-1][n][l][s]==NULL );
1728  OP_Helike_Xsectn[i2-1][n][l][s] = ((double*)MALLOC( (size_t)(i7)*sizeof(double)));
1729  OP_Helike_Energy[i2-1][n][l][s] = ((double*)MALLOC( (size_t)(i7)*sizeof(double)));
1730 
1731  // now begin reading in lines
1732  // there must be exactly i7 ordered pairs,
1733  for( long k=0; k<i7; k++ )
1734  {
1735  i=1;
1736  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
1737  OP_Helike_Energy[i2-1][n][l][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1738  OP_Helike_Xsectn[i2-1][n][l][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1739 
1740  // make sure data is well-behaved
1741  if( k > 0 )
1742  {
1743  ASSERT( OP_Helike_Energy[i2-1][n][l][s][k] > OP_Helike_Energy[i2-1][n][l][s][k-1] );
1744  ASSERT( OP_Helike_Xsectn[i2-1][n][l][s][k] < OP_Helike_Xsectn[i2-1][n][l][s][k-1] );
1745  }
1746 
1747  // try to read one more item off the line and verify that lgEOL is true
1748  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
1749  ASSERT( lgEOL );
1750  }
1751  }
1752 
1753  fclose( ioDATA );
1754 
1755  return;
1756 }
1757 
1758 t_yield::t_yield()
1759 {
1760  /* preset two arrays that will hold auger data
1761  * set to very large values so that code will blow
1762  * if not set properly below */
1763  for( int nelem=0; nelem < LIMELM; nelem++ )
1764  {
1765  for( int ion=0; ion < LIMELM; ion++ )
1766  {
1767  for( int ns=0; ns < 7; ns++ )
1768  {
1769  n_elec_eject[nelem][ion][ns] = LONG_MAX;
1770  for( int nelec=0; nelec < 10; nelec++ )
1771  {
1772  frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
1773  }
1774  }
1775  }
1776  }
1777 
1778  lgKillAuger = false;
1779 }
1780 
1781 void t_yield::init_yield()
1782 {
1783  char chLine[FILENAME_PATH_LENGTH_2];
1784  const char* chFilename;
1785 
1786  /* following is double for sscanf to work */
1787  double temp[15];
1788 
1789  DEBUG_ENTRY( "init_yield()" );
1790 
1791  /*************************************************************
1792  * *
1793  * get the Auger electron yield data set *
1794  * *
1795  *************************************************************/
1796 
1797  /* NB NB -- This test of Heavy.nsShells remains here to assure
1798  * that it contains meaningful values since needed below, once
1799  * t_Heavy is a Singleton, it can be removed !!! */
1800  ASSERT( Heavy.nsShells[2][0] > 0 );
1801 
1802  /* hydrogen and helium will not be done below, so set yields here*/
1803  n_elec_eject[ipHYDROGEN][0][0] = 1;
1804  n_elec_eject[ipHELIUM][0][0] = 1;
1805  n_elec_eject[ipHELIUM][1][0] = 1;
1806 
1807  frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
1808  frac_elec_eject[ipHELIUM][0][0][0] = 1;
1809  frac_elec_eject[ipHELIUM][1][0][0] = 1;
1810 
1811  /* open file auger.dat, yield data file that came from
1812  * >>refer all auger Kaastra, J. S., & Mewe, R. 1993, A&AS, 97, 443-482 */
1813  chFilename = "mewe_nelectron.dat";
1814 
1815  if( trace.lgTrace )
1816  fprintf( ioQQQ, " init_yield reading %s\n", chFilename );
1817 
1818  FILE *ioDATA;
1819 
1820  /* open the file */
1821  ioDATA = open_data( chFilename, "r" );
1822 
1823  /* now read in all auger yields
1824  * will do elements from li on up,
1825  * skip any line starting with #
1826  * this loop goes from lithium to Zn */
1827  for( int nelem=2; nelem < LIMELM; nelem++ )
1828  {
1829  /* nelem is on the shifted C scale, so 2 is Li */
1830  for( int ion=0; ion <= nelem; ion++ )
1831  {
1832  for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1833  {
1834  char ch1 = '#';
1835  /* the * is the old comment char, accept it, but really want # */
1836  while( ch1 == '#' || ch1 == '*' )
1837  {
1838  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1839  {
1840  fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns );
1841  cdEXIT(EXIT_FAILURE);
1842  }
1843  ch1 = chLine[0];
1844  }
1845  /*printf("string is %s\n",chLine );*/
1846  sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
1847  &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1848  &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
1849  &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] );
1850  n_elec_eject[nelem][ion][ns] = (long int)temp[3];
1851 
1852  ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 );
1853 
1854  /* can pop off up to 10 auger electrons, these are the probabilities*/
1855  for( int j=0; j < 10; j++ )
1856  {
1857  frac_elec_eject[nelem][ion][ns][j] = (realnum)temp[j+4];
1858  ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. );
1859  }
1860  }
1861  }
1862  /* activate this print statement to get yield for k-shell of all atoms */
1863  /*fprintf(ioQQQ,"yyyield\t%li\t%.4e\n", nelem+1, frac_elec_eject[nelem][0][0][0] );*/
1864  }
1865 
1866  fclose( ioDATA );
1867 
1868  if( trace.lgTrace )
1869  {
1870  /* this is set with no auger command */
1871  if( lgKillAuger )
1872  fprintf( ioQQQ, " Auger yields will be killed.\n");
1873  fprintf( ioQQQ, " reading %s OK\n", chFilename );
1874  }
1875 
1876  /* open file mewe_fluor.dat, yield data file that came from
1877  * >>refer all auger Kaastra, J. S., & Mewe, R. 1993, 97, 443-482 */
1878  chFilename = "mewe_fluor.dat";
1879 
1880  if( trace.lgTrace )
1881  fprintf( ioQQQ, " init_yield reading %s\n", chFilename );
1882 
1883  /* open the file */
1884  ioDATA = open_data( chFilename, "r" );
1885 
1886  /* now read in all auger yields
1887  * will do elements from li on up,
1888  * skip any line starting with # */
1889  do
1890  {
1891  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1892  {
1893  fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 );
1894  cdEXIT(EXIT_FAILURE);
1895  }
1896  }
1897  while( chLine[0] == '#' );
1898 
1899  bool lgEOL = false;
1900 
1901  nfl_lines = 0;
1902  do
1903  {
1904  const int NKM = 10;
1905  int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
1906  int nAuger;
1907 
1908  if( nfl_lines >= MEWE_FLUOR )
1909  TotalInsanity();
1910 
1911  /*printf("string is %s\n",chLine );*/
1912  sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf",
1913  &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1914  &temp[5], &temp[6] );
1915 
1916  /* the atomic number, C is 5 */
1917  nfl_nelem[nfl_lines] = (int)temp[0]-1;
1918  ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM );
1919 
1920  /* the ion stage for target, atom is 0 */
1921  nfl_ion[nfl_lines] = (int)temp[1]-1;
1922  ASSERT( nfl_ion[nfl_lines] >= 0 && nfl_ion[nfl_lines] <= nfl_nelem[nfl_lines]+1 );
1923 
1924  /* the target's shell */
1925  nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1];
1926  ASSERT( nfl_nshell[nfl_lines] >= 0 &&
1927  /* nsShells is shell number, where K is 1 */
1928  nfl_nshell[nfl_lines] < Heavy.nsShells[nfl_nelem[nfl_lines]][nfl_ion[nfl_lines]]-1 );
1929 
1930  /* this is the number of Auger electrons ejected */
1931  nAuger = (int)temp[3];
1932  /* so this is the spectrum of the photons that are emitted */
1933  nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1;
1934  /* must be gt 0 since at least photoelectron comes off */
1935  ASSERT( nfl_ion_emit[nfl_lines] > 0 && nfl_ion_emit[nfl_lines] <= nfl_nelem[nfl_lines]+1);
1936 
1937  /* this is the type of line as defined in their paper */
1938  nfl_nLine[nfl_lines] = (int)temp[4];
1939 
1940  /* energy in Ryd */
1941  fl_energy[nfl_lines] = (realnum)temp[5] / (realnum)EVRYD;
1942  ASSERT( fl_energy[nfl_lines] > 0. );
1943 
1944  /* fluor yield */
1945  fl_yield[nfl_lines] = (realnum)temp[6];
1946  /* NB cannot assert <=1 since data file has yields around 1.3 - 1.4 */
1947  ASSERT( fl_yield[nfl_lines] >= 0 );
1948 
1949  ++nfl_lines;
1950 
1951  do
1952  {
1953  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1954  lgEOL = true;
1955  }
1956  while( chLine[0]=='#' && !lgEOL );
1957  }
1958  while( !lgEOL );
1959 
1960  fclose( ioDATA );
1961  if( trace.lgTrace )
1962  fprintf( ioQQQ, " reading %s OK\n", chFilename );
1963 }
1964 
1965 // read autoionization data from Badnell data file
1966 STATIC void ReadBadnellAIData(const string& fnam, // filename containing the Badnell data
1967  long nelem, // nelem is on C scale
1968  long ion, // ion is on C scale
1969  TransitionList& UTA, // UTA lines will be pushed on this stack
1970  bitset<IS_TOP> Skip) // option to skip transitions from a particular shell
1971 {
1972  DEBUG_ENTRY( "ReadBadnellAIData()" );
1973 
1974  if( trace.lgTrace )
1975  fprintf( ioQQQ," ReadBadnellAIData reading %s\n", fnam.c_str() );
1976 
1977  fstream ioDATA;
1978  open_data( ioDATA, fnam.c_str(), mode_r );
1979 
1980  string line;
1981  getline( ioDATA, line );
1982  ASSERT( line.substr(0,4) == "SEQ=" );
1983  getline( ioDATA, line );
1984  getline( ioDATA, line );
1985  // we don't need the parent level data, so we will skip it...
1986  ASSERT( line.substr(3,21) == "PARENT LEVEL INDEXING" );
1987  int nParent;
1988  istringstream iss( line.substr(65,4) );
1989  iss >> nParent;
1990  // data lines containing data for all levels of the parent ion will span nMulti lines
1991  int nMulti = (nParent+5)/6;
1992  for( int i=0; i < nParent+5; ++i )
1993  getline( ioDATA, line );
1994 
1995  // here starts the header for the level data we need
1996  ASSERT( line.substr(3,26) == "IC RESOLVED LEVEL INDEXING" );
1997  int nLevel;
1998  istringstream iss2( line.substr(63,6) );
1999  iss2 >> nLevel;
2000  // skip rest of header
2001  for( int i=0; i < 3; ++i )
2002  getline( ioDATA, line );
2003 
2004  // now get the level data
2005  vector<t_BadnellLevel> level( nLevel );
2006  for( int i=0; i < nLevel; ++i )
2007  {
2008  getline( ioDATA, line );
2009  istringstream iss3( line );
2010  int indx, irsl;
2011  iss3 >> indx >> irsl;
2012  level[indx-1].irsl = irsl;
2013  level[indx-1].config = line.substr(16,20);
2014  istringstream iss4( line.substr(37,1) );
2015  iss4 >> level[indx-1].S;
2016  istringstream iss5( line.substr(39,1) );
2017  iss5 >> level[indx-1].L;
2018  istringstream iss6( line.substr(41,4) );
2019  double J;
2020  iss6 >> J;
2021  level[indx-1].g = nint(2.*J + 1.);
2022  istringstream iss7( line.substr(46,11) );
2023  iss7 >> level[indx-1].energy;
2024 
2025  // which inner shell has been excited?
2026  level[indx-1].lgAutoIonizing = ( line[57] == '*' );
2027  if( level[indx-1].lgAutoIonizing )
2028  {
2029  if( level[indx-1].config.find( "1S1" ) != string::npos )
2030  level[indx-1].WhichShell = IS_K_SHELL;
2031  else if( level[indx-1].config.find( "2S1" ) != string::npos )
2032  level[indx-1].WhichShell = IS_L1_SHELL;
2033  else if( level[indx-1].config.find( "2P5" ) != string::npos )
2034  level[indx-1].WhichShell = IS_L2_SHELL;
2035  else
2036  TotalInsanity();
2037  }
2038  else
2039  {
2040  level[indx-1].WhichShell = IS_NONE;
2041  }
2042  }
2043 
2044  // levels are done, now move on to the lines
2045  // first search for start of the header
2046  while( getline( ioDATA, line ) )
2047  {
2048  if( line.find( "IRSL IRSL" ) != string::npos )
2049  break;
2050  }
2051  // skip rest of the header
2052  for( int i=0; i < nMulti-1; ++i )
2053  getline( ioDATA, line );
2054 
2055  // blank line that will be pushed on the UTA line stack
2056  qList BlankStates(1);
2057  TransitionList BlankList("BlankList",&BlankStates,1);
2058  TransitionList::iterator BlankLine = BlankList.begin();
2059  (*BlankLine).Junk();
2060 
2061  // start reading the line data
2062  while( getline( ioDATA, line ) )
2063  {
2064  // have we reached the end of the line data?
2065  if( line.size() < 10 )
2066  break;
2067 
2068  // test if there is an autoionization rate on this line
2069  // if not, it only contains radiative data and we skip it...
2070  if( line.size() < 50 )
2071  continue;
2072 
2073  // there may be an asterisk here; wipe it out since we don't need it
2074  line[19] = ' ';
2075 
2076  int irsl_lo, irsl_hi, dum;
2077  double edif, Bij, Rji, Aai;
2078  istringstream iss8( line );
2079  // irsl_lo: index for lower level of transition
2080  // irsl_hi: index for upper level of transition
2081  // edif: energy difference between levels in Ryd
2082  // Bij: UPWARD Einstein A for transition irsl_lo -> irsl_hi
2083  // Rji: sum of Aji for all radiative transitions to lower levels
2084  // Aai: autoionization rate from level irsl_hi
2085  iss8 >> irsl_lo >> irsl_hi >> dum >> dum >> edif >> Bij >> Rji >> Aai;
2086  // ind_lo and ind_hi are on C scale
2087  int ind_lo = irsl2ind( level, irsl_lo );
2088  int ind_hi = irsl2ind( level, irsl_hi );
2089  ASSERT( level[ind_hi].lgAutoIonizing );
2090  // skip rest of partial autoionization rates
2091  for( int i=0; i < nMulti-1; ++i )
2092  getline( ioDATA, line );
2093  // skip this transition if it does not originate from ground level
2094  // or if the user requested to skip excitations from this inner shell
2095  if( ind_lo == 0 && !Skip[level[ind_hi].WhichShell] )
2096  {
2097  UTA.push_back( *BlankLine );
2098  InitTransition( UTA.back() );
2099 
2100  // t_emission has nelem and ion on fortran scale...
2101  (*UTA.back().Hi()).nelem() = nelem+1;
2102  (*UTA.back().Hi()).IonStg() = ion+1;
2103 
2104  (*UTA.back().Hi()).g() = (realnum)level[ind_hi].g;
2105  (*UTA.back().Lo()).g() = (realnum)level[ind_lo].g;
2106 
2107  double WavNum = edif*RYD_INF;
2108 
2109  /* wavelength in Angstroms */
2110  UTA.back().WLAng() = (realnum)(1e8/WavNum);
2111  UTA.back().EnergyWN() = (realnum)WavNum;
2112 
2113  /* store branching ratio for autoionization */
2114  double frac_ioniz = Aai/(Rji + Aai);
2115  ASSERT( frac_ioniz >= 0. && frac_ioniz <= 1. );
2116  UTA.back().Emis().AutoIonizFrac() = (realnum)frac_ioniz;
2117 
2118  /* this is true spontaneous rate for doubly excited state to ground
2119  * and is used for pumping, and also relaxing back to inner shell */
2120  /* Badnell gives UPWARD transition rate Alu, an unusual notation,
2121  * convert it here to the normal downward transition rate Aul */
2122  UTA.back().Emis().Aul() = (realnum)(Bij*(*UTA.back().Lo()).g()/(*UTA.back().Hi()).g());
2123  UTA.back().Emis().gf() =
2124  (realnum)GetGF( UTA.back().Emis().Aul(), UTA.back().EnergyWN(), (*UTA.back().Hi()).g() );
2125 
2126  UTA.back().Emis().iRedisFun() = ipPRD;
2127 
2128  UTA.back().Emis().dampXvel() = (realnum)((Rji+Aai)/UTA.back().EnergyWN()/PI4);
2129 
2130  ASSERT( UTA.back().Emis().dampXvel() > 0. );
2131 
2132  // remove this line if it is too weak
2133  if( UTA.back().Emis().gf() < level[ind_lo].g * f_cutoff )
2134  UTA.pop_back();
2135  }
2136  }
2137 
2138  // perform a sanity check on the tail of the file
2139  getline( ioDATA, line );
2140  ASSERT( line.substr(3,7) == "NRSLMX=" );
2141 
2142  ioDATA.close();
2143 
2144  if( trace.lgTrace )
2145  fprintf( ioQQQ, " reading %s OK\n", fnam.c_str() );
2146 }
2147 
2148 inline void InitTransition(const TransitionProxy& t)
2149 {
2150  t.AddHiState();
2151  t.AddLoState();
2152  t.AddLine2Stack();
2153 }
2154 
2155 inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl)
2156 {
2157  for( unsigned int i=0; i < level.size(); ++i )
2158  {
2159  if( level[i].irsl == irsl )
2160  return (int)i;
2161  }
2162  // we should never get here...
2163  TotalInsanity();
2164 }
ipT4368
long ipT4368
Definition: atmdat_readin.cpp:50
ipTOI46
long ipTOI46
Definition: atmdat_readin.cpp:51
ipT3934
long ipT3934
Definition: atmdat_readin.cpp:81
ipT834
long ipT834
Definition: atmdat_readin.cpp:48
ipT990
long ipT990
Definition: atmdat_readin.cpp:43
t_struc::coolstr
double * coolstr
Definition: struc.h:78
ReadBadnellAIData
STATIC void ReadBadnellAIData(const string &fnam, long nelem, long ion, TransitionList &UTA, bitset< IS_TOP > Skip)
Definition: atmdat_readin.cpp:1966
ipS4_1405
long ipS4_1405
Definition: atmdat_readin.cpp:90
t_atmdat::lgLamdaOn
bool lgLamdaOn
Definition: atmdat.h:239
ipT8542
long ipT8542
Definition: atmdat_readin.cpp:82
ipfsNa490
long ipfsNa490
Definition: atmdat_readin.cpp:57
t_struc::gas_phase
realnum ** gas_phase
Definition: struc.h:75
atmdat_2phot_setSplineCoefs
void atmdat_2phot_setSplineCoefs()
Definition: atmdat_2photon.cpp:213
read_SH98_He1_cross_sections
STATIC void read_SH98_He1_cross_sections(void)
Definition: atmdat_readin.cpp:1446
ipTS19
long ipTS19
Definition: atmdat_readin.cpp:77
ipT770
long ipT770
Definition: atmdat_readin.cpp:54
ipP713
long ipP713
Definition: atmdat_readin.cpp:75
t_yield::nfl_nLine
int nfl_nLine[MEWE_FLUOR]
Definition: yield.h:34
ipT63
long ipT63
Definition: atmdat_readin.cpp:53
ipTFe16
long ipTFe16
Definition: atmdat_readin.cpp:87
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
ipT1037
long ipT1037
Definition: atmdat_readin.cpp:54
t_BadnellLevel::WhichShell
exc_type WhichShell
Definition: atmdat_readin.cpp:126
ipAl529
long ipAl529
Definition: atmdat_readin.cpp:66
lines.h
yield.h
ipNi1_11m
long ipNi1_11m
Definition: atmdat_readin.cpp:97
ipTNe24
long ipTNe24
Definition: atmdat_readin.cpp:56
ipAl09204
long ipAl09204
Definition: atmdat_readin.cpp:67
ipT1808
long ipT1808
Definition: atmdat_readin.cpp:69
cs1_flag_lev2
realnum * cs1_flag_lev2
Definition: taulines.cpp:40
t_struc::xIonDense
realnum *** xIonDense
Definition: struc.h:64
t_yield::nfl_nshell
int nfl_nshell[MEWE_FLUOR]
Definition: yield.h:32
ipSi2_2350
long ipSi2_2350
Definition: atmdat_readin.cpp:94
ipTFe34
long ipTFe34
Definition: atmdat_readin.cpp:87
ipTr48
long ipTr48
Definition: atmdat_readin.cpp:87
EmissionProxy::AutoIonizFrac
realnum & AutoIonizFrac() const
Definition: emission.h:583
ipS4_1407
long ipS4_1407
Definition: atmdat_readin.cpp:90
ipTuv3
long ipTuv3
Definition: atmdat_readin.cpp:86
struc.h
dense
t_dense dense
Definition: dense.cpp:24
NHSDIM
#define NHSDIM
Definition: atmdat.h:123
elementnames.h
ipTS11
long ipTS11
Definition: atmdat_readin.cpp:77
ipTCa3
long ipTCa3
Definition: atmdat_readin.cpp:83
ipxMg52417
long ipxMg52417
Definition: atmdat_readin.cpp:60
ipSc13264
long ipSc13264
Definition: atmdat_readin.cpp:84
t_atmdat::ntemp
long int ntemp[2][HS_NZ]
Definition: atmdat.h:188
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:175
ipTSi25
long ipTSi25
Definition: atmdat_readin.cpp:73
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
ipO4_1400
long ipO4_1400
Definition: atmdat_readin.cpp:91
atmdat_readin
void atmdat_readin(void)
Definition: atmdat_readin.cpp:149
database_readin
void database_readin(void)
Definition: species.cpp:42
geometry.h
nLevel1
long int nLevel1
Definition: taulines.cpp:28
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
ipT1304
long ipT1304
Definition: atmdat_readin.cpp:49
ipAlI3090
long ipAlI3090
Definition: atmdat_readin.cpp:64
ipCl973
long ipCl973
Definition: atmdat_readin.cpp:78
ipNA_LIKE
const int ipNA_LIKE
Definition: iso.h:72
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
ipTFe46
long ipTFe46
Definition: atmdat_readin.cpp:88
ipSii2215
long ipSii2215
Definition: atmdat_readin.cpp:69
UNUSED
#define UNUSED
Definition: cpu.h:14
ipT671
long ipT671
Definition: atmdat_readin.cpp:45
realnum
float realnum
Definition: cddefines.h:103
iterations
t_iterations iterations
Definition: iterations.cpp:5
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
TransitionProxy::AddHiState
void AddHiState() const
Definition: transition.cpp:687
lines_setup
void lines_setup(void)
Definition: atmdat_lines_setup.cpp:19
N
static const int N
Definition: thirdparty.cpp:2814
nWindLine
long nWindLine
Definition: cdinit.cpp:19
ipC2_2325
long ipC2_2325
Definition: atmdat_readin.cpp:93
ipO4_1397
long ipO4_1397
Definition: atmdat_readin.cpp:91
STATIC
#define STATIC
Definition: cddefines.h:97
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:307
nUTA
long int nUTA
Definition: taulines.cpp:26
TransitionList::size
size_t size(void) const
Definition: transition.h:297
ipT1661
long ipT1661
Definition: atmdat_readin.cpp:48
ipSii2518
long ipSii2518
Definition: atmdat_readin.cpp:68
ipTSi41
long ipTSi41
Definition: atmdat_readin.cpp:72
ipTSi499
long ipTSi499
Definition: atmdat_readin.cpp:72
mole.h
ipNI_pumpIndirect
long ipNI_pumpIndirect
Definition: atmdat_readin.cpp:99
NI_NDP
const int NI_NDP
Definition: taulines.h:90
TransitionList::resize
void resize(size_t newsize)
Definition: transition.h:285
t_struc::testr
realnum * testr
Definition: struc.h:25
ipT374g
long ipT374g
Definition: atmdat_readin.cpp:44
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:403
ipT270
long ipT270
Definition: atmdat_readin.cpp:41
ipT2804
long ipT2804
Definition: atmdat_readin.cpp:59
ipfsMg755
long ipfsMg755
Definition: atmdat_readin.cpp:64
t_yield::nfl_nelem
int nfl_nelem[MEWE_FLUOR]
Definition: yield.h:30
ipTO88
long ipTO88
Definition: atmdat_readin.cpp:53
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
ipS4_1398
long ipS4_1398
Definition: atmdat_readin.cpp:90
ipT1122
long ipT1122
Definition: atmdat_readin.cpp:88
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
ipFeI3884
long ipFeI3884
Definition: atmdat_readin.cpp:85
t_struc::drad
realnum * drad
Definition: struc.h:53
ipIRON
const int ipIRON
Definition: cddefines.h:330
t_struc::drad_last
realnum * drad_last
Definition: struc.h:58
ipSi2_2336
long ipSi2_2336
Definition: atmdat_readin.cpp:94
ipT52
long ipT52
Definition: atmdat_readin.cpp:53
ipAl6912
long ipAl6912
Definition: atmdat_readin.cpp:66
InitTransition
void InitTransition(const TransitionProxy &t)
Definition: atmdat_readin.cpp:2148
NUM_HS98_DATA_POINTS
#define NUM_HS98_DATA_POINTS
Definition: atmdat.h:127
t_BadnellLevel::irsl
int irsl
Definition: atmdat_readin.cpp:120
trace.h
t_struc::StatesElem
realnum **** StatesElem
Definition: struc.h:67
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:101
ipS1_56m
long ipS1_56m
Definition: atmdat_readin.cpp:95
ipT780
long ipT780
Definition: atmdat_readin.cpp:54
TransitionList::begin
iterator begin(void)
Definition: transition.h:305
ipT310
long ipT310
Definition: atmdat_readin.cpp:40
ipTAr7
long ipTAr7
Definition: atmdat_readin.cpp:79
ipfsNa373
long ipfsNa373
Definition: atmdat_readin.cpp:57
ProxyIterator
Definition: proxy_iterator.h:58
t_struc::drad_x_fillfac
realnum * drad_x_fillfac
Definition: struc.h:27
t_struc::hiistr
realnum * hiistr
Definition: struc.h:29
t_yield::nfl_lines
long int nfl_lines
Definition: yield.h:41
read_Helike_cross_sections
STATIC void read_Helike_cross_sections(void)
Definition: atmdat_readin.cpp:1591
ipT370
long ipT370
Definition: atmdat_readin.cpp:42
t_atmdat::ncut
long int ncut[2][HS_NZ]
Definition: atmdat.h:188
lines_service.h
ipT1198
long ipT1198
Definition: atmdat_readin.cpp:76
ipTMg6
long ipTMg6
Definition: atmdat_readin.cpp:63
ipNi1_7m
long ipNi1_7m
Definition: atmdat_readin.cpp:97
ipFeI3729
long ipFeI3729
Definition: atmdat_readin.cpp:85
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
t_BadnellLevel::t_BadnellLevel
t_BadnellLevel()
Definition: atmdat_readin.cpp:127
ipTMg14
long ipTMg14
Definition: atmdat_readin.cpp:63
dynamics.h
strstr_s
const char * strstr_s(const char *haystack, const char *needle)
Definition: cddefines.h:1429
t_yield::t_yield
t_yield()
Definition: atmdat_readin.cpp:1758
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
ipTAr22
long ipTAr22
Definition: atmdat_readin.cpp:79
ipAlI3957
long ipAlI3957
Definition: atmdat_readin.cpp:64
t_yield::nfl_ion
int nfl_ion[MEWE_FLUOR]
Definition: yield.h:31
struc
t_struc struc
Definition: struc.cpp:6
iso.h
ipT58
long ipT58
Definition: atmdat_readin.cpp:63
PI4
const UNUSED double PI4
Definition: physconst.h:35
ipC2_2327
long ipC2_2327
Definition: atmdat_readin.cpp:93
t_BadnellLevel::energy
realnum energy
Definition: atmdat_readin.cpp:124
ipAL_LIKE
const int ipAL_LIKE
Definition: iso.h:74
t_atmdat::nDensity
long int nDensity[2][HS_NZ]
Definition: atmdat.h:188
ipTAl48
long ipTAl48
Definition: atmdat_readin.cpp:68
ipP0260
long ipP0260
Definition: atmdat_readin.cpp:74
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_struc::H2_abund
realnum * H2_abund
Definition: struc.h:72
ipxMg72569
long ipxMg72569
Definition: atmdat_readin.cpp:62
t_yield::nfl_ion_emit
int nfl_ion_emit[MEWE_FLUOR]
Definition: yield.h:33
t_struc::depth
realnum * depth
Definition: struc.h:51
TransitionList::push_back
void push_back(const TransitionProxy &tr)
Definition: transition.h:313
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
atmdat.h
HS_NZ
#define HS_NZ
Definition: atmdat.h:125
t_yield::lgKillAuger
bool lgKillAuger
Definition: yield.h:44
ipTOI29
long ipTOI29
Definition: atmdat_readin.cpp:51
ipS1_25m
long ipS1_25m
Definition: atmdat_readin.cpp:95
ipCl1_11m
long ipCl1_11m
Definition: atmdat_readin.cpp:95
MIN2
#define MIN2
Definition: cddefines.h:761
TransitionProxy
Definition: transition.h:23
t_struc::nzlim
long int nzlim
Definition: struc.h:19
ipT1335
long ipT1335
Definition: atmdat_readin.cpp:38
ipT1486
long ipT1486
Definition: atmdat_readin.cpp:43
ipfsNa421
long ipfsNa421
Definition: atmdat_readin.cpp:57
ipT1656
long ipT1656
Definition: atmdat_readin.cpp:37
ipT1895
long ipT1895
Definition: atmdat_readin.cpp:70
ipT209
long ipT209
Definition: atmdat_readin.cpp:46
ipT1305
long ipT1305
Definition: atmdat_readin.cpp:71
ipTNe36
long ipTNe36
Definition: atmdat_readin.cpp:56
ipT3969
long ipT3969
Definition: atmdat_readin.cpp:82
t_BadnellLevel::g
int g
Definition: atmdat_readin.cpp:123
ipT1207
long ipT1207
Definition: atmdat_readin.cpp:70
ipT7291
long ipT7291
Definition: atmdat_readin.cpp:83
mode_r
const ios_base::openmode mode_r
Definition: cpu.h:212
ipT1039
long ipT1039
Definition: atmdat_readin.cpp:50
ipT895
long ipT895
Definition: atmdat_readin.cpp:55
MeweCoef
t_MeweCoef MeweCoef
Definition: mewecoef.cpp:5
t_atmdat::lgStoutOn
bool lgStoutOn
Definition: atmdat.h:241
ipT315
long ipT315
Definition: atmdat_readin.cpp:45
ipSi10143
long ipSi10143
Definition: atmdat_readin.cpp:72
HS_He1_Energy
double **** HS_He1_Energy
Definition: atmdat.cpp:9
t_struc::windv
realnum * windv
Definition: struc.h:35
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
ipT4561
long ipT4561
Definition: atmdat_readin.cpp:60
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
ipxK03462
long ipxK03462
Definition: atmdat_readin.cpp:80
t_atmdat::lgChiantiHybrid
bool lgChiantiHybrid
Definition: atmdat.h:233
ipMgI2853
long ipMgI2853
Definition: atmdat_readin.cpp:58
ipfsMg790
long ipfsMg790
Definition: atmdat_readin.cpp:64
ipT1243
long ipT1243
Definition: atmdat_readin.cpp:43
ipT1260
long ipT1260
Definition: atmdat_readin.cpp:71
ipT9830
long ipT9830
Definition: atmdat_readin.cpp:37
ipT5577
long ipT5577
Definition: atmdat_readin.cpp:47
mewecoef.h
dense.h
ipP0233
long ipP0233
Definition: atmdat_readin.cpp:74
ipTAl568
long ipTAl568
Definition: atmdat_readin.cpp:68
eina
double eina(double gf, double enercm, double gup)
Definition: lines_service.cpp:84
ipT1194
long ipT1194
Definition: atmdat_readin.cpp:75
OP_Helike_Energy
double ***** OP_Helike_Energy
Definition: atmdat.cpp:11
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
ipxNa6143
long ipxNa6143
Definition: atmdat_readin.cpp:58
t_atmdat::lgChiantiOn
bool lgChiantiOn
Definition: atmdat.h:231
ipT2670
long ipT2670
Definition: atmdat_readin.cpp:65
trace
t_trace trace
Definition: trace.cpp:5
ipT280
long ipT280
Definition: atmdat_readin.cpp:40
TransitionList::pop_back
void pop_back(void)
Definition: transition.h:301
ipSc05231
long ipSc05231
Definition: atmdat_readin.cpp:84
ipxMg08303
long ipxMg08303
Definition: atmdat_readin.cpp:62
cddefines.h
qList
Definition: quantumstate.h:40
ipT324
long ipT324
Definition: atmdat_readin.cpp:45
ipSi2_2334
long ipSi2_2334
Definition: atmdat_readin.cpp:94
ipCALCIUM
const int ipCALCIUM
Definition: cddefines.h:324
t_struc::DenMass
realnum * DenMass
Definition: struc.h:49
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
ipSi10_606
long ipSi10_606
Definition: atmdat_readin.cpp:49
t_mole_global::num_total
int num_total
Definition: mole.h:314
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
ipT2140
long ipT2140
Definition: atmdat_readin.cpp:45
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_struc::volstr
realnum * volstr
Definition: struc.h:26
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
ipTO1025
long ipTO1025
Definition: atmdat_readin.cpp:51
ipT333
long ipT333
Definition: atmdat_readin.cpp:46
t_struc::histr
realnum * histr
Definition: struc.h:28
nint
long nint(double x)
Definition: cddefines.h:719
ipCl04117
long ipCl04117
Definition: atmdat_readin.cpp:78
ipN3_1752
long ipN3_1752
Definition: atmdat_readin.cpp:92
hyperfine.h
ipT146
long ipT146
Definition: atmdat_readin.cpp:52
FeIICreate
void FeIICreate(void)
Definition: atom_feii.cpp:219
ipT705
long ipT705
Definition: atmdat_readin.cpp:60
ipTAr13
long ipTAr13
Definition: atmdat_readin.cpp:79
ipTOI11
long ipTOI11
Definition: atmdat_readin.cpp:51
t_struc::xLyman_depth
realnum * xLyman_depth
Definition: struc.h:55
t_yield::fl_yield
realnum fl_yield[MEWE_FLUOR]
Definition: yield.h:37
t_struc::depth_last
realnum * depth_last
Definition: struc.h:57
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
irsl2ind
int irsl2ind(vector< t_BadnellLevel > &level, int irsl)
Definition: atmdat_readin.cpp:2155
ipTFe56
long ipTFe56
Definition: atmdat_readin.cpp:88
t_struc::pres_radiation_lines_curr
realnum * pres_radiation_lines_curr
Definition: struc.h:43
ipCaI4228
long ipCaI4228
Definition: atmdat_readin.cpp:81
t_struc::hden
realnum * hden
Definition: struc.h:45
ipTSi65
long ipTSi65
Definition: atmdat_readin.cpp:73
ipT8727
long ipT8727
Definition: atmdat_readin.cpp:38
ipC2_2328
long ipC2_2328
Definition: atmdat_readin.cpp:93
heavy.h
ipMgI2026
long ipMgI2026
Definition: atmdat_readin.cpp:59
t_BadnellLevel::lgAutoIonizing
bool lgAutoIonizing
Definition: atmdat_readin.cpp:125
ipFeI3457
long ipFeI3457
Definition: atmdat_readin.cpp:85
ipxMg52855
long ipxMg52855
Definition: atmdat_readin.cpp:61
t_BadnellLevel::S
int S
Definition: atmdat_readin.cpp:121
ipT7324
long ipT7324
Definition: atmdat_readin.cpp:83
ipT8498
long ipT8498
Definition: atmdat_readin.cpp:82
DynaCreateArrays
void DynaCreateArrays(void)
Definition: dynamics.cpp:1400
MEWE_FLUOR
#define MEWE_FLUOR
Definition: yield.h:10
ipTOI13
long ipTOI13
Definition: atmdat_readin.cpp:50
ipT2796
long ipT2796
Definition: atmdat_readin.cpp:59
MAX2
#define MAX2
Definition: cddefines.h:782
ipT374x
long ipT374x
Definition: atmdat_readin.cpp:44
IS_K_SHELL
@ IS_K_SHELL
Definition: atmdat_readin.cpp:115
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
ipTMg625
long ipTMg625
Definition: atmdat_readin.cpp:62
ipfsCl233
long ipfsCl233
Definition: atmdat_readin.cpp:77
LIMELM
const int LIMELM
Definition: cddefines.h:258
ipT1403
long ipT1403
Definition: atmdat_readin.cpp:70
RYD_INF
const UNUSED double RYD_INF
Definition: physconst.h:115
t_BadnellLevel::L
int L
Definition: atmdat_readin.cpp:122
ipc31175
long ipc31175
Definition: atmdat_readin.cpp:40
ipN3_1754
long ipN3_1754
Definition: atmdat_readin.cpp:92
ipC2_2324
long ipC2_2324
Definition: atmdat_readin.cpp:93
t_yield::nelem
int nelem(long n) const
Definition: yield.h:59
ipTAl550
long ipTAl550
Definition: atmdat_readin.cpp:68
ipS4_1417
long ipS4_1417
Definition: atmdat_readin.cpp:90
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:854
UTALines
TransitionList UTALines("UTALines", &AnonStates)
NLINEHS
#define NLINEHS
Definition: atmdat.h:124
ipFe1_24m
long ipFe1_24m
Definition: atmdat_readin.cpp:96
OP_Helike_NumPts
long **** OP_Helike_NumPts
Definition: atmdat.cpp:12
ipTSi35
long ipTSi35
Definition: atmdat_readin.cpp:73
ipKI7745
long ipKI7745
Definition: atmdat_readin.cpp:80
ipTAr8
long ipTAr8
Definition: atmdat_readin.cpp:80
IS_L1_SHELL
@ IS_L1_SHELL
Definition: atmdat_readin.cpp:115
ipT122
long ipT122
Definition: atmdat_readin.cpp:46
TauLines
TransitionList TauLines("TauLines", &AnonStates)
ipTFe26
long ipTFe26
Definition: atmdat_readin.cpp:87
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
t_struc::pressure
realnum * pressure
Definition: struc.h:33
fudge
double fudge(long int ipnt)
Definition: service.cpp:481
ipTSi521
long ipTSi521
Definition: atmdat_readin.cpp:72
ipSi1_130m
long ipSi1_130m
Definition: atmdat_readin.cpp:97
ipxK04598
long ipxK04598
Definition: atmdat_readin.cpp:80
ipT157
long ipT157
Definition: atmdat_readin.cpp:42
ipPRD
const int ipPRD
Definition: cddefines.h:290
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
ipAl6366
long ipAl6366
Definition: atmdat_readin.cpp:66
ipT835
long ipT835
Definition: atmdat_readin.cpp:48
ipTFe35
long ipTFe35
Definition: atmdat_readin.cpp:88
ipAr06453
long ipAr06453
Definition: atmdat_readin.cpp:80
ipFe1_54m
long ipFe1_54m
Definition: atmdat_readin.cpp:96
IS_TOP
@ IS_TOP
Definition: atmdat_readin.cpp:115
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:664
INTBIG
const long INTBIG
Definition: atmdat_readin.cpp:30
TransitionList
Definition: transition.h:274
ionbal.h
t_yield::ion
int ion(long n) const
Definition: yield.h:60
t_ionbal::lgInnerShell_Kisielius
bool lgInnerShell_Kisielius
Definition: ionbal.h:177
t_BadnellLevel
Definition: atmdat_readin.cpp:117
t_struc::heatstr
double * heatstr
Definition: struc.h:79
ipT1909
long ipT1909
Definition: atmdat_readin.cpp:38
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
ipT6363
long ipT6363
Definition: atmdat_readin.cpp:47
ipxMg72261
long ipxMg72261
Definition: atmdat_readin.cpp:61
t_MeweCoef::g
realnum g[210][4]
Definition: mewecoef.h:13
ipT8446
long ipT8446
Definition: atmdat_readin.cpp:50
t_struc::DenParticles
realnum * DenParticles
Definition: struc.h:47
physconst.h
ipTSi3
long ipTSi3
Definition: atmdat_readin.cpp:74
ipO4_1401
long ipO4_1401
Definition: atmdat_readin.cpp:91
HS_He1_Xsectn
double **** HS_He1_Xsectn
Definition: atmdat.cpp:8
t_struc::ednstr
realnum * ednstr
Definition: struc.h:30
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
t_ionbal::lgInnerShell_Gu06
bool lgInnerShell_Gu06
Definition: ionbal.h:180
ipFe1_111m
long ipFe1_111m
Definition: atmdat_readin.cpp:96
atmdat_outer_shell
void atmdat_outer_shell(long int iz, long int in, long int *imax, long int *ig0, long int *ig1)
Definition: atmdat_outer_shell.cpp:8
ipT88
long ipT88
Definition: atmdat_readin.cpp:55
ipT765
long ipT765
Definition: atmdat_readin.cpp:43
ipT1085
long ipT1085
Definition: atmdat_readin.cpp:42
fixit
void fixit(void)
Definition: service.cpp:991
ipT977
long ipT977
Definition: atmdat_readin.cpp:39
ipxMg71190
long ipxMg71190
Definition: atmdat_readin.cpp:61
ipT312
long ipT312
Definition: atmdat_readin.cpp:41
t_struc::GasPressure
realnum * GasPressure
Definition: struc.h:41
ipxK04154
long ipxK04154
Definition: atmdat_readin.cpp:81
t_struc::AccelGravity
realnum * AccelGravity
Definition: struc.h:39
ipFeI3021
long ipFeI3021
Definition: atmdat_readin.cpp:86
t_atmdat::lgStoutHybrid
bool lgStoutHybrid
Definition: atmdat.h:243
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
ipAl8370
long ipAl8370
Definition: atmdat_readin.cpp:67
ipT191
long ipT191
Definition: atmdat_readin.cpp:89
t_yield::n_elec_eject
long int n_elec_eject[30][30][7]
Definition: yield.h:26
taulines.h
ipT386
long ipT386
Definition: atmdat_readin.cpp:39
ipT1527
long ipT1527
Definition: atmdat_readin.cpp:71
ipT205
long ipT205
Definition: atmdat_readin.cpp:46
t_yield::frac_elec_eject
realnum frac_elec_eject[30][30][7][10]
Definition: yield.h:25
ipAl8575
long ipAl8575
Definition: atmdat_readin.cpp:66
OP_Helike_Xsectn
double ***** OP_Helike_Xsectn
Definition: atmdat.cpp:10
dprintf
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1009
ipO4_1405
long ipO4_1405
Definition: atmdat_readin.cpp:91
ipCo11527
long ipCo11527
Definition: atmdat_readin.cpp:89
ipT1863
long ipT1863
Definition: atmdat_readin.cpp:65
ipxNa6862
long ipxNa6862
Definition: atmdat_readin.cpp:58
geometry
t_geometry geometry
Definition: geometry.cpp:5
t_struc::AccelTotalOutward
realnum * AccelTotalOutward
Definition: struc.h:37
t_atmdat::ElecTemp
double ElecTemp[2][HS_NZ][NHSDIM]
Definition: atmdat.h:182
ipTAr9
long ipTAr9
Definition: atmdat_readin.cpp:79
ipT150
long ipT150
Definition: atmdat_readin.cpp:52
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
t_struc::o3str
realnum * o3str
Definition: struc.h:31
t_atmdat::Emiss
double Emiss[2][HS_NZ][NHSDIM][NHSDIM][NLINEHS]
Definition: atmdat.h:184
ipS4_1424
long ipS4_1424
Definition: atmdat_readin.cpp:90
ipTMg610
long ipTMg610
Definition: atmdat_readin.cpp:62
HyperfineCreate
void HyperfineCreate(void)
Definition: atom_hyperfine.cpp:392
ipT786
long ipT786
Definition: atmdat_readin.cpp:76
ipTSi4
long ipTSi4
Definition: atmdat_readin.cpp:74
iterations.h
ipSi1_68m
long ipSi1_68m
Definition: atmdat_readin.cpp:97
ipN3_1747
long ipN3_1747
Definition: atmdat_readin.cpp:92
ipTMg4
long ipTMg4
Definition: atmdat_readin.cpp:63
ipT1032
long ipT1032
Definition: atmdat_readin.cpp:54
ipT1256
long ipT1256
Definition: atmdat_readin.cpp:75
ipP0318
long ipP0318
Definition: atmdat_readin.cpp:75
ipT789
long ipT789
Definition: atmdat_readin.cpp:49
ipC2_2329
long ipC2_2329
Definition: atmdat_readin.cpp:93
ipMG_LIKE
const int ipMG_LIKE
Definition: iso.h:73
ipLI_LIKE
const int ipLI_LIKE
Definition: iso.h:64
ipT1666
long ipT1666
Definition: atmdat_readin.cpp:48
ipSi619
long ipSi619
Definition: atmdat_readin.cpp:71
t_BadnellLevel::config
string config
Definition: atmdat_readin.cpp:119
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
ipT57
long ipT57
Definition: atmdat_readin.cpp:47
ipCl04203
long ipCl04203
Definition: atmdat_readin.cpp:78
ipT639
long ipT639
Definition: atmdat_readin.cpp:67
t_iterations::iter_malloc
long int iter_malloc
Definition: iterations.h:29
ipT274
long ipT274
Definition: atmdat_readin.cpp:41
h2.h
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
ipxNe0676
long ipxNe0676
Definition: atmdat_readin.cpp:55
ipxNa0746
long ipxNa0746
Definition: atmdat_readin.cpp:58
t_struc::molecules
realnum ** molecules
Definition: struc.h:71
nSpecies
long int nSpecies
Definition: taulines.cpp:21
ipTNe13
long ipTNe13
Definition: atmdat_readin.cpp:55
ipT1214
long ipT1214
Definition: atmdat_readin.cpp:52
TransitionList::back
const TransitionProxy back(void)
Definition: transition.h:317
ipTNe16
long ipTNe16
Definition: atmdat_readin.cpp:56
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
ipT1548
long ipT1548
Definition: atmdat_readin.cpp:39
ipxK07319
long ipxK07319
Definition: atmdat_readin.cpp:81
ipFeI2966
long ipFeI2966
Definition: atmdat_readin.cpp:86
ipT1200
long ipT1200
Definition: atmdat_readin.cpp:44
TransitionList::reserve
void reserve(size_t newsize)
Definition: transition.h:289
NISO
const int NISO
Definition: cddefines.h:261
t_yield::init_yield
void init_yield()
Definition: atmdat_readin.cpp:1781
ipNI_pumpDirect
long ipNI_pumpDirect[NI_NDP]
Definition: atmdat_readin.cpp:98
t_atmdat::Density
double Density[2][HS_NZ][NHSDIM]
Definition: atmdat.h:181
ipT291
long ipT291
Definition: atmdat_readin.cpp:40
ipSi2_2344
long ipSi2_2344
Definition: atmdat_readin.cpp:94
ipT1855
long ipT1855
Definition: atmdat_readin.cpp:65
input.h
ipTNe14
long ipTNe14
Definition: atmdat_readin.cpp:56
ipFe1_35m
long ipFe1_35m
Definition: atmdat_readin.cpp:96
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
ipN3_1749
long ipN3_1749
Definition: atmdat_readin.cpp:92
ipT1239
long ipT1239
Definition: atmdat_readin.cpp:44
TransitionProxy::AddLoState
void AddLoState() const
Definition: transition.cpp:676
ipN3_1751
long ipN3_1751
Definition: atmdat_readin.cpp:92
ipT304
long ipT304
Definition: atmdat_readin.cpp:52
IS_L2_SHELL
@ IS_L2_SHELL
Definition: atmdat_readin.cpp:115
f_cutoff
const realnum f_cutoff
Definition: atmdat_readin.cpp:147
ipT5895
long ipT5895
Definition: atmdat_readin.cpp:57
atomfeii.h
IS_NONE
@ IS_NONE
Definition: atmdat_readin.cpp:115
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipT630
long ipT630
Definition: atmdat_readin.cpp:49
exc_type
exc_type
Definition: atmdat_readin.cpp:115
ipxMg51325
long ipxMg51325
Definition: atmdat_readin.cpp:60
ipTS1720
long ipTS1720
Definition: atmdat_readin.cpp:76
ipT1394
long ipT1394
Definition: atmdat_readin.cpp:70
ipTS34
long ipTS34
Definition: atmdat_readin.cpp:77
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
ipT6300
long ipT6300
Definition: atmdat_readin.cpp:47
ipT26
long ipT26
Definition: atmdat_readin.cpp:53
ipT610
long ipT610
Definition: atmdat_readin.cpp:41
ipSi2_2329
long ipSi2_2329
Definition: atmdat_readin.cpp:94
t_geometry::nend
long int * nend
Definition: geometry.h:80
ipT1550
long ipT1550
Definition: atmdat_readin.cpp:39
ipT8662
long ipT8662
Definition: atmdat_readin.cpp:83
ipfsCl214
long ipfsCl214
Definition: atmdat_readin.cpp:77
ipO4_1407
long ipO4_1407
Definition: atmdat_readin.cpp:91
t_yield::fl_energy
realnum fl_energy[MEWE_FLUOR]
Definition: yield.h:35
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17