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 */
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!
116
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
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
132
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
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");
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 */
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 */
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");
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");
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");
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 );
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) )
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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" );
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");
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 );
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" );
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 );
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 );
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
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 );
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 );
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");
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() );
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 */
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)
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" );
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 );
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 );
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" );
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 );
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
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
1430 // read cross sections for some he-like ions
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
1442
1443 return;
1444}
1445
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
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
1616 OP_Helike_Xsectn[ipHELIUM] = NULL;
1617 OP_Helike_Energy[ipHELIUM] = NULL;
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" );
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" );
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
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
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 );
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 );
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;
1919
1920 /* the ion stage for target, atom is 0 */
1921 nfl_ion[nfl_lines] = (int)temp[1]-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 */
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 */
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}
double ***** OP_Helike_Xsectn
Definition atmdat.cpp:10
double **** HS_He1_Xsectn
Definition atmdat.cpp:8
double **** HS_He1_Energy
Definition atmdat.cpp:9
long **** OP_Helike_NumPts
Definition atmdat.cpp:12
t_atmdat atmdat
Definition atmdat.cpp:6
double ***** OP_Helike_Energy
Definition atmdat.cpp:11
void atmdat_outer_shell(long int iz, long int in, long int *imax, long int *ig0, long int *ig1)
#define NLINEHS
Definition atmdat.h:124
#define NHSDIM
Definition atmdat.h:123
#define NUM_HS98_DATA_POINTS
Definition atmdat.h:127
#define HS_NZ
Definition atmdat.h:125
void atmdat_2phot_setSplineCoefs()
void lines_setup(void)
long ipT157
long ipxMg51325
void atmdat_readin(void)
long ipTAr13
long ipP0233
long ipTMg610
long ipT8727
long ipCl1_11m
long ipCo11527
long ipTSi25
long ipT770
long ipTAl48
long ipT1085
long ipT8498
long ipTNe14
long ipC2_2329
long ipT386
long ipAl6912
long ipTS34
long ipT8662
long ipTNe13
long ipTAr22
long ipT977
long ipSi2_2334
long ipTNe24
long ipTFe34
long ipMgI2026
long ipfsMg790
long ipxNa0746
long ipT6363
long ipxK03462
long ipT1661
long ipO4_1401
long ipSi10143
long ipAl8370
long ipAr06453
long ipT324
long ipT2140
const realnum f_cutoff
long ipAl8575
long ipTAl568
long ipT205
long ipO4_1405
long ipAl09204
long ipT1666
long ipSi2_2336
const long INTBIG
long ipC2_2327
long ipCl973
long ipNi1_7m
long ipTSi65
long ipO4_1397
long ipTFe26
long ipTCa3
long ipSi2_2329
long ipT146
long ipT8542
long ipT26
long ipFe1_24m
long ipT280
long ipTNe36
long ipT7324
long ipT1550
long ipT1194
long ipFe1_35m
long ipAlI3957
long ipxMg52855
long ipO4_1407
long ipT1394
long ipT1207
long ipT780
long ipSi1_68m
long ipTSi41
STATIC void ReadBadnellAIData(const string &fnam, long nelem, long ion, TransitionList &UTA, bitset< IS_TOP > Skip)
long ipCl04203
long ipTO88
long ipT786
long ipT1305
long ipC2_2328
long ipTS19
long ipC2_2324
long ipT990
long ipC2_2325
long ipNi1_11m
long ipT1656
void InitTransition(const TransitionProxy &t)
long ipTOI46
long ipc31175
long ipP713
long ipT1037
long ipT88
long ipN3_1752
long ipT5577
long ipAl6366
long ipT1243
long ipT1863
long ipS1_25m
long ipCaI4228
long ipT304
long ipTMg625
long ipTMg14
long ipT209
int irsl2ind(vector< t_BadnellLevel > &level, int irsl)
long ipTAl550
long ipT5895
long ipxNa6143
long ipT7291
long ipSi2_2350
long ipT1122
long ipT765
long ipN3_1754
long ipT1909
long ipT9830
long ipAl529
long ipP0318
long ipT374x
long ipTO1025
long ipN3_1751
long ipT1486
long ipfsNa373
long ipfsMg755
long ipT639
long ipTOI11
long ipxMg52417
long ipP0260
long ipTAr8
long ipxMg71190
long ipSii2215
long ipT370
long ipT191
long ipxK04598
long ipT270
long ipN3_1747
long ipT122
long ipTOI13
@ IS_L1_SHELL
@ IS_NONE
@ IS_L2_SHELL
@ IS_K_SHELL
@ IS_TOP
long ipKI7745
long ipT1256
long ipT58
long ipT52
long ipT2804
long ipfsCl214
long ipT3969
long ipxMg08303
long ipT4561
long ipxNe0676
long ipT1527
long ipO4_1400
long ipT8446
long ipSc05231
long ipT1198
long ipT4368
long ipMgI2853
long ipT835
long ipFe1_111m
long ipfsNa421
long ipxMg72261
long ipTAr7
long ipTS11
long ipT310
long ipN3_1749
long ipfsCl233
long ipT333
long ipT895
long ipT150
long ipT705
long ipTMg4
long ipTuv3
long ipT834
long ipT1304
long ipT1855
long ipTOI29
long ipT1214
long ipFe1_54m
long ipSc13264
long ipT274
long ipTFe46
long ipS4_1407
long ipTFe56
long ipT1403
long ipTNe16
long ipT3934
long ipCl04117
long ipTS1720
long ipTSi521
long ipT57
long ipTSi4
long ipTAr9
STATIC void read_SH98_He1_cross_sections(void)
long ipT2796
STATIC void read_Helike_cross_sections(void)
long ipT315
long ipFeI2966
long ipT630
long ipT1039
long ipS1_56m
long ipT1260
long ipT1335
long ipxK07319
long ipFeI3021
long ipSii2518
long ipSi619
long ipT1239
long ipSi2_2344
long ipAlI3090
long ipT671
long ipS4_1398
long ipFeI3729
long ipxK04154
long ipNI_pumpIndirect
long ipT1895
long ipSi10_606
long ipT2670
long ipTFe35
long ipFeI3884
long ipS4_1405
long ipTSi35
long ipT1548
long ipTSi3
long ipxMg72569
long ipfsNa490
long ipT6300
long ipFeI3457
long ipT312
long ipS4_1424
long ipT1200
long ipT291
long ipTMg6
long ipSi1_130m
long ipxNa6862
long ipTSi499
long ipTr48
long ipNI_pumpDirect[NI_NDP]
long ipT63
long ipS4_1417
long ipT789
long ipTFe16
long ipT610
long ipT1808
long ipT1032
long ipT374g
void FeIICreate(void)
void HyperfineCreate(void)
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipIRON
Definition cddefines.h:330
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
double fudge(long int ipnt)
Definition service.cpp:481
const int LIMELM
Definition cddefines.h:258
const int ipLITHIUM
Definition cddefines.h:307
#define EXIT_SUCCESS
Definition cddefines.h:138
int dprintf(FILE *fp, const char *format,...)
Definition service.cpp:1009
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const char * strstr_s(const char *haystack, const char *needle)
Definition cddefines.h:1429
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long nint(double x)
Definition cddefines.h:719
const int ipCALCIUM
Definition cddefines.h:324
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition cddefines.h:854
void fixit(void)
Definition service.cpp:991
const int ipPRD
Definition cddefines.h:290
long nWindLine
Definition cdinit.cpp:19
int & iRedisFun() const
Definition emission.h:403
realnum & gf() const
Definition emission.h:513
realnum & AutoIonizFrac() const
Definition emission.h:583
realnum & Aul() const
Definition emission.h:613
realnum & dampXvel() const
Definition emission.h:553
static t_yield & Inst()
Definition cddefines.h:175
TransitionProxy::iterator iterator
Definition transition.h:280
void push_back(const TransitionProxy &tr)
Definition transition.h:313
void pop_back(void)
Definition transition.h:301
const TransitionProxy back(void)
Definition transition.h:317
iterator begin(void)
Definition transition.h:305
void AddLine2Stack() const
realnum & WLAng() const
Definition transition.h:429
realnum & EnergyWN() const
Definition transition.h:438
qList::iterator Lo() const
Definition transition.h:392
void AddHiState() const
void AddLoState() const
qList::iterator Hi() const
Definition transition.h:396
EmissionList::reference Emis() const
Definition transition.h:408
void init_yield()
int nfl_ion_emit[MEWE_FLUOR]
Definition yield.h:33
long int n_elec_eject[30][30][7]
Definition yield.h:26
int nfl_nshell[MEWE_FLUOR]
Definition yield.h:32
int nfl_ion[MEWE_FLUOR]
Definition yield.h:31
int nfl_nLine[MEWE_FLUOR]
Definition yield.h:34
realnum fl_yield[MEWE_FLUOR]
Definition yield.h:37
long int nfl_lines
Definition yield.h:41
realnum fl_energy[MEWE_FLUOR]
Definition yield.h:35
int nfl_nelem[MEWE_FLUOR]
Definition yield.h:30
realnum frac_elec_eject[30][30][7][10]
Definition yield.h:25
int nelem(long n) const
Definition yield.h:59
bool lgKillAuger
Definition yield.h:44
int ion(long n) const
Definition yield.h:60
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const ios_base::openmode mode_r
Definition cpu.h:212
#define UNUSED
Definition cpu.h:14
t_dense dense
Definition dense.cpp:24
void DynaCreateArrays(void)
t_elementnames elementnames
t_geometry geometry
Definition geometry.cpp:5
t_Heavy Heavy
Definition heavy.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipLI_LIKE
Definition iso.h:64
const int ipAL_LIKE
Definition iso.h:74
const int ipMG_LIKE
Definition iso.h:73
const int ipNA_LIKE
Definition iso.h:72
t_iterations iterations
Definition iterations.cpp:5
double eina(double gf, double enercm, double gup)
double GetGF(double trans_prob, double enercm, double gup)
t_MeweCoef MeweCoef
Definition mewecoef.cpp:5
t_mole_global mole_global
Definition mole.cpp:6
#define S(I_, J_)
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double RYD_INF
Definition physconst.h:115
UNUSED const double EVRYD
Definition physconst.h:189
static double * g
Definition species2.cpp:28
void database_readin(void)
Definition species.cpp:42
t_struc struc
Definition struc.cpp:6
long int nSpecies
Definition taulines.cpp:21
TransitionList UTALines("UTALines", &AnonStates)
TransitionList TauLine2("TauLine2", &AnonStates)
long int nUTA
Definition taulines.cpp:26
long int nLevel1
Definition taulines.cpp:28
TransitionList TauLines("TauLines", &AnonStates)
realnum * cs1_flag_lev2
Definition taulines.cpp:40
vector< TransitionList > AllTransitions
Definition taulines.cpp:8
const int NI_NDP
Definition taulines.h:90
static const int N
t_trace trace
Definition trace.cpp:5
#define MEWE_FLUOR
Definition yield.h:10