cloudy trunk
Loading...
Searching...
No Matches
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! */
30const 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 */
37long 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;
90long 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!
115typedef enum { IS_NONE, IS_K_SHELL, IS_L1_SHELL, IS_L2_SHELL, IS_TOP } exc_type;
116
117struct 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
131STATIC void read_SH98_He1_cross_sections(void);
132
133STATIC void read_Helike_cross_sections(void);
134
135// read autoionization data from Badnell data file
136STATIC 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
142inline void InitTransition(const TransitionProxy& t);
143inline 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
147const realnum f_cutoff = 1.e-4f;
148
149void 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
1446STATIC 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
1591STATIC 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
1758t_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
1781void 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
1966STATIC 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
2148inline void InitTransition(const TransitionProxy& t)
2149{
2150 t.AddHiState();
2151 t.AddLoState();
2152 t.AddLine2Stack();
2153}
2154
2155inline 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}
OP_Helike_Xsectn
double ***** OP_Helike_Xsectn
Definition: atmdat.cpp:10
HS_He1_Xsectn
double **** HS_He1_Xsectn
Definition: atmdat.cpp:8
HS_He1_Energy
double **** HS_He1_Energy
Definition: atmdat.cpp:9
OP_Helike_NumPts
long **** OP_Helike_NumPts
Definition: atmdat.cpp:12
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
OP_Helike_Energy
double ***** OP_Helike_Energy
Definition: atmdat.cpp:11
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
NLINEHS
#define NLINEHS
Definition: atmdat.h:124
NHSDIM
#define NHSDIM
Definition: atmdat.h:123
NUM_HS98_DATA_POINTS
#define NUM_HS98_DATA_POINTS
Definition: atmdat.h:127
HS_NZ
#define HS_NZ
Definition: atmdat.h:125
atmdat_2phot_setSplineCoefs
void atmdat_2phot_setSplineCoefs()
Definition: atmdat_2photon.cpp:213
lines_setup
void lines_setup(void)
Definition: atmdat_lines_setup.cpp:19
ipT157
long ipT157
Definition: atmdat_readin.cpp:42
ipxMg51325
long ipxMg51325
Definition: atmdat_readin.cpp:60
atmdat_readin
void atmdat_readin(void)
Definition: atmdat_readin.cpp:149
ipTAr13
long ipTAr13
Definition: atmdat_readin.cpp:79
ipP0233
long ipP0233
Definition: atmdat_readin.cpp:74
ipTMg610
long ipTMg610
Definition: atmdat_readin.cpp:62
ipT8727
long ipT8727
Definition: atmdat_readin.cpp:38
ipCl1_11m
long ipCl1_11m
Definition: atmdat_readin.cpp:95
ipCo11527
long ipCo11527
Definition: atmdat_readin.cpp:89
ipTSi25
long ipTSi25
Definition: atmdat_readin.cpp:73
ipT770
long ipT770
Definition: atmdat_readin.cpp:54
ipTAl48
long ipTAl48
Definition: atmdat_readin.cpp:68
ipT1085
long ipT1085
Definition: atmdat_readin.cpp:42
ipT8498
long ipT8498
Definition: atmdat_readin.cpp:82
ipTNe14
long ipTNe14
Definition: atmdat_readin.cpp:56
ipC2_2329
long ipC2_2329
Definition: atmdat_readin.cpp:93
ipT386
long ipT386
Definition: atmdat_readin.cpp:39
ipAl6912
long ipAl6912
Definition: atmdat_readin.cpp:66
ipTS34
long ipTS34
Definition: atmdat_readin.cpp:77
ipT8662
long ipT8662
Definition: atmdat_readin.cpp:83
ipTNe13
long ipTNe13
Definition: atmdat_readin.cpp:55
ipTAr22
long ipTAr22
Definition: atmdat_readin.cpp:79
ipT977
long ipT977
Definition: atmdat_readin.cpp:39
ipSi2_2334
long ipSi2_2334
Definition: atmdat_readin.cpp:94
ipTNe24
long ipTNe24
Definition: atmdat_readin.cpp:56
ipTFe34
long ipTFe34
Definition: atmdat_readin.cpp:87
ipMgI2026
long ipMgI2026
Definition: atmdat_readin.cpp:59
ipfsMg790
long ipfsMg790
Definition: atmdat_readin.cpp:64
ipxNa0746
long ipxNa0746
Definition: atmdat_readin.cpp:58
ipT6363
long ipT6363
Definition: atmdat_readin.cpp:47
ipxK03462
long ipxK03462
Definition: atmdat_readin.cpp:80
ipT1661
long ipT1661
Definition: atmdat_readin.cpp:48
ipO4_1401
long ipO4_1401
Definition: atmdat_readin.cpp:91
ipSi10143
long ipSi10143
Definition: atmdat_readin.cpp:72
ipAl8370
long ipAl8370
Definition: atmdat_readin.cpp:67
ipAr06453
long ipAr06453
Definition: atmdat_readin.cpp:80
ipT324
long ipT324
Definition: atmdat_readin.cpp:45
ipT2140
long ipT2140
Definition: atmdat_readin.cpp:45
f_cutoff
const realnum f_cutoff
Definition: atmdat_readin.cpp:147
ipAl8575
long ipAl8575
Definition: atmdat_readin.cpp:66
ipTAl568
long ipTAl568
Definition: atmdat_readin.cpp:68
ipT205
long ipT205
Definition: atmdat_readin.cpp:46
ipO4_1405
long ipO4_1405
Definition: atmdat_readin.cpp:91
ipAl09204
long ipAl09204
Definition: atmdat_readin.cpp:67
ipT1666
long ipT1666
Definition: atmdat_readin.cpp:48
ipSi2_2336
long ipSi2_2336
Definition: atmdat_readin.cpp:94
INTBIG
const long INTBIG
Definition: atmdat_readin.cpp:30
ipC2_2327
long ipC2_2327
Definition: atmdat_readin.cpp:93
ipCl973
long ipCl973
Definition: atmdat_readin.cpp:78
ipNi1_7m
long ipNi1_7m
Definition: atmdat_readin.cpp:97
ipTSi65
long ipTSi65
Definition: atmdat_readin.cpp:73
ipO4_1397
long ipO4_1397
Definition: atmdat_readin.cpp:91
ipTFe26
long ipTFe26
Definition: atmdat_readin.cpp:87
ipTCa3
long ipTCa3
Definition: atmdat_readin.cpp:83
ipSi2_2329
long ipSi2_2329
Definition: atmdat_readin.cpp:94
ipT146
long ipT146
Definition: atmdat_readin.cpp:52
ipT8542
long ipT8542
Definition: atmdat_readin.cpp:82
ipT26
long ipT26
Definition: atmdat_readin.cpp:53
ipFe1_24m
long ipFe1_24m
Definition: atmdat_readin.cpp:96
ipT280
long ipT280
Definition: atmdat_readin.cpp:40
ipTNe36
long ipTNe36
Definition: atmdat_readin.cpp:56
ipT7324
long ipT7324
Definition: atmdat_readin.cpp:83
ipT1550
long ipT1550
Definition: atmdat_readin.cpp:39
ipT1194
long ipT1194
Definition: atmdat_readin.cpp:75
ipFe1_35m
long ipFe1_35m
Definition: atmdat_readin.cpp:96
ipAlI3957
long ipAlI3957
Definition: atmdat_readin.cpp:64
ipxMg52855
long ipxMg52855
Definition: atmdat_readin.cpp:61
ipO4_1407
long ipO4_1407
Definition: atmdat_readin.cpp:91
ipT1394
long ipT1394
Definition: atmdat_readin.cpp:70
ipT1207
long ipT1207
Definition: atmdat_readin.cpp:70
ipT780
long ipT780
Definition: atmdat_readin.cpp:54
ipSi1_68m
long ipSi1_68m
Definition: atmdat_readin.cpp:97
ipTSi41
long ipTSi41
Definition: atmdat_readin.cpp:72
ReadBadnellAIData
STATIC void ReadBadnellAIData(const string &fnam, long nelem, long ion, TransitionList &UTA, bitset< IS_TOP > Skip)
Definition: atmdat_readin.cpp:1966
ipCl04203
long ipCl04203
Definition: atmdat_readin.cpp:78
ipTO88
long ipTO88
Definition: atmdat_readin.cpp:53
ipT786
long ipT786
Definition: atmdat_readin.cpp:76
ipT1305
long ipT1305
Definition: atmdat_readin.cpp:71
ipC2_2328
long ipC2_2328
Definition: atmdat_readin.cpp:93
ipTS19
long ipTS19
Definition: atmdat_readin.cpp:77
ipC2_2324
long ipC2_2324
Definition: atmdat_readin.cpp:93
ipT990
long ipT990
Definition: atmdat_readin.cpp:43
ipC2_2325
long ipC2_2325
Definition: atmdat_readin.cpp:93
ipNi1_11m
long ipNi1_11m
Definition: atmdat_readin.cpp:97
ipT1656
long ipT1656
Definition: atmdat_readin.cpp:37
InitTransition
void InitTransition(const TransitionProxy &t)
Definition: atmdat_readin.cpp:2148
ipTOI46
long ipTOI46
Definition: atmdat_readin.cpp:51
ipc31175
long ipc31175
Definition: atmdat_readin.cpp:40
ipP713
long ipP713
Definition: atmdat_readin.cpp:75
ipT1037
long ipT1037
Definition: atmdat_readin.cpp:54
ipT88
long ipT88
Definition: atmdat_readin.cpp:55
ipN3_1752
long ipN3_1752
Definition: atmdat_readin.cpp:92
ipT5577
long ipT5577
Definition: atmdat_readin.cpp:47
ipAl6366
long ipAl6366
Definition: atmdat_readin.cpp:66
ipT1243
long ipT1243
Definition: atmdat_readin.cpp:43
ipT1863
long ipT1863
Definition: atmdat_readin.cpp:65
ipS1_25m
long ipS1_25m
Definition: atmdat_readin.cpp:95
ipCaI4228
long ipCaI4228
Definition: atmdat_readin.cpp:81
ipT304
long ipT304
Definition: atmdat_readin.cpp:52
ipTMg625
long ipTMg625
Definition: atmdat_readin.cpp:62
ipTMg14
long ipTMg14
Definition: atmdat_readin.cpp:63
ipT209
long ipT209
Definition: atmdat_readin.cpp:46
irsl2ind
int irsl2ind(vector< t_BadnellLevel > &level, int irsl)
Definition: atmdat_readin.cpp:2155
ipTAl550
long ipTAl550
Definition: atmdat_readin.cpp:68
ipT5895
long ipT5895
Definition: atmdat_readin.cpp:57
ipxNa6143
long ipxNa6143
Definition: atmdat_readin.cpp:58
ipT7291
long ipT7291
Definition: atmdat_readin.cpp:83
ipSi2_2350
long ipSi2_2350
Definition: atmdat_readin.cpp:94
ipT1122
long ipT1122
Definition: atmdat_readin.cpp:88
ipT765
long ipT765
Definition: atmdat_readin.cpp:43
ipN3_1754
long ipN3_1754
Definition: atmdat_readin.cpp:92
ipT1909
long ipT1909
Definition: atmdat_readin.cpp:38
ipT9830
long ipT9830
Definition: atmdat_readin.cpp:37
ipAl529
long ipAl529
Definition: atmdat_readin.cpp:66
ipP0318
long ipP0318
Definition: atmdat_readin.cpp:75
ipT374x
long ipT374x
Definition: atmdat_readin.cpp:44
ipTO1025
long ipTO1025
Definition: atmdat_readin.cpp:51
ipN3_1751
long ipN3_1751
Definition: atmdat_readin.cpp:92
ipT1486
long ipT1486
Definition: atmdat_readin.cpp:43
ipfsNa373
long ipfsNa373
Definition: atmdat_readin.cpp:57
ipfsMg755
long ipfsMg755
Definition: atmdat_readin.cpp:64
ipT639
long ipT639
Definition: atmdat_readin.cpp:67
ipTOI11
long ipTOI11
Definition: atmdat_readin.cpp:51
ipxMg52417
long ipxMg52417
Definition: atmdat_readin.cpp:60
ipP0260
long ipP0260
Definition: atmdat_readin.cpp:74
ipTAr8
long ipTAr8
Definition: atmdat_readin.cpp:80
ipxMg71190
long ipxMg71190
Definition: atmdat_readin.cpp:61
ipSii2215
long ipSii2215
Definition: atmdat_readin.cpp:69
ipT370
long ipT370
Definition: atmdat_readin.cpp:42
ipT191
long ipT191
Definition: atmdat_readin.cpp:89
ipxK04598
long ipxK04598
Definition: atmdat_readin.cpp:80
ipT270
long ipT270
Definition: atmdat_readin.cpp:41
ipN3_1747
long ipN3_1747
Definition: atmdat_readin.cpp:92
ipT122
long ipT122
Definition: atmdat_readin.cpp:46
ipTOI13
long ipTOI13
Definition: atmdat_readin.cpp:50
exc_type
exc_type
Definition: atmdat_readin.cpp:115
IS_L1_SHELL
@ IS_L1_SHELL
Definition: atmdat_readin.cpp:115
IS_NONE
@ IS_NONE
Definition: atmdat_readin.cpp:115
IS_L2_SHELL
@ IS_L2_SHELL
Definition: atmdat_readin.cpp:115
IS_K_SHELL
@ IS_K_SHELL
Definition: atmdat_readin.cpp:115
IS_TOP
@ IS_TOP
Definition: atmdat_readin.cpp:115
ipKI7745
long ipKI7745
Definition: atmdat_readin.cpp:80
ipT1256
long ipT1256
Definition: atmdat_readin.cpp:75
ipT58
long ipT58
Definition: atmdat_readin.cpp:63
ipT52
long ipT52
Definition: atmdat_readin.cpp:53
ipT2804
long ipT2804
Definition: atmdat_readin.cpp:59
ipfsCl214
long ipfsCl214
Definition: atmdat_readin.cpp:77
ipT3969
long ipT3969
Definition: atmdat_readin.cpp:82
ipxMg08303
long ipxMg08303
Definition: atmdat_readin.cpp:62
ipT4561
long ipT4561
Definition: atmdat_readin.cpp:60
ipxNe0676
long ipxNe0676
Definition: atmdat_readin.cpp:55
ipT1527
long ipT1527
Definition: atmdat_readin.cpp:71
ipO4_1400
long ipO4_1400
Definition: atmdat_readin.cpp:91
ipT8446
long ipT8446
Definition: atmdat_readin.cpp:50
ipSc05231
long ipSc05231
Definition: atmdat_readin.cpp:84
ipT1198
long ipT1198
Definition: atmdat_readin.cpp:76
ipT4368
long ipT4368
Definition: atmdat_readin.cpp:50
ipMgI2853
long ipMgI2853
Definition: atmdat_readin.cpp:58
ipT835
long ipT835
Definition: atmdat_readin.cpp:48
ipFe1_111m
long ipFe1_111m
Definition: atmdat_readin.cpp:96
ipfsNa421
long ipfsNa421
Definition: atmdat_readin.cpp:57
ipxMg72261
long ipxMg72261
Definition: atmdat_readin.cpp:61
ipTAr7
long ipTAr7
Definition: atmdat_readin.cpp:79
ipTS11
long ipTS11
Definition: atmdat_readin.cpp:77
ipT310
long ipT310
Definition: atmdat_readin.cpp:40
ipN3_1749
long ipN3_1749
Definition: atmdat_readin.cpp:92
ipfsCl233
long ipfsCl233
Definition: atmdat_readin.cpp:77
ipT333
long ipT333
Definition: atmdat_readin.cpp:46
ipT895
long ipT895
Definition: atmdat_readin.cpp:55
ipT150
long ipT150
Definition: atmdat_readin.cpp:52
ipT705
long ipT705
Definition: atmdat_readin.cpp:60
ipTMg4
long ipTMg4
Definition: atmdat_readin.cpp:63
ipTuv3
long ipTuv3
Definition: atmdat_readin.cpp:86
ipT834
long ipT834
Definition: atmdat_readin.cpp:48
ipT1304
long ipT1304
Definition: atmdat_readin.cpp:49
ipT1855
long ipT1855
Definition: atmdat_readin.cpp:65
ipTOI29
long ipTOI29
Definition: atmdat_readin.cpp:51
ipT1214
long ipT1214
Definition: atmdat_readin.cpp:52
ipFe1_54m
long ipFe1_54m
Definition: atmdat_readin.cpp:96
ipSc13264
long ipSc13264
Definition: atmdat_readin.cpp:84
ipT274
long ipT274
Definition: atmdat_readin.cpp:41
ipTFe46
long ipTFe46
Definition: atmdat_readin.cpp:88
ipS4_1407
long ipS4_1407
Definition: atmdat_readin.cpp:90
ipTFe56
long ipTFe56
Definition: atmdat_readin.cpp:88
ipT1403
long ipT1403
Definition: atmdat_readin.cpp:70
ipTNe16
long ipTNe16
Definition: atmdat_readin.cpp:56
ipT3934
long ipT3934
Definition: atmdat_readin.cpp:81
ipCl04117
long ipCl04117
Definition: atmdat_readin.cpp:78
ipTS1720
long ipTS1720
Definition: atmdat_readin.cpp:76
ipTSi521
long ipTSi521
Definition: atmdat_readin.cpp:72
ipT57
long ipT57
Definition: atmdat_readin.cpp:47
ipTSi4
long ipTSi4
Definition: atmdat_readin.cpp:74
ipTAr9
long ipTAr9
Definition: atmdat_readin.cpp:79
read_SH98_He1_cross_sections
STATIC void read_SH98_He1_cross_sections(void)
Definition: atmdat_readin.cpp:1446
ipT2796
long ipT2796
Definition: atmdat_readin.cpp:59
read_Helike_cross_sections
STATIC void read_Helike_cross_sections(void)
Definition: atmdat_readin.cpp:1591
ipT315
long ipT315
Definition: atmdat_readin.cpp:45
ipFeI2966
long ipFeI2966
Definition: atmdat_readin.cpp:86
ipT630
long ipT630
Definition: atmdat_readin.cpp:49
ipT1039
long ipT1039
Definition: atmdat_readin.cpp:50
ipS1_56m
long ipS1_56m
Definition: atmdat_readin.cpp:95
ipT1260
long ipT1260
Definition: atmdat_readin.cpp:71
ipT1335
long ipT1335
Definition: atmdat_readin.cpp:38
ipxK07319
long ipxK07319
Definition: atmdat_readin.cpp:81
ipFeI3021
long ipFeI3021
Definition: atmdat_readin.cpp:86
ipSii2518
long ipSii2518
Definition: atmdat_readin.cpp:68
ipSi619
long ipSi619
Definition: atmdat_readin.cpp:71
ipT1239
long ipT1239
Definition: atmdat_readin.cpp:44
ipSi2_2344
long ipSi2_2344
Definition: atmdat_readin.cpp:94
ipAlI3090
long ipAlI3090
Definition: atmdat_readin.cpp:64
ipT671
long ipT671
Definition: atmdat_readin.cpp:45
ipS4_1398
long ipS4_1398
Definition: atmdat_readin.cpp:90
ipFeI3729
long ipFeI3729
Definition: atmdat_readin.cpp:85
ipxK04154
long ipxK04154
Definition: atmdat_readin.cpp:81
ipNI_pumpIndirect
long ipNI_pumpIndirect
Definition: atmdat_readin.cpp:99
ipT1895
long ipT1895
Definition: atmdat_readin.cpp:70
ipSi10_606
long ipSi10_606
Definition: atmdat_readin.cpp:49
ipT2670
long ipT2670
Definition: atmdat_readin.cpp:65
ipTFe35
long ipTFe35
Definition: atmdat_readin.cpp:88
ipFeI3884
long ipFeI3884
Definition: atmdat_readin.cpp:85
ipS4_1405
long ipS4_1405
Definition: atmdat_readin.cpp:90
ipTSi35
long ipTSi35
Definition: atmdat_readin.cpp:73
ipT1548
long ipT1548
Definition: atmdat_readin.cpp:39
ipTSi3
long ipTSi3
Definition: atmdat_readin.cpp:74
ipxMg72569
long ipxMg72569
Definition: atmdat_readin.cpp:62
ipfsNa490
long ipfsNa490
Definition: atmdat_readin.cpp:57
ipT6300
long ipT6300
Definition: atmdat_readin.cpp:47
ipFeI3457
long ipFeI3457
Definition: atmdat_readin.cpp:85
ipT312
long ipT312
Definition: atmdat_readin.cpp:41
ipS4_1424
long ipS4_1424
Definition: atmdat_readin.cpp:90
ipT1200
long ipT1200
Definition: atmdat_readin.cpp:44
ipT291
long ipT291
Definition: atmdat_readin.cpp:40
ipTMg6
long ipTMg6
Definition: atmdat_readin.cpp:63
ipSi1_130m
long ipSi1_130m
Definition: atmdat_readin.cpp:97
ipxNa6862
long ipxNa6862
Definition: atmdat_readin.cpp:58
ipTSi499
long ipTSi499
Definition: atmdat_readin.cpp:72
ipTr48
long ipTr48
Definition: atmdat_readin.cpp:87
ipNI_pumpDirect
long ipNI_pumpDirect[NI_NDP]
Definition: atmdat_readin.cpp:98
ipT63
long ipT63
Definition: atmdat_readin.cpp:53
ipS4_1417
long ipS4_1417
Definition: atmdat_readin.cpp:90
ipT789
long ipT789
Definition: atmdat_readin.cpp:49
ipTFe16
long ipTFe16
Definition: atmdat_readin.cpp:87
ipT610
long ipT610
Definition: atmdat_readin.cpp:41
ipT1808
long ipT1808
Definition: atmdat_readin.cpp:69
ipT1032
long ipT1032
Definition: atmdat_readin.cpp:54
ipT374g
long ipT374g
Definition: atmdat_readin.cpp:44
FeIICreate
void FeIICreate(void)
Definition: atom_feii.cpp:219
HyperfineCreate
void HyperfineCreate(void)
Definition: atom_hyperfine.cpp:392
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
STATIC
#define STATIC
Definition: cddefines.h:97
ipIRON
const int ipIRON
Definition: cddefines.h:330
MIN2
#define MIN2
Definition: cddefines.h:761
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
fudge
double fudge(long int ipnt)
Definition: service.cpp:481
LIMELM
const int LIMELM
Definition: cddefines.h:258
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:307
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
dprintf
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1009
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
NISO
const int NISO
Definition: cddefines.h:261
strstr_s
const char * strstr_s(const char *haystack, const char *needle)
Definition: cddefines.h:1429
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
realnum
float realnum
Definition: cddefines.h:103
MAX2
#define MAX2
Definition: cddefines.h:782
nint
long nint(double x)
Definition: cddefines.h:719
ipCALCIUM
const int ipCALCIUM
Definition: cddefines.h:324
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:854
fixit
void fixit(void)
Definition: service.cpp:991
ipPRD
const int ipPRD
Definition: cddefines.h:290
nWindLine
long nWindLine
Definition: cdinit.cpp:19
EmissionProxy::iRedisFun
int & iRedisFun() const
Definition: emission.h:403
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
EmissionProxy::AutoIonizFrac
realnum & AutoIonizFrac() const
Definition: emission.h:583
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:175
TransitionList::resize
void resize(size_t newsize)
Definition: transition.h:285
TransitionList::push_back
void push_back(const TransitionProxy &tr)
Definition: transition.h:313
TransitionList::size
size_t size(void) const
Definition: transition.h:297
TransitionList::pop_back
void pop_back(void)
Definition: transition.h:301
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
TransitionList::reserve
void reserve(size_t newsize)
Definition: transition.h:289
TransitionList::back
const TransitionProxy back(void)
Definition: transition.h:317
TransitionList::begin
iterator begin(void)
Definition: transition.h:305
TransitionProxy::AddLine2Stack
void AddLine2Stack() const
Definition: transition.cpp:664
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
TransitionProxy::AddHiState
void AddHiState() const
Definition: transition.cpp:687
TransitionProxy::AddLoState
void AddLoState() const
Definition: transition.cpp:676
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
t_ionbal::lgInnerShell_Kisielius
bool lgInnerShell_Kisielius
Definition: ionbal.h:177
t_ionbal::lgInnerShell_Gu06
bool lgInnerShell_Gu06
Definition: ionbal.h:180
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
t_mole_global::num_total
int num_total
Definition: mole.h:314
t_yield::init_yield
void init_yield()
Definition: atmdat_readin.cpp:1781
t_yield::nfl_ion_emit
int nfl_ion_emit[MEWE_FLUOR]
Definition: yield.h:33
t_yield::n_elec_eject
long int n_elec_eject[30][30][7]
Definition: yield.h:26
t_yield::nfl_nshell
int nfl_nshell[MEWE_FLUOR]
Definition: yield.h:32
t_yield::t_yield
t_yield()
Definition: atmdat_readin.cpp:1758
t_yield::nfl_ion
int nfl_ion[MEWE_FLUOR]
Definition: yield.h:31
t_yield::nfl_nLine
int nfl_nLine[MEWE_FLUOR]
Definition: yield.h:34
t_yield::fl_yield
realnum fl_yield[MEWE_FLUOR]
Definition: yield.h:37
t_yield::nfl_lines
long int nfl_lines
Definition: yield.h:41
t_yield::fl_energy
realnum fl_energy[MEWE_FLUOR]
Definition: yield.h:35
t_yield::nfl_nelem
int nfl_nelem[MEWE_FLUOR]
Definition: yield.h:30
t_yield::frac_elec_eject
realnum frac_elec_eject[30][30][7][10]
Definition: yield.h:25
t_yield::nelem
int nelem(long n) const
Definition: yield.h:59
t_yield::lgKillAuger
bool lgKillAuger
Definition: yield.h:44
t_yield::ion
int ion(long n) const
Definition: yield.h:60
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
mode_r
const ios_base::openmode mode_r
Definition: cpu.h:212
UNUSED
#define UNUSED
Definition: cpu.h:14
dense
t_dense dense
Definition: dense.cpp:24
DynaCreateArrays
void DynaCreateArrays(void)
Definition: dynamics.cpp:1400
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
geometry
t_geometry geometry
Definition: geometry.cpp:5
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
ipLI_LIKE
const int ipLI_LIKE
Definition: iso.h:64
ipAL_LIKE
const int ipAL_LIKE
Definition: iso.h:74
ipMG_LIKE
const int ipMG_LIKE
Definition: iso.h:73
ipNA_LIKE
const int ipNA_LIKE
Definition: iso.h:72
iterations
t_iterations iterations
Definition: iterations.cpp:5
eina
double eina(double gf, double enercm, double gup)
Definition: lines_service.cpp:84
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:101
lines_service.h
MeweCoef
t_MeweCoef MeweCoef
Definition: mewecoef.cpp:5
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
PI4
UNUSED const double PI4
Definition: physconst.h:35
RYD_INF
UNUSED const double RYD_INF
Definition: physconst.h:115
EVRYD
UNUSED const double EVRYD
Definition: physconst.h:189
g
static double * g
Definition: species2.cpp:28
database_readin
void database_readin(void)
Definition: species.cpp:42
struc
t_struc struc
Definition: struc.cpp:6
t_BadnellLevel
Definition: atmdat_readin.cpp:118
t_BadnellLevel::S
int S
Definition: atmdat_readin.cpp:121
t_BadnellLevel::g
int g
Definition: atmdat_readin.cpp:123
t_BadnellLevel::energy
realnum energy
Definition: atmdat_readin.cpp:124
t_BadnellLevel::config
string config
Definition: atmdat_readin.cpp:119
t_BadnellLevel::t_BadnellLevel
t_BadnellLevel()
Definition: atmdat_readin.cpp:127
t_BadnellLevel::lgAutoIonizing
bool lgAutoIonizing
Definition: atmdat_readin.cpp:125
t_BadnellLevel::L
int L
Definition: atmdat_readin.cpp:122
t_BadnellLevel::irsl
int irsl
Definition: atmdat_readin.cpp:120
t_BadnellLevel::WhichShell
exc_type WhichShell
Definition: atmdat_readin.cpp:126
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_MeweCoef::g
realnum g[210][4]
Definition: mewecoef.h:13
t_atmdat::lgChiantiHybrid
bool lgChiantiHybrid
Definition: atmdat.h:233
t_atmdat::ElecTemp
double ElecTemp[2][HS_NZ][NHSDIM]
Definition: atmdat.h:182
t_atmdat::lgStoutHybrid
bool lgStoutHybrid
Definition: atmdat.h:243
t_atmdat::nDensity
long int nDensity[2][HS_NZ]
Definition: atmdat.h:188
t_atmdat::lgLamdaOn
bool lgLamdaOn
Definition: atmdat.h:239
t_atmdat::lgStoutOn
bool lgStoutOn
Definition: atmdat.h:241
t_atmdat::Emiss
double Emiss[2][HS_NZ][NHSDIM][NHSDIM][NLINEHS]
Definition: atmdat.h:184
t_atmdat::ntemp
long int ntemp[2][HS_NZ]
Definition: atmdat.h:188
t_atmdat::lgChiantiOn
bool lgChiantiOn
Definition: atmdat.h:231
t_atmdat::ncut
long int ncut[2][HS_NZ]
Definition: atmdat.h:188
t_atmdat::Density
double Density[2][HS_NZ][NHSDIM]
Definition: atmdat.h:181
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_geometry::nend
long int * nend
Definition: geometry.h:80
t_iterations::iter_malloc
long int iter_malloc
Definition: iterations.h:29
t_struc::depth_last
realnum * depth_last
Definition: struc.h:57
t_struc::AccelGravity
realnum * AccelGravity
Definition: struc.h:39
t_struc::depth
realnum * depth
Definition: struc.h:51
t_struc::xIonDense
realnum *** xIonDense
Definition: struc.h:64
t_struc::molecules
realnum ** molecules
Definition: struc.h:71
t_struc::drad_x_fillfac
realnum * drad_x_fillfac
Definition: struc.h:27
t_struc::nzlim
long int nzlim
Definition: struc.h:19
t_struc::windv
realnum * windv
Definition: struc.h:35
t_struc::heatstr
double * heatstr
Definition: struc.h:79
t_struc::drad_last
realnum * drad_last
Definition: struc.h:58
t_struc::testr
realnum * testr
Definition: struc.h:25
t_struc::volstr
realnum * volstr
Definition: struc.h:26
t_struc::drad
realnum * drad
Definition: struc.h:53
t_struc::histr
realnum * histr
Definition: struc.h:28
t_struc::gas_phase
realnum ** gas_phase
Definition: struc.h:75
t_struc::coolstr
double * coolstr
Definition: struc.h:78
t_struc::o3str
realnum * o3str
Definition: struc.h:31
t_struc::pressure
realnum * pressure
Definition: struc.h:33
t_struc::H2_abund
realnum * H2_abund
Definition: struc.h:72
t_struc::StatesElem
realnum **** StatesElem
Definition: struc.h:67
t_struc::ednstr
realnum * ednstr
Definition: struc.h:30
t_struc::hden
realnum * hden
Definition: struc.h:45
t_struc::GasPressure
realnum * GasPressure
Definition: struc.h:41
t_struc::xLyman_depth
realnum * xLyman_depth
Definition: struc.h:55
t_struc::pres_radiation_lines_curr
realnum * pres_radiation_lines_curr
Definition: struc.h:43
t_struc::DenMass
realnum * DenMass
Definition: struc.h:49
t_struc::AccelTotalOutward
realnum * AccelTotalOutward
Definition: struc.h:37
t_struc::hiistr
realnum * hiistr
Definition: struc.h:29
t_struc::DenParticles
realnum * DenParticles
Definition: struc.h:47
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
nSpecies
long int nSpecies
Definition: taulines.cpp:21
UTALines
TransitionList UTALines("UTALines", &AnonStates)
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
nUTA
long int nUTA
Definition: taulines.cpp:26
nLevel1
long int nLevel1
Definition: taulines.cpp:28
TauLines
TransitionList TauLines("TauLines", &AnonStates)
cs1_flag_lev2
realnum * cs1_flag_lev2
Definition: taulines.cpp:40
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
NI_NDP
const int NI_NDP
Definition: taulines.h:90
N
static const int N
Definition: thirdparty.cpp:2814
trace
t_trace trace
Definition: trace.cpp:5
MEWE_FLUOR
#define MEWE_FLUOR
Definition: yield.h:10
Generated by doxygen 1.9.5