cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_readin.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*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 "taulines.h"
7 #include "mewecoef.h"
8 #include "iterations.h"
9 #include "heavy.h"
10 #include "yield.h"
11 #include "trace.h"
12 #include "lines.h"
13 #include "struc.h"
14 #include "dynamics.h"
15 #include "elementnames.h"
16 #include "hyperfine.h"
17 #include "atmdat.h"
18 #include "iso.h"
19 #include "save.h"
20 #include "mole.h"
21 #include "two_photon.h"
22 #include "dense.h"
23 #include "lines_service.h"
24 
25 /* definition for whether level 2 lines are enabled, will be set to -1
26  * with no level2 command */
27 /*long nWindLine = NWINDDIM;*/
28 /*realnum TauLine2[NWINDDIM][NTA];*/
29 /*realnum **TauLine2;*/
30 
31 // NB NB - IS_TOP should always be the last entry of this enum!
33 
35 {
36  string config;
37  int irsl;
38  int S;
39  int L;
40  int g; // 2*J+1
41  realnum energy; // in cm^-1
44  t_BadnellLevel() : irsl(0), S(0), L(0), g(0), energy(0.f), lgAutoIonizing(false), WhichShell(IS_NONE) {}
45 };
46 
47 STATIC void init_struc();
49 STATIC void read_mewe_gbar();
51 
52 // read Hummer and Storey 98 He1 photoionization cross-section data
54 
56 
57 // read autoionization data from Badnell data file
58 STATIC void ReadBadnellAIData(const string& fnam, // filename containing the Badnell data
59  long nelem, // nelem is on C scale
60  long ion, // ion is on C scale
61  TransitionList& UTA, // UTA lines will be pushed on this stack
62  bitset<IS_TOP> Skip); // option to skip transitions from a particular shell
63 
64 STATIC void validate_magic_number_1arg( const char *chFilename, FILE *ioFile,
65  const long magicExp );
66 STATIC void validate_magic_number_3arg( const char *chFilename, FILE *ioFile,
67  const long yearExp, const long monthExp, const long dayExp );
68 STATIC void read_UTA_lines();
69 
70 // simple helper functions for ReadBadnellAIData
71 inline void InitTransition(const TransitionProxy& t);
72 inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl);
73 
74 // UTA lines below this absorption oscillator strength value will be ignored -
75 // F+13 paper plotted f not gf
76 const realnum f_cutoff = 1.e-4f;
77 
78 void atmdat_readin(void)
79 {
80  static bool lgFirstCall = true;
81 
82  DEBUG_ENTRY( "atmdat_readin()" );
83 
84  /* do nothing if not first call */
85  if( !lgFirstCall )
86  {
87  /* do not do anything, but make sure that number of zones has not increased */
88  bool lgTooBig = false;
89  for( long j=0; j < iterations.iter_malloc; j++ )
90  {
91  if( iterations.nend[j]>=struc.nzlim )
92  lgTooBig = true;
93  }
94  if( lgTooBig )
95  {
96  fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
97  fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
98  fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
99  fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
100  fprintf(ioQQQ," Sorry.\n");
102  }
103  return;
104  }
105 
106  lgFirstCall = false; /* do not reevaluate again */
107 
108  /* make sure that molecules have been initialized - this will fail
109  * if this routine is called before size of molecular network is known */
110  if( !mole_global.num_total )
111  {
112  /* mole_global.num_comole_calc can't be zero */
113  TotalInsanity();
114  }
115 
116  init_struc();
117 
118  /* allocate space for some arrays used by dynamics routines, and zero out vars */
119  DynaCreateArrays( );
120 
121  /*************************************************************
122  * *
123  * get the level 2 line, opacity project, data set *
124  * *
125  *************************************************************/
126 
127  /* nWindLine is initialized to the dimension of the vector when it is
128  * initialized in the definition at the start of this file.
129  * it is set to -1 with the "no level2" command, which
130  * stops us from trying to establish this vector */
131  if( nWindLine > 0 )
132  {
134  }
135 
136  /* the UTA line sets - trace will print summary of various data sources */
138  read_UTA_lines();
139 
140  /* read in data for the set of hyperfine structure lines, and allocate
141  * space for the transition HFLines[nHFLines] structure */
142  HyperfineCreate();
143 
144  /* Make sure that if hybrid is on, then Stout/Chianti are on */
146  {
147  TotalInsanity();
148  }
150  {
151  TotalInsanity();
152  }
153 
154  /* read in atomic and molecular models from third-party databases */
156  database_readin();
157  else
158  nSpecies = 0;
159 
160  /* initialize the large block of level 1 real lines, and OP level 2 lines */
161  lines_setup();
162 
163  /* mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat ========*/
164  /* read in g-bar data taken from
165  *>>refer all gbar Mewe, R., Gronenschild, E. H. B. M., van den Oord, G. H. J. 1985, A&AS, 62, 197 */
166  /* open file with Mewe coefficients */
167 
168  read_mewe_gbar();
169 
170  /* This is what remains of the t_yield initialization
171  * this should not be in the constructor of t_yield
172  * since it initializes a different struct */
173  for( long nelem=0; nelem < LIMELM; nelem++ )
174  for( long ion=0; ion < LIMELM; ion++ )
175  Heavy.nsShells[nelem][ion] = LONG_MAX;
176 
177  /* now read in all auger yields
178  * will do elements from li on up,
179  * skip any line starting with #
180  * this loop goes from lithium to Zn */
181  for( long nelem=2; nelem < LIMELM; nelem++ )
182  {
183  /* nelem is on the shifted C scale, so 2 is Li */
184  for( long ion=0; ion <= nelem; ion++ )
185  {
186  /* number of bound electrons, = atomic number for neutral */
187  long nelec = nelem - ion + 1;
188  /* one of dima's routines to determine the number of electrons
189  * for this species, nelem +1 to shift to physical number */
190  /* subroutine atmdat_outer_shell(iz,in,imax,ig0,ig1)
191  * iz - atomic number from 1 to 30 (integer)
192  * in - number of electrons from 1 to iz (integer)
193  * Output: imax - number of the outer shell
194  */
195  long imax,
196  ig0,
197  ig1;
198  atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);
199 
200  ASSERT( imax > 0 && imax <= 10 );
201 
202  /* nsShells[nelem][ion] is outer shell number for ion with nelec electrons
203  * on physics scale, with K shell being 1 */
204  Heavy.nsShells[nelem][ion] = imax;
205  }
206  }
207 
208  /*************************************************************
209  * *
210  * get the Auger electron yield data set *
211  * *
212  *************************************************************/
213 
215 
216  /****************************************************************
217  * *
218  * get the Hummer and Storey model case A, B results, these are *
219  * the two data files e1b.dat and e2b.dat, for H and He *
220  * *
221  ****************************************************************/
222 
224 
225  // read cross sections for neutral helium
227  // read cross sections for some he-like ions
229 
230  // set up spline coefficients for two-photon continua
232 
233  return;
234 }
235 
237 {
238  DEBUG_ENTRY( "init_struc()" );
239 
240  /* create space for the structure variables */
241  /* nzlim will be limit, and number allocated */
242  /* >>chng 01 jul 28, define this var, do all following mallocs */
243  for( long j=0; j < iterations.iter_malloc; j++ )
244  {
246  }
247 
248  /* sloppy, but add one extra for safely */
249  ++struc.nzlim;
250 
251  struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
252 
253  struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
254 
255  struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
256 
257  struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
258 
259  struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
260 
261  struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
262 
263  struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
264 
265  struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
266 
267  struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
268 
269  struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
270 
271  struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
272 
273  struc.AccelTotalOutward = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
274 
275  struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
276 
277  struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
278 
279  struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
280 
281  struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
282 
283  struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
284 
285  struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
286 
287  struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
288 
289  struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
290 
291  struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
292 
293  struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
294 
295  struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
296 
297  if (mole_global.num_calc != 0)
298  {
299  struc.molecules = ((realnum**)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum* )));
300  struc.molecules[0] = ((realnum*)MALLOC( (size_t)(struc.nzlim*mole_global.num_calc)*sizeof(realnum )));
301  for(long nz=1;nz<struc.nzlim;nz++)
302  {
304  }
305  }
306  else
307  {
308  struc.molecules = NULL;
309  }
310 
311 
312  struc.H2_abund = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
313 
314  /* create space for gas phase abundances array, first create space for the elements */
315  struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)struc.nzlim );
316  struc.gas_phase[0] =
317  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM*struc.nzlim ) );
318  for(long nz=1;nz<struc.nzlim;++nz)
319  {
320  struc.gas_phase[nz] = struc.gas_phase[nz-1]+LIMELM;
321  }
322 
323  /* create space for struc.xIonDense array, first create space for the zones */
324  struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) );
325  struc.xIonDense[0] = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim*LIMELM) );
326  struc.xIonDense[0][0] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim*LIMELM*(LIMELM+1)) );
327 
328  for( long nz=0; nz<struc.nzlim; ++nz)
329  {
330  if (nz > 0)
331  {
332  struc.xIonDense[nz] = struc.xIonDense[nz-1]+LIMELM;
333  struc.xIonDense[nz][0] = struc.xIonDense[nz-1][0]+LIMELM*(LIMELM+1);
334  }
335  for( long ipZ=0; ipZ<LIMELM;++ipZ )
336  {
337  if ( ipZ>0 )
338  struc.xIonDense[nz][ipZ] = struc.xIonDense[nz][ipZ-1]+(LIMELM+1);
339  /* space for the ionization stage */
340  }
341  }
342 
343  struc.StatesElem = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(struc.nzlim) );
344  for( long nz=0; nz<struc.nzlim; ++nz)
345  {
346  struc.StatesElem[nz] = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM) );
347  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
348  {
349  if( dense.lgElmtOn[nelem] )
350  {
351  struc.StatesElem[nz][nelem] = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+1) );
352  for( long ion=0; ion<nelem+1; ion++ )
353  {
354  long ipISO = nelem-ion;
355  if( ipISO < NISO )
356  {
357  struc.StatesElem[nz][nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)iso_sp[ipISO][nelem].numLevels_max);
358  }
359  else
360  {
361  fixit("for now, point non-iso ions to NULL");
362  struc.StatesElem[nz][nelem][ion] = NULL;
363  }
364  }
365  }
366  }
367  }
368 
369  /* some structure variables */
370  for( long i=0; i < struc.nzlim; i++ )
371  {
372  struc.testr[i] = 0.;
373  struc.volstr[i] = 0.;
374  struc.drad_x_fillfac[i] = 0.;
375  struc.histr[i] = 0.;
376  struc.hiistr[i] = 0.;
377  struc.ednstr[i] = 0.;
378  struc.o3str[i] = 0.;
379  struc.heatstr[i] = 0.;
380  struc.coolstr[i] = 0.;
381  struc.pressure[i] = 0.;
383  struc.GasPressure[i] = 0.;
384  struc.DenParticles[i] = 0.;
385  struc.depth[i] = 0.;
386  for( long mol=0;mol<mole_global.num_calc;mol++ )
387  {
388  struc.molecules[i][mol] = 0.;
389  }
390  struc.H2_abund[i] = 0.;
391  }
392 }
393 
395 {
396  DEBUG_ENTRY( "read_level2_lines()" );
397 
398  /* begin level2 level 2 wind block */
399  /* open file with level 2 line data */
400 
401  /* create the TauLine2 emline array */
403  AllTransitions.push_back(TauLine2);
404  cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum )));
405 
406  /* first initialize entire array to dangerously large negative numbers */
407  for( long i=0; i< nWindLine; ++i )
408  {
409  /* >>chng 99 jul 16, from setting all t[] to flt_max, to call for
410  * following, each member of structure set to own type of impossible value */
411  TauLine2[i].Junk();
412 
413  TauLine2[i].AddHiState();
414  TauLine2[i].AddLoState();
415  TauLine2[i].AddLine2Stack();
416  }
417 
418  if( trace.lgTrace )
419  fprintf( ioQQQ," atmdat_readin reading level2.dat\n");
420 
421  FILE *ioDATA = open_data( "level2.dat", "r" );
422 
423  validate_magic_number_3arg( "level2.dat", ioDATA, 9, 11, 18 );
424 
425  char chLine[FILENAME_PATH_LENGTH_2] = { 0 };
426 
427  /* now get the actual data */
428  for( long i=0; i < nWindLine; ++i )
429  {
430  do
431  {
432  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
433  {
434  fprintf( ioQQQ, " level2.dat error getting line %li\n", i );
436  }
437  /* skip any line starting with # */
438  } while ( chLine[0]=='#' );
439 
440  /* this must be double for sscanf to work below */
441  double tt[7] = { 0 };
442  /*printf("string is %s\n",chLine );*/
443  sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " ,
444  &tt[0] ,
445  &tt[1] ,
446  &tt[2] ,
447  &tt[3] ,
448  &tt[4] ,
449  &tt[5] ,
450  &tt[6] );
451  /* these are readjusted into their final form in the structure
452  * in routine lines_setup*/
453  (*TauLine2[i].Hi()).nelem() = (int)tt[0];
454  (*TauLine2[i].Hi()).IonStg() = (int)tt[1];
455  (*TauLine2[i].Lo()).g() = (realnum)tt[2];
456  (*TauLine2[i].Hi()).g() = (realnum)tt[3];
457  TauLine2[i].Emis().gf() = (realnum)tt[4];
458  TauLine2[i].EnergyWN() = (realnum)tt[5];
459  cs1_flag_lev2[i] = (realnum)tt[6];
460  }
461 
462  /* get magic number off last line */
463  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
464  {
465  fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
467  }
468  long magic;
469  sscanf( chLine , "%ld" , &magic );
470  if( 999 != magic )
471  {
472  fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n",
473  magic );
475  }
476  fclose( ioDATA );
477  if( trace.lgTrace )
478  fprintf( ioQQQ," reading level2.dat OK\n");
479 }
480 
482 {
483  DEBUG_ENTRY( "read_mewe_gbar()" );
484 
485  if( trace.lgTrace )
486  fprintf( ioQQQ," atmdat_readin reading mewe_gbar.dat\n");
487 
488  FILE *ioDATA = open_data( "mewe_gbar.dat", "r" );
489 
490  validate_magic_number_1arg( "mewe_gbar.dat", ioDATA, 9101 );
491 
492  char chLine[FILENAME_PATH_LENGTH_2] = { 0 };
493 
494  /* now get the actual data, indices are correct for c, in Fort went from 2 to 210 */
495  for( long i=1; i < 210; i++ )
496  {
497  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
498  {
499  fprintf( ioQQQ, " mewe_gbar.dat error getting line %li\n", i );
501  }
502  /*printf("%s\n",chLine);*/
503  double help[4];
504  sscanf( chLine, "%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] );
505  for( int l=0; l < 4; ++l )
506  MeweCoef.g[i][l] = (realnum)help[l];
507  }
508 
509  /* get magic number off last line */
510  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
511  {
512  fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
514  }
515 
516  long magic;
517  sscanf( chLine , "%ld" , &magic );
518 
519  if( magic != 9101 )
520  {
521  fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n",
522  magic );
524  }
525 
526  fclose( ioDATA );
527 
528  if( trace.lgTrace )
529  fprintf( ioQQQ," reading mewe_gbar.dat OK \n");
530 }
531 
533 {
534  DEBUG_ENTRY( "read_Hummer_Storey()" );
535 
536  /* now get Hummer and Storey case b data, loop over H, He */
537  /* >>chng 01 aug 08, generalized this to both case A and B, and all elements
538  * up to HS_NZ, now 8 */
539  for( long ipZ=0; ipZ<HS_NZ; ++ipZ )
540  {
541  /* don't do the minor elements, Li-B */
542  if( ipZ>1 && ipZ<5 ) continue;
543 
544  for( long iCase=0; iCase<2; ++iCase )
545  {
546  /* open Case A or B data */
547  /* >>chng 01 aug 08, add HS_ to start of file names to indicate Hummer Storey origin */
548  /* create file name for this charge
549  * first character after e is charge number for this element,
550  * then follows whether this is case A or case B data */
551  char chFilename[FILENAME_PATH_LENGTH_2] = { 0 };
552 
553  sprintf( chFilename, "HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ? 'a' : 'b' );
554 
555  if( trace.lgTrace )
556  fprintf( ioQQQ," atmdat_readin reading Hummer Storey emission file %s\n",chFilename );
557 
558  /* open the file */
559  FILE *ioDATA = open_data( chFilename, "r" );
560 
561  /* read in the number of temperatures and densities*/
562  {
563  int nread = fscanf( ioDATA, "%li %li ",
564  &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );
565 
566  if (nread != 2)
567  {
568  fprintf(ioQQQ, "atmdat_readin: bad input file format\n");
570  }
571  }
572 
573  /* check that ntemp and nDensity are below NHSDIM,
574  * set to 15 in atmdat_HS_caseB.h */
575  assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
576  assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM);
577 
578  /* loop reading in line emissivities for all temperatures*/
579  for( long ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
580  {
581  for( long ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
582  {
583  char cha;
584  long int junk, junk2;
585  {
586  int nread = fscanf( ioDATA, " %lf %li %lf %c %li %ld ",
587  &atmdat.Density[iCase][ipZ][ipDens], &junk ,
588  &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 ,
589  /* highest principal quantum number in table, usually 25 */
590  &atmdat.ncut[iCase][ipZ] );
591 
592  if (nread != 6)
593  {
594  fprintf(ioQQQ, "atmdat_readin: bad input file format\n");
596  }
597  }
598 
599  long ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
600  ASSERT( ne<=NLINEHS );
601  for( long j=0; j < ne; j++ )
602  {
603  int nread = fscanf( ioDATA, "%lf ",
604  &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
605 
606  if (nread != 1)
607  {
608  fprintf(ioQQQ, "atmdat_readin: bad input file format\n");
610  }
611  }
612  }
613  }
614 
615  /*this is end of read-in loop */
616  fclose(ioDATA);
617  if( trace.lgTrace )
618  fprintf( ioQQQ," reading %s OK\n", chFilename );
619 
620 # if 0
621  /* print list of densities and temperatures */
622  for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
623  {
624  fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]);
625  }
626  fprintf(ioQQQ,"\n");
627  for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
628  {
629  fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]);
630  }
631  fprintf(ioQQQ,"\n");
632 # endif
633  }
634  }
635 }
636 
638 {
639  DEBUG_ENTRY( "read_SH98_He1_cross_sections()" );
640 
641  FILE *ioDATA;
642  bool lgEOL;
643 
644  char chPath[FILENAME_PATH_LENGTH_2],
645  chDirectory[FILENAME_PATH_LENGTH_2],
646  chLine[FILENAME_PATH_LENGTH_2];
647 
648  const int ipNUM_FILES = 10;
649 
650  char chFileNames[ipNUM_FILES][10] = {
651  "p0202.3se",
652  "p0202.3po",
653  "p0202.3ge",
654  "p0202.3fo",
655  "p0202.3de",
656  "p0202.1se",
657  "p0202.1po",
658  "p0202.1ge",
659  "p0202.1fo",
660  "p0202.1de" };
661 
662  HS_He1_Xsectn = ((double****)MALLOC( (size_t)(26)*sizeof(double***)));
663  HS_He1_Energy = ((double****)MALLOC( (size_t)(26)*sizeof(double***)));
664 
665  // point these to NULL and use real quantum numbers rather than starting at 0
666  HS_He1_Xsectn[0] = NULL;
667  HS_He1_Energy[0] = NULL;
668 
669  for( long in = 1; in<=25; in++ )
670  {
671  // malloc n values of angular momentum, but not more than 5
672  HS_He1_Xsectn[in] = ((double***)MALLOC( (size_t)(MIN2(5,in))*sizeof(double**)));
673  HS_He1_Energy[in] = ((double***)MALLOC( (size_t)(MIN2(5,in))*sizeof(double**)));
674  for( long il = 0; il<MIN2(5,in); il++ )
675  {
676  // malloc two values of spin
677  HS_He1_Xsectn[in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
678  HS_He1_Energy[in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
679  HS_He1_Xsectn[in][il][0] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
680  HS_He1_Energy[in][il][0] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
681  HS_He1_Xsectn[in][il][1] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
682  HS_He1_Energy[in][il][1] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
683  }
684  }
685 
686 # ifdef _WIN32
687  strcpy( chDirectory, "sh98_he1\\pi\\" );
688 # else
689  strcpy( chDirectory, "sh98_he1/pi/" );
690 # endif
691 
692  //HS_He1_data[25][<=4][2]
693 
694  for( long ipFile=0; ipFile<ipNUM_FILES; ipFile++ )
695  {
696  long S, L, index, N=0;
697  long UNUSED P;
698 
699  strcpy( chPath, chDirectory );
700  strcat( chPath, chFileNames[ipFile] );
701  ioDATA = open_data( chPath, "r" );
702 
703  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
704  {
705  long i=1, s;
706  long i1, i2, i3;
707  long numDataPoints;
708 
709  // first line (read above in while) is not needed except that "0 0 0" marks EOF
710  i1 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
711  i2 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
712  i3 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
713  if( i1==0 && i2==0 && i3==0 )
714  break;
715 
716  // don't need next two lines in each set
717  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
718  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
719 
720  i=1;
721  // 4th line in each set identifies the quantum level
722  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
723  S = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
724  L = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
725  P = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
726  index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
727 
728  //indices start with unity
729  ASSERT( index >= 1 );
730 
731  // index is energy order in that series and is therefore related
732  // to principal quantum number. For triplet S must add one because
733  // series starts at n=2.
734  if( L==0 && S==3 )
735  N = L + index + 1;
736  else
737  N = L + index;
738 
739  // data go up to n=25
740  ASSERT( N<=25 );
741 
742  if( S==1 )
743  s=0;
744  else if( S==3 )
745  s=1;
746  else
747  TotalInsanity();
748 
749  i=1;
750  // 5th line in each set contains the number of energies
751  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
752  //first value is not needed
753  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
754  numDataPoints = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
755  // each set has exactly 811 data points, might as well assert this
756  if( numDataPoints != NUM_HS98_DATA_POINTS )
757  {
758  fprintf( ioQQQ, "wrong number of data points!\n" );
760  }
761 
762  // don't need 6th line either
763  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
764 
765  // now begin reading in lines
766  // there must be exactly numDataPoints ordered pairs,
767  // throw an assert for any deviation
768  for( long k=0; k<NUM_HS98_DATA_POINTS; k++ )
769  {
770  i=1;
771  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
772  HS_He1_Energy[N][L][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
773  HS_He1_Xsectn[N][L][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
774  }
775  }
776 
777  // we reached end of file, assert last quantum number was 25
778  ASSERT( N==25 );
779 
780  fclose( ioDATA );
781  }
782 
783  return;
784 }
785 
787 {
788  DEBUG_ENTRY( "read_Helike_cross_sections()" );
789 
790  FILE *ioDATA;
791  bool lgEOL;
792 
793  char chLine[FILENAME_PATH_LENGTH_2];
794  char chFileName[23] = "helike_pcs_topbase.dat";
795 
796  // the data only go as high as n=10
797  const int MaxN = 10;
798  long last_i1=0;
799 
800  // data will be used up to calcium (nelem=19)
801  // data exists up to iron but we will not use it because it has strong resonance
802  // features, but is nonetheless very nearly hydrogenic
803  OP_Helike_Xsectn = ((double*****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(double****)));
804  OP_Helike_Energy = ((double*****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(double****)));
805  OP_Helike_NumPts = ((long****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(long***)));
806 
807  // point these to NULL because we will never use them
811  OP_Helike_Xsectn[ipHELIUM] = NULL;
812  OP_Helike_Energy[ipHELIUM] = NULL;
814 
815  for( long nelem = ipLITHIUM; nelem<= ipCALCIUM; nelem++ )
816  {
817  // malloc principal quantum number
818  OP_Helike_Xsectn[nelem] = ((double****)MALLOC( (size_t)(MaxN+1)*sizeof(double***)));
819  OP_Helike_Energy[nelem] = ((double****)MALLOC( (size_t)(MaxN+1)*sizeof(double***)));
820  OP_Helike_NumPts[nelem] = ((long***)MALLOC( (size_t)(MaxN+1)*sizeof(long**)));
821 
822  for( long in = 1; in<=MaxN; in++ )
823  {
824  // malloc angular momentum
825  OP_Helike_Xsectn[nelem][in] = ((double***)MALLOC( (size_t)(in)*sizeof(double**)));
826  OP_Helike_Energy[nelem][in] = ((double***)MALLOC( (size_t)(in)*sizeof(double**)));
827  OP_Helike_NumPts[nelem][in] = ((long**)MALLOC( (size_t)(in)*sizeof(long*)));
828  for( long il = 0; il<in; il++ )
829  {
830  // malloc two values of spin
831  OP_Helike_Xsectn[nelem][in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
832  OP_Helike_Energy[nelem][in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
833  OP_Helike_NumPts[nelem][in][il] = ((long*)MALLOC( (size_t)(2)*sizeof(long)));
834  // point these to NULL for now, won't know how many we need until we begin parsing
835  OP_Helike_Xsectn[nelem][in][il][0] = NULL;
836  OP_Helike_Energy[nelem][in][il][0] = NULL;
837  OP_Helike_NumPts[nelem][in][il][0] = 0;
838  OP_Helike_Xsectn[nelem][in][il][1] = NULL;
839  OP_Helike_Energy[nelem][in][il][1] = NULL;
840  OP_Helike_NumPts[nelem][in][il][1] = 0;
841  }
842  }
843  }
844 
845  ioDATA = open_data( chFileName, "r" );
846 
847  // Header looks like this:
848  // ================================================
849  // I NZ NE ISLP ILV E(RYD) NP
850  // ================================================
851  // 1 3 2 100 1 -5.53159E+00 55
852 
853  // so skip the first three lines.
854  for( long i=0; i<3; i++)
855  {
856  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
857  {
858  fprintf( ioQQQ,"PROBLEM corruption in TOPbase Helike pcs datafile.\nSorry\n" );
860  }
861  }
862 
863  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
864  {
865  long i=1;
866  long n, l, s;
867  long i1, i2, i3, i4, i5, i7;
868  double i6;
869 
870  i1 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
871  i2 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
872  i3 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
873  i4 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
874  i5 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
875  i6 = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
876  i7 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
877 
878  if( lgEOL )
879  {
880  fprintf( ioQQQ,"PROBLEM corruption in TOPbase Helike pcs datafile.\nSorry\n" );
882  }
883 
884  // this marks end of data.
885  if( i1==i2 && i1==i3 && i1==i4 && i1==i5 && i1==i7 && i1==-1 )
886  break;
887 
888  // first parameter is level index (overall, not just for series or charge)
889  // check that it is as expected
890  if( ! ( i1 > 0 && i1 == (last_i1 + 1) && i1 <= 795 ) )
891  {
892  fprintf( ioQQQ, "syntax error found in %s\n", chFileName );
894  }
895  last_i1 = i1;
896  // second parameter is nuclear charge
897  ASSERT( (i2-1)<=ipCALCIUM && (i2-1)>=ipLITHIUM );
898  // third parameter must be 2 for helike.
899  ASSERT( i3==2 );
900  // fourth parameter is (2S+1)*100+L*10+P
901  ASSERT( i4>=100 && i4<400 );
902  if( i4 >= 300 )
903  s=1;
904  else
905  s=0;
906  l = (i4 - (2*s+1)*100)/10;
907  // data only goes up to l=2
908  ASSERT( l<=2 );
909  // fifth is index in the series, related to principal quantum number
910  ASSERT( i5>=1 && i5<=10 );
911  if( s==1 && l==0 )
912  n = i5 + 1;
913  else
914  n = i5 + l;
915  ASSERT( l<=MaxN );
916  // sixth is threshhold energy, don't need but assert negative
917  // \todo 3 save this and renorm cross-section with ratio of actual to recorded Eth?
918  if( i6 >= 0. )
919  {
920  fprintf( ioQQQ, "invalid threshold energy in %s\n", chFileName );
922  }
923  // seventh parameter is number of data points, can be zero, use to MALLOC otherwise
924  OP_Helike_NumPts[i2-1][n][l][s] = i7;
925  if( i7==0 )
926  continue;
927 
928  ASSERT( i7 > 0 );
929  ASSERT( OP_Helike_Xsectn[i2-1][n][l][s]==NULL );
930  ASSERT( OP_Helike_Energy[i2-1][n][l][s]==NULL );
931  OP_Helike_Xsectn[i2-1][n][l][s] = ((double*)MALLOC( (size_t)(i7)*sizeof(double)));
932  OP_Helike_Energy[i2-1][n][l][s] = ((double*)MALLOC( (size_t)(i7)*sizeof(double)));
933 
934  // now begin reading in lines
935  // there must be exactly i7 ordered pairs,
936  for( long k=0; k<i7; k++ )
937  {
938  i=1;
939  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
940  OP_Helike_Energy[i2-1][n][l][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
941  OP_Helike_Xsectn[i2-1][n][l][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
942 
943  // make sure data is well-behaved
944  if( k > 0 )
945  {
946  ASSERT( OP_Helike_Energy[i2-1][n][l][s][k] > OP_Helike_Energy[i2-1][n][l][s][k-1] );
947  ASSERT( OP_Helike_Xsectn[i2-1][n][l][s][k] < OP_Helike_Xsectn[i2-1][n][l][s][k-1] );
948  }
949 
950  // try to read one more item off the line and verify that lgEOL is true
951  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
952  ASSERT( lgEOL );
953  }
954  }
955 
956  fclose( ioDATA );
957 
958  return;
959 }
960 
962 {
963  /* preset two arrays that will hold auger data
964  * set to very large values so that code will blow
965  * if not set properly below */
966  for( int nelem=0; nelem < LIMELM; nelem++ )
967  {
968  for( int ion=0; ion < LIMELM; ion++ )
969  {
970  for( int ns=0; ns < 7; ns++ )
971  {
972  n_elec_eject[nelem][ion][ns] = LONG_MAX;
973  for( int nelec=0; nelec < 10; nelec++ )
974  {
975  frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
976  }
977  }
978  }
979  }
980 
981  lgKillAuger = false;
982 }
983 
985 {
986  char chLine[FILENAME_PATH_LENGTH_2];
987  const char* chFilename;
988 
989  /* following is double for sscanf to work */
990  double temp[14];
991 
992  DEBUG_ENTRY( "init_yield()" );
993 
994  /*************************************************************
995  * *
996  * get the Auger electron yield data set *
997  * *
998  *************************************************************/
999 
1000  /* NB NB -- This test of Heavy.nsShells remains here to assure
1001  * that it contains meaningful values since needed below, once
1002  * t_Heavy is a Singleton, it can be removed !!! */
1003  ASSERT( Heavy.nsShells[2][0] > 0 );
1004 
1005  /* hydrogen and helium will not be done below, so set yields here*/
1006  n_elec_eject[ipHYDROGEN][0][0] = 1;
1007  n_elec_eject[ipHELIUM][0][0] = 1;
1008  n_elec_eject[ipHELIUM][1][0] = 1;
1009 
1010  frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
1011  frac_elec_eject[ipHELIUM][0][0][0] = 1;
1012  frac_elec_eject[ipHELIUM][1][0][0] = 1;
1013 
1014  /* open file auger.dat, yield data file that came from
1015  * >>refer all auger Kaastra, J. S., & Mewe, R. 1993, A&AS, 97, 443-482 */
1016  chFilename = "mewe_nelectron.dat";
1017 
1018  if( trace.lgTrace )
1019  fprintf( ioQQQ, " init_yield reading %s\n", chFilename );
1020 
1021  FILE *ioDATA;
1022 
1023  /* open the file */
1024  ioDATA = open_data( chFilename, "r" );
1025 
1026  /* now read in all auger yields
1027  * will do elements from li on up,
1028  * skip any line starting with #
1029  * this loop goes from lithium to Zn */
1030  for( int nelem=2; nelem < LIMELM; nelem++ )
1031  {
1032  /* nelem is on the shifted C scale, so 2 is Li */
1033  for( int ion=0; ion <= nelem; ion++ )
1034  {
1035  for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1036  {
1037  char ch1 = '#';
1038  /* the * is the old comment char, accept it, but really want # */
1039  while( ch1 == '#' || ch1 == '*' )
1040  {
1041  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1042  {
1043  fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns );
1045  }
1046  ch1 = chLine[0];
1047  }
1048 
1049  if( sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
1050  &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1051  &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
1052  &temp[10],&temp[11],&temp[12],&temp[13] ) != 14 )
1053  {
1054  fprintf( ioQQQ, "failed to read correct number of arguments in %s\n",
1055  chFilename );
1057  }
1058  n_elec_eject[nelem][ion][ns] = (long int)temp[3];
1059 
1060  ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 );
1061 
1062  /* can pop off up to 10 auger electrons, these are the probabilities*/
1063  for( int j=0; j < 10; j++ )
1064  {
1065  frac_elec_eject[nelem][ion][ns][j] = (realnum)temp[j+4];
1066  ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. );
1067  }
1068  }
1069  }
1070  /* activate this print statement to get yield for k-shell of all atoms */
1071  /*fprintf(ioQQQ,"yyyield\t%li\t%.4e\n", nelem+1, frac_elec_eject[nelem][0][0][0] );*/
1072  }
1073 
1074  fclose( ioDATA );
1075 
1076  if( trace.lgTrace )
1077  {
1078  /* this is set with no auger command */
1079  if( lgKillAuger )
1080  fprintf( ioQQQ, " Auger yields will be killed.\n");
1081  fprintf( ioQQQ, " reading %s OK\n", chFilename );
1082  }
1083 
1084  /* open file mewe_fluor.dat, yield data file that came from
1085  * >>refer all auger Kaastra, J. S., & Mewe, R. 1993, 97, 443-482 */
1086  chFilename = "mewe_fluor.dat";
1087 
1088  if( trace.lgTrace )
1089  fprintf( ioQQQ, " init_yield reading %s\n", chFilename );
1090 
1091  /* open the file */
1092  ioDATA = open_data( chFilename, "r" );
1093 
1094  /* now read in all auger yields
1095  * will do elements from li on up,
1096  * skip any line starting with # */
1097  do
1098  {
1099  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1100  {
1101  fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 );
1103  }
1104  }
1105  while( chLine[0] == '#' );
1106 
1107  bool lgEOL = false;
1108 
1109  nfl_lines = 0;
1110  do
1111  {
1112  const int NKM = 10;
1113  int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
1114  int nAuger;
1115 
1116  if( nfl_lines >= MEWE_FLUOR )
1117  TotalInsanity();
1118 
1119  /*printf("string is %s\n",chLine );*/
1120  sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf",
1121  &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1122  &temp[5], &temp[6] );
1123 
1124  /* the atomic number, C is 5 */
1125  nfl_nelem[nfl_lines] = (int)temp[0]-1;
1126  ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM );
1127 
1128  /* the ion stage for target, atom is 0 */
1129  nfl_ion[nfl_lines] = (int)temp[1]-1;
1131 
1132  /* the target's shell */
1133  nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1];
1134  ASSERT( nfl_nshell[nfl_lines] >= 0 &&
1135  /* nsShells is shell number, where K is 1 */
1137 
1138  /* this is the number of Auger electrons ejected */
1139  nAuger = (int)temp[3];
1140  /* so this is the spectrum of the photons that are emitted */
1141  nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1;
1142  /* must be gt 0 since at least photoelectron comes off */
1144 
1145  /* this is the type of line as defined in their paper */
1146  nfl_nLine[nfl_lines] = (int)temp[4];
1147 
1148  /* energy in Ryd */
1149  fl_energy[nfl_lines] = (realnum)temp[5] / (realnum)EVRYD;
1150  ASSERT( fl_energy[nfl_lines] > 0. );
1151 
1152  /* fluor yield */
1153  fl_yield[nfl_lines] = (realnum)temp[6];
1154  /* NB cannot assert <=1 since data file has yields around 1.3 - 1.4 */
1155  ASSERT( fl_yield[nfl_lines] >= 0 );
1156 
1157  ++nfl_lines;
1158 
1159  do
1160  {
1161  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1162  lgEOL = true;
1163  }
1164  while( chLine[0]=='#' && !lgEOL );
1165  }
1166  while( !lgEOL );
1167 
1168  fclose( ioDATA );
1169  if( trace.lgTrace )
1170  fprintf( ioQQQ, " reading %s OK\n", chFilename );
1171 }
1172 
1173 // read autoionization data from Badnell data file
1174 STATIC void ReadBadnellAIData(const string& fnam, // filename containing the Badnell data
1175  long nelem, // nelem is on C scale
1176  long ion, // ion is on C scale
1177  TransitionList& UTA, // UTA lines will be pushed on this stack
1178  bitset<IS_TOP> Skip) // option to skip transitions from a particular shell
1179 {
1180  DEBUG_ENTRY( "ReadBadnellAIData()" );
1181 
1182  if( trace.lgTrace )
1183  fprintf( ioQQQ," ReadBadnellAIData reading %s\n", fnam.c_str() );
1184 
1185  fstream ioDATA;
1186  open_data( ioDATA, fnam.c_str(), mode_r );
1187 
1188  string line;
1189  getline( ioDATA, line );
1190  ASSERT( line.substr(0,4) == "SEQ=" );
1191  getline( ioDATA, line );
1192  getline( ioDATA, line );
1193  // we don't need the parent level data, so we will skip it...
1194  ASSERT( line.substr(3,21) == "PARENT LEVEL INDEXING" );
1195  int nParent;
1196  istringstream iss( line.substr(65,4) );
1197  iss >> nParent;
1198  // data lines containing data for all levels of the parent ion will span nMulti lines
1199  int nMulti = (nParent+5)/6;
1200  for( int i=0; i < nParent+5; ++i )
1201  getline( ioDATA, line );
1202 
1203  // here starts the header for the level data we need
1204  ASSERT( line.substr(3,26) == "IC RESOLVED LEVEL INDEXING" );
1205  int nLevel;
1206  istringstream iss2( line.substr(63,6) );
1207  iss2 >> nLevel;
1208  // skip rest of header
1209  for( int i=0; i < 3; ++i )
1210  getline( ioDATA, line );
1211 
1212  // now get the level data
1213  vector<t_BadnellLevel> level( nLevel );
1214  for( int i=0; i < nLevel; ++i )
1215  {
1216  getline( ioDATA, line );
1217  istringstream iss3( line );
1218  int indx, irsl;
1219  iss3 >> indx >> irsl;
1220  level[indx-1].irsl = irsl;
1221  level[indx-1].config = line.substr(16,20);
1222  istringstream iss4( line.substr(37,1) );
1223  iss4 >> level[indx-1].S;
1224  istringstream iss5( line.substr(39,1) );
1225  iss5 >> level[indx-1].L;
1226  istringstream iss6( line.substr(41,4) );
1227  double J;
1228  iss6 >> J;
1229  level[indx-1].g = nint(2.*J + 1.);
1230  istringstream iss7( line.substr(46,11) );
1231  iss7 >> level[indx-1].energy;
1232 
1233  // which inner shell has been excited?
1234  level[indx-1].lgAutoIonizing = ( line[57] == '*' );
1235  if( level[indx-1].lgAutoIonizing )
1236  {
1237  if( level[indx-1].config.find( "1S1" ) != string::npos )
1238  level[indx-1].WhichShell = IS_K_SHELL;
1239  else if( level[indx-1].config.find( "2S1" ) != string::npos )
1240  level[indx-1].WhichShell = IS_L1_SHELL;
1241  else if( level[indx-1].config.find( "2P5" ) != string::npos )
1242  level[indx-1].WhichShell = IS_L2_SHELL;
1243  else
1244  TotalInsanity();
1245  }
1246  else
1247  {
1248  level[indx-1].WhichShell = IS_NONE;
1249  }
1250  }
1251 
1252  // levels are done, now move on to the lines
1253  // first search for start of the header
1254  while( getline( ioDATA, line ) )
1255  {
1256  if( line.find( "IRSL IRSL" ) != string::npos )
1257  break;
1258  }
1259  // skip rest of the header
1260  for( int i=0; i < nMulti-1; ++i )
1261  getline( ioDATA, line );
1262 
1263  // blank line that will be pushed on the UTA line stack
1264  qList BlankStates("BlankStates",1);
1265  TransitionList BlankList("BlankList",&BlankStates,1);
1266  TransitionList::iterator BlankLine = BlankList.begin();
1267  (*BlankLine).Junk();
1268 
1269  // start reading the line data
1270  while( getline( ioDATA, line ) )
1271  {
1272  // have we reached the end of the line data?
1273  if( line.size() < 10 )
1274  break;
1275 
1276  // test if there is an autoionization rate on this line
1277  // if not, it only contains radiative data and we skip it...
1278  if( line.size() < 50 )
1279  continue;
1280 
1281  // there may be an asterisk here; wipe it out since we don't need it
1282  line[19] = ' ';
1283 
1284  int irsl_lo, irsl_hi, dum;
1285  double edif, Bij, Rji, Aai;
1286  istringstream iss8( line );
1287  // irsl_lo: index for lower level of transition
1288  // irsl_hi: index for upper level of transition
1289  // edif: energy difference between levels in Ryd
1290  // Bij: UPWARD Einstein A for transition irsl_lo -> irsl_hi
1291  // Rji: sum of Aji for all radiative transitions to lower levels
1292  // Aai: autoionization rate from level irsl_hi
1293  iss8 >> irsl_lo >> irsl_hi >> dum >> dum >> edif >> Bij >> Rji >> Aai;
1294  // ind_lo and ind_hi are on C scale
1295  int ind_lo = irsl2ind( level, irsl_lo );
1296  int ind_hi = irsl2ind( level, irsl_hi );
1297  ASSERT( level[ind_hi].lgAutoIonizing );
1298  // skip rest of partial autoionization rates
1299  for( int i=0; i < nMulti-1; ++i )
1300  getline( ioDATA, line );
1301  // skip this transition if it does not originate from ground level
1302  // or if the user requested to skip excitations from this inner shell
1303  if( ind_lo == 0 && !Skip[level[ind_hi].WhichShell] )
1304  {
1305  UTA.push_back( *BlankLine );
1306  InitTransition( UTA.back() );
1307 
1308  // t_emission has nelem and ion on fortran scale...
1309  (*UTA.back().Hi()).nelem() = nelem+1;
1310  (*UTA.back().Hi()).IonStg() = ion+1;
1311 
1312  (*UTA.back().Hi()).g() = (realnum)level[ind_hi].g;
1313  (*UTA.back().Lo()).g() = (realnum)level[ind_lo].g;
1314 
1315  double WavNum = edif*RYD_INF;
1316 
1317  /* wavelength in Angstroms */
1318  UTA.back().WLAng() = (realnum)(1e8/WavNum);
1319  UTA.back().EnergyWN() = (realnum)WavNum;
1320 
1321  /* store branching ratio for autoionization */
1322  double frac_ioniz = Aai/(Rji + Aai);
1323  ASSERT( frac_ioniz >= 0. && frac_ioniz <= 1. );
1324  UTA.back().Emis().AutoIonizFrac() = (realnum)frac_ioniz;
1325 
1326  /* this is true spontaneous rate for doubly excited state to ground
1327  * and is used for pumping, and also relaxing back to inner shell */
1328  /* Badnell gives UPWARD transition rate Alu, an unusual notation,
1329  * convert it here to the normal downward transition rate Aul */
1330  UTA.back().Emis().Aul() = (realnum)(Bij*(*UTA.back().Lo()).g()/(*UTA.back().Hi()).g());
1331  UTA.back().Emis().gf() =
1332  (realnum)GetGF( UTA.back().Emis().Aul(), UTA.back().EnergyWN(), (*UTA.back().Hi()).g() );
1333 
1334  UTA.back().Emis().iRedisFun() = ipPRD;
1335 
1336  UTA.back().Emis().dampXvel() = (realnum)((Rji+Aai)/UTA.back().EnergyWN()/PI4);
1337 
1338  ASSERT( UTA.back().Emis().dampXvel() > 0. );
1339 
1340  // remove this line if it is too weak
1341  if( UTA.back().Emis().gf() < level[ind_lo].g * f_cutoff )
1342  UTA.pop_back();
1343  }
1344  }
1345 
1346  // perform a sanity check on the tail of the file
1347  getline( ioDATA, line );
1348  ASSERT( line.substr(3,7) == "NRSLMX=" );
1349 
1350  ioDATA.close();
1351 
1352  if( trace.lgTrace )
1353  fprintf( ioQQQ, " reading %s OK\n", fnam.c_str() );
1354 }
1355 
1356 inline void InitTransition(const TransitionProxy& t)
1357 {
1358  t.AddHiState();
1359  t.AddLoState();
1360  t.AddLine2Stack();
1361 }
1362 
1363 inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl)
1364 {
1365  for( unsigned int i=0; i < level.size(); ++i )
1366  {
1367  if( level[i].irsl == irsl )
1368  return (int)i;
1369  }
1370  // we should never get here...
1371  TotalInsanity();
1372 }
1373 
1374 STATIC void validate_magic_number_1arg( const char *chFilename, FILE *ioFile, const long magicExp )
1375 {
1376  DEBUG_ENTRY( "validate_magic_number_1arg()" );
1377 
1378  char chLine[FILENAME_PATH_LENGTH_2] = { 0 };
1379  while( read_whole_line( chLine , (int)sizeof(chLine) , ioFile ) != NULL )
1380  {
1381  // skip any #
1382  if( chLine[0] != '#')
1383  break;
1384  }
1385 
1386  long magic = -1;
1387  sscanf( chLine, "%ld", &magic );
1388 
1389  if( magic != magicExp )
1390  {
1391  fprintf( ioQQQ,
1392  " atmdat_readin: the version of '%s' is not the current version.\n",
1393  chFilename );
1394  fprintf( ioQQQ,
1395  " I expected to find the number %ld and got %ld instead.\n" ,
1396  magicExp, magic );
1397  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1399  }
1400 }
1401 
1402 STATIC void validate_magic_number_3arg( const char *chFilename, FILE *ioFile,
1403  const long yearExp, const long monthExp, const long dayExp )
1404 {
1405  DEBUG_ENTRY( "validate_magic_number_3arg()" );
1406 
1407  long year,
1408  month,
1409  day;
1410 
1411  char chLine[FILENAME_PATH_LENGTH_2] = { 0 };
1412  while( read_whole_line( chLine , (int)sizeof(chLine) , ioFile ) != NULL )
1413  {
1414  // skip any #
1415  if( chLine[0] != '#')
1416  break;
1417  }
1418 
1419  sscanf( chLine, "%ld %ld %ld", &year, &month, &day );
1420 
1421  if( year != yearExp || month != monthExp || day != dayExp )
1422  {
1423  fprintf( ioQQQ,
1424  " atmdat_readin: the version of '%s' is not the current version.\n",
1425  chFilename );
1426  fprintf( ioQQQ,
1427  " I expected to find the number %ld %ld %ld and got %ld %ld %ld instead.\n" ,
1428  yearExp, monthExp, dayExp, year, month, day );
1429  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1431  }
1432 }
1433 
1435 {
1436  DEBUG_ENTRY( "read_UTA_lines()" );
1437 
1438  /* reserve space for all data sets, no worries if too small though... */
1439  UTALines.reserve( 4000 );
1440  AllTransitions.push_back(UTALines);
1441 
1442  // version of element symbols in lower case and without spaces....
1443  const char* chElmSymLC[] =
1444  { "h", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al", "si", "p",
1445  "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn", "fe", "co", "ni", "cu", "zn" };
1446 
1447  // save cite for UTA sources, insure no double counting
1448  char chUTA_ref[LIMELM][LIMELM][5];
1449  for( long nelem=0; nelem < LIMELM; ++nelem )
1450  for( long ion=0; ion <= nelem; ++ion )
1451  strcpy( chUTA_ref[nelem][ion] , "" );
1452 
1453  /* first read in the Badnell data */
1454  for( long ipISO=ipLI_LIKE; ipISO < ipAL_LIKE; ++ipISO )
1455  {
1456  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
1457  {
1458  // ion = 0 for neutral atom
1459  long ion = nelem - ipISO;
1460  strcpy( chUTA_ref[nelem][ion] , "B" );
1461 
1462  bitset<IS_TOP> Skip;
1463  if( ipISO < ipNA_LIKE )
1464  {
1465  // construct file name
1466  ostringstream oss;
1467  // Badnell calls Li-like series He-like, etc...
1468  oss << "UTA/nrb00_" << chElmSymLC[ipISO-1] << "_";
1469  // Badnell uses ion = 1 for neutral atom
1470  oss << chElmSymLC[nelem] << ion+1 << "ic1-2.dat";
1471  // now read the data...
1472  ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
1473  }
1474  else
1475  {
1476  // from Na-like onwards both K-shell (ic1-3) and L-shell (ic2-3)
1477  // excitation are treated in two separate files
1478  ostringstream oss;
1479  oss << "UTA/nrb00_" << chElmSymLC[ipISO-1] << "_";
1480  oss << chElmSymLC[nelem] << ion+1 << "ic1-3.dat";
1481  ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
1482 
1483  // Kisielius L1, L2-shell data sets take precedence for Na, Mg, and Al-like iron
1484  if( atmdat.lgInnerShell_Kisielius && nelem == ipIRON &&
1485  ipISO >= ipNA_LIKE && ipISO <= ipAL_LIKE )
1486  {
1487  Skip[IS_L1_SHELL] = 1;
1488  Skip[IS_L2_SHELL] = 1;
1489  }
1490 
1491  ostringstream oss2;
1492  oss2 << "UTA/nrb00_" << chElmSymLC[ipISO-1] << "_";
1493  oss2 << chElmSymLC[nelem] << ion+1 << "ic2-3.dat";
1494  ReadBadnellAIData( oss2.str(), nelem, ion, UTALines, Skip );
1495  }
1496  }
1497  }
1498 
1499  /* these are the statistical weights for the ground levels of all ions of iron */
1500  const realnum StatWeightGroundLevelIron[] =
1501  { 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,
1502  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 };
1503 
1504  // blank line that will be pushed on the UTA line stack
1505  qList BlankStates("BlankStates",1);
1506  TransitionList BlankList("BlankList",&BlankStates,1);
1507  TransitionList::iterator BlankLine = BlankList.begin();
1508  (*BlankLine).Junk();
1509 
1510  /* next read in the Gu file */
1511  {
1512  /* read the Gu et al. (2006) data
1513  * >>refer Fe UTA Gu, M. F., Holczer T., Behar E., & Kahn S. M. 2006, ApJ 641, 1227-1232 */
1514  if( trace.lgTrace )
1515  fprintf( ioQQQ," atmdat_readin reading UTA_Gu06.dat\n");
1516 
1517  FILE *ioGU06 = open_data( "UTA/UTA_Gu06.dat", "r" );
1518 
1519  validate_magic_number_3arg( "UTA_Gu06.dat", ioGU06, 2007, 1, 23 );
1520 
1521  int nelemGu =-1, ionGu=-1;
1522  char chLine[FILENAME_PATH_LENGTH_2] = { 0 };
1523 
1524  while( read_whole_line( chLine, (int)sizeof(chLine), ioGU06 ) != NULL )
1525  {
1526  if( chLine[0] != '#' )
1527  {
1528  long ion,
1529  i2;
1530  double EnergyAng, Aul, oscill, Aauto;
1531 
1532  sscanf( chLine, "%4li%5li%8lf%13lf%12lf",
1533  &ion, &i2, &EnergyAng, &Aul, &Aauto );
1534  sscanf( &chLine[54], "%13lf", &oscill );
1535 
1536  // avoid duplication of ions: anything upto and including the Mg-like
1537  // series is covered by the Badnell data set, Al-like is covered by
1538  // the Kisielius data set if it is turned on.
1539  int ipISO = ipIRON - ion + 1;
1540  int ipThres = atmdat.lgInnerShell_Kisielius ? ipAL_LIKE : ipMG_LIKE;
1541  if( ipISO <= ipThres )
1542  continue;
1543 
1544  UTALines.push_back( *BlankLine );
1546 
1547  /* all these are iron, first number was ion stage with 0 the atom */
1548  (*UTALines.back().Hi()).nelem() = ipIRON+1;
1549 
1550  /* now do stage of ionization */
1551  ASSERT( ion > 0 && ion <= ipIRON );
1552  /* the ion stage - 1 is atom - this data file has different
1553  * format from other two - others are 0 for atom */
1554  (*UTALines.back().Hi()).IonStg() = ion;
1555  if( ipIRON!=nelemGu || ion!=ionGu )
1556  {
1557  // one label per ion
1558  nelemGu = ipIRON;
1559  ionGu = ion;
1560  strcpy( chUTA_ref[ipIRON][ion-1] , "G" );
1561  }
1562 
1563  /* these are the statistical weights
1564  * lower levels are not included in the original data file */
1565  if( strstr_s( chLine, "(J=1/2)" ) != NULL )
1566  (*UTALines.back().Hi()).g() = 2.f;
1567  else if( strstr_s( chLine, "(J=1)" ) != NULL )
1568  (*UTALines.back().Hi()).g() = 3.f;
1569  else if( strstr_s( chLine, "(J=3/2)" ) != NULL )
1570  (*UTALines.back().Hi()).g() = 4.f;
1571  else if( strstr_s( chLine, "(J=2)" ) != NULL )
1572  (*UTALines.back().Hi()).g() = 5.f;
1573  else if( strstr_s( chLine, "(J=5/2)" ) != NULL )
1574  (*UTALines.back().Hi()).g() = 6.f;
1575  else if( strstr_s( chLine, "(J=3)" ) != NULL )
1576  (*UTALines.back().Hi()).g() = 7.f;
1577  else if( strstr_s( chLine, "(J=7/2)" ) != NULL )
1578  (*UTALines.back().Hi()).g() = 8.f;
1579  else if( strstr_s( chLine, "(J=4)" ) != NULL )
1580  (*UTALines.back().Hi()).g() = 9.f;
1581  else if( strstr_s( chLine, "(J=9/2)" ) != NULL )
1582  (*UTALines.back().Hi()).g() = 10.f;
1583  else if( strstr_s( chLine, "(J=5)" ) != NULL )
1584  (*UTALines.back().Hi()).g() = 11.f;
1585  else if( strstr_s( chLine, "(J=11/2)" ) != NULL )
1586  (*UTALines.back().Hi()).g() = 12.f;
1587  else
1588  TotalInsanity();
1589  (*UTALines.back().Lo()).g() = StatWeightGroundLevelIron[ion-1];
1590 
1591  /* wavelength in Angstroms */
1592  double fenergyWN = 1e8/EnergyAng;
1593  UTALines.back().EnergyWN() = fenergyWN;
1594  UTALines.back().WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1595 
1596  /* store branching ratio for autoionization */
1597  double frac_ioniz = Aauto/(Aul + Aauto);
1598  ASSERT( frac_ioniz >= 0. && frac_ioniz <= 1. );
1599  UTALines.back().Emis().AutoIonizFrac() = (realnum)frac_ioniz;
1600 
1601  /* save gf scanned from line */
1602  UTALines.back().Emis().gf() = (*UTALines.back().Lo()).g() * (realnum)oscill;
1603 
1604  /* this is true spontaneous rate for doubly excited state to inner
1605  * shell UTA, and is used for pumping, and also relaxing back to inner shell */
1606  UTALines.back().Emis().Aul() =
1607  (realnum)eina( UTALines.back().Emis().gf(),
1608  UTALines.back().EnergyWN(), (*UTALines.back().Hi()).g() );
1609 
1610  ASSERT( fp_equal_tol( (realnum)Aul, UTALines.back().Emis().Aul(), 1.e-3f*(realnum)Aul ) );
1611 
1612  UTALines.back().Emis().iRedisFun() = ipPRD;
1613 
1614  UTALines.back().Emis().dampXvel() = (realnum)(
1615  (UTALines.back().Emis().Aul()+Aauto) /
1616  UTALines.back().EnergyWN()/PI4);
1617  ASSERT( UTALines.back().Emis().dampXvel()>0. );
1618 
1619  // ignore line if too weak
1620  if( UTALines.back().Emis().gf() < StatWeightGroundLevelIron[ion-1] * f_cutoff )
1621  UTALines.pop_back();
1622  }
1623  }
1624 
1625  fclose( ioGU06 );
1626 
1627  if( trace.lgTrace )
1628  fprintf( ioQQQ, " reading UTA_Gu06.dat OK\n" );
1629  }
1630 
1632  {
1633  /* last read in the Romas Kisielius data
1634  *>>refer Fe UTA Kisielius, R., Hibbert, A.. Ferland, G. J., et al. 2003, MNRAS, 344, 696 */
1635  if( trace.lgTrace )
1636  fprintf( ioQQQ," atmdat_readin reading UTA_Kisielius.dat\n");
1637 
1638  FILE *ioROMAS = open_data( "UTA/UTA_Kisielius.dat", "r" );
1639 
1640  validate_magic_number_3arg( "UTA_Kisielius.dat", ioROMAS, 11, 8, 25 );
1641 
1642  long int nRomasUsed = 0 , nRomasTotal = 0;
1643  int nelemRomas=-1 , ionRomas=-1;
1644  FILE *ioROMASused;
1645  bool lgSaveRomasUsed = false;
1646  if( lgSaveRomasUsed )
1647  {
1648  if( (ioROMASused = open_data("RomasUsed.txt","w")) == NULL )
1649  {
1650  fprintf(ioQQQ,"could not open RomasUsed.txt\n");
1652  }
1653  }
1654 
1655  char chLine[FILENAME_PATH_LENGTH_2] = { 0 };
1656 
1657  while( read_whole_line( chLine, (int)sizeof(chLine), ioROMAS ) != NULL )
1658  {
1659  /* only look at lines without '#' in first col */
1660  if( chLine[0] != '#' )
1661  {
1662  long int i1, i2, i3;
1663  double f1, f2, oscill;
1664  double frac_relax;
1665 
1666  ++nRomasTotal;
1667  sscanf( chLine, "%li\t%li\t%li\t%lf\t%lf\t%lf\t%lf",
1668  &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
1669 
1670  // skip line if too weak
1671  // Fe+13 has 25 000 lines, most are vastly weak and do not add to integrated rate
1672  // the cutoff of 1e-4 reduces the number of Romas UTA lines from 84224 to 573
1673  if( oscill < f_cutoff )
1674  continue;
1675 
1676  if( lgSaveRomasUsed )
1677  fprintf(ioROMASused , "%s" , chLine);
1678 
1679  /* For Fe XIV both levels of the 2P* ground term are present in the data.
1680  * following true, ignore the transitions from the excited level
1681  * false, assume ground term populated by statistical weight if
1682  * "fudge 0" also appears in input stream */
1683  const bool lgAllowSplitFe14 = false;
1684  if( lgAllowSplitFe14 || i2 == StatWeightGroundLevelIron[i1] )
1685  {
1686  ++nRomasUsed;
1687  UTALines.push_back( *BlankLine );
1689 
1690  /* all these are iron, first number was ion stage with 0 the atom */
1691  (*UTALines.back().Hi()).nelem() = ipIRON+1;
1692 
1693  /* now do stage of ionization */
1694  ASSERT( i1 >= 0 && i1 < ipIRON );
1695  /* NB - the plus one is because 1 will be subtracted when used,
1696  * in the original data file i1 is 0 for the atom */
1697  (*UTALines.back().Hi()).IonStg() = i1 + 1;
1698 
1699  realnum facpop = 1.;
1700  // allow low or high density limit for level populations in ground term of Fe XIV
1701  if( i1 == 13 )
1702  {
1703  // Fe XIV - fudge -1 returns number of fudge parameters, 0 if not specified
1704  if( lgAllowSplitFe14 )
1705  {
1706  if( i2 == StatWeightGroundLevelIron[i1] )
1707  facpop= 0.333;
1708  else
1709  facpop = 0.667;
1710  }
1711  else
1712  {
1713  if( i2 == StatWeightGroundLevelIron[i1] )
1714  facpop= 1.;
1715  else
1716  facpop = 1e-10;
1717  }
1718  }
1719  if( ipIRON!=nelemRomas || i1!=ionRomas )
1720  {
1721  // one label per ion
1722  nelemRomas = ipIRON;
1723  ionRomas = i1;
1724  strcpy( chUTA_ref[ipIRON][i1] , "K" );
1725  }
1726 
1727  /* these were the statistical weights */
1728  (*UTALines.back().Hi()).g() = (realnum)i3;
1729  (*UTALines.back().Lo()).g() = (realnum)i2;
1730 
1731  UTALines.back().WLAng() = (realnum)f1;
1732  UTALines.back().EnergyWN() = 1.e8f/(realnum)f1;
1733  // print transitions contributing to 15.5A UTA feature
1734 # if 0
1735  if( i1==13 && f1>15.35 && f1<15.55)
1736  {
1737  fprintf(ioQQQ,"DEBUG %li\t%.5f\t%.3e\n",i2, f1 , oscill * facpop);
1738  }
1739 # endif
1740 
1741  /* this is true spontaneous rate for doubly excited state to inner shell,
1742  * and is used for pumping, and also relaxing back to inner shell */
1743  UTALines.back().Emis().gf() = (*UTALines.back().Lo()).g() * (realnum)oscill*facpop;
1744  UTALines.back().Emis().Aul() =
1745  (realnum)eina( UTALines.back().Emis().gf(),
1746  UTALines.back().EnergyWN(),
1747  (*UTALines.back().Hi()).g() );
1748 
1749  UTALines.back().Emis().iRedisFun() = ipPRD;
1750 
1751  /* store branching ratio for autoionization */
1752  ASSERT( frac_relax >= 0.f && frac_relax <= 1.f );
1753  UTALines.back().Emis().AutoIonizFrac() = 1.f-(realnum)frac_relax;
1754 
1755  // Romas data do not have autoionization rates so take typical
1756  // value of 1e15 s-1, suggested by Badnell OP data
1757  UTALines.back().Emis().dampXvel() = (realnum)(
1758  (UTALines.back().Emis().Aul()+1e15) /
1759  UTALines.back().EnergyWN()/PI4);
1760  ASSERT( UTALines.back().Emis().dampXvel()>0. );
1761  //fprintf(ioQQQ,"DEBUGGG %li %.3e\n", nRomasUsed, UTALines.back().Emis().dampXvel() );
1762  }
1763  }
1764  }
1765 
1766  fclose( ioROMAS );
1767  if( lgSaveRomasUsed )
1768  fclose( ioROMASused );
1769 
1770  if( trace.lgTrace )
1771  fprintf( ioQQQ, " reading UTA_Kisielius.dat OK,used %li lines from a total of %li\n" , nRomasUsed , nRomasTotal );
1772  }
1773 
1774  /* option to dump UTA lines, either save file, or in main output */
1776  {
1777  FILE *ioUTA = ioQQQ;
1778  if( save.lgSDSOn )
1779  ioUTA = save.ipSDSFile;
1780 
1781  fprintf(ioUTA, "##################################################\n");
1782  fprintf(ioUTA,"UTA data sources; B=Badnell 05; G==Gu 06, K=Kisielius 03, 13\n");
1783  fprintf(ioUTA," ion ");
1784  for( long ion=0; ion<=LIMELM; ++ion )
1785  fprintf(ioUTA,"%4li",ion);
1786  fprintf(ioUTA,"\n");
1787  for( long nelem=0; nelem<LIMELM; ++nelem )
1788  {
1789  fprintf(ioUTA,"%4s ", elementnames.chElementSym[nelem] );
1790  for( long ion=0; ion<=nelem; ++ion )
1791  {
1792  fprintf(ioUTA,"%4s",chUTA_ref[nelem][ion] );
1793  }
1794  fprintf(ioUTA,"\n");
1795  }
1796  fprintf(ioUTA," ion ");
1797  for( long ion=0; ion<=LIMELM; ++ion )
1798  fprintf(ioUTA,"%4li",ion);
1799  fprintf(ioUTA,"\n");
1800  fprintf(ioUTA,"Badnell 05=2005MNRAS.360..458B; Gu 06=2006ApJ...641.1227G; Kisielius 03, 13= 2003MNRAS.344..696K, 2013ApJ...767..123F\n");
1801  fprintf(ioUTA, "##################################################\n\n");
1802  }
1803 
1804  if( false )
1805  {
1806  for( size_t iu=0; iu < UTALines.size(); ++iu )
1807  {
1808  string chLab = chIonLbl( UTALines[iu] );
1809  dprintf( ioQQQ, "%5ld %s wavl %7.3f glo %2g gup %2g Aul %.2e gf %.2e ai branch %.3f\n",
1810  (long)iu,
1811  chLab.c_str(),
1812  UTALines[iu].WLAng(),
1813  (*UTALines[iu].Lo()).g(),
1814  (*UTALines[iu].Hi()).g(),
1815  UTALines[iu].Emis().Aul(),
1816  UTALines[iu].Emis().gf(),
1817  UTALines[iu].Emis().AutoIonizFrac() );
1818  }
1820  }
1821 }
realnum * hden
Definition: struc.h:25
realnum fl_energy[MEWE_FLUOR]
Definition: yield.h:35
int ion(long n) const
Definition: yield.h:60
t_mole_global mole_global
Definition: mole.cpp:7
STATIC void read_UTA_lines()
long int iter_malloc
Definition: iterations.h:40
const int ipAL_LIKE
Definition: iso.h:76
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:751
t_atmdat atmdat
Definition: atmdat.cpp:6
int nelem(long n) const
Definition: yield.h:59
realnum fl_yield[MEWE_FLUOR]
Definition: yield.h:37
int & iRedisFun() const
Definition: emission.h:438
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:230
bool lgStoutHybrid
Definition: atmdat.h:404
STATIC void read_level2_lines()
size_t size(void) const
Definition: transition.h:331
TransitionList UTALines("UTALines",&AnonStates)
realnum * windv
Definition: struc.h:25
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:296
int num_total
Definition: mole.h:362
realnum * depth_last
Definition: struc.h:25
exc_type WhichShell
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
#define NHSDIM
Definition: atmdat.h:265
int irsl2ind(vector< t_BadnellLevel > &level, int irsl)
int num_calc
Definition: mole.h:362
t_struc struc
Definition: struc.cpp:6
t_Heavy Heavy
Definition: heavy.cpp:5
void DynaCreateArrays(void)
Definition: dynamics.cpp:1398
bool lgStoutOn
Definition: atmdat.h:402
int nfl_nshell[MEWE_FLUOR]
Definition: yield.h:32
const int NISO
Definition: cddefines.h:310
double eina(double gf, double enercm, double gup)
vector< long int > nend
Definition: iterations.h:71
static const int N
void reserve(size_t newsize)
Definition: transition.h:323
void lines_setup(void)
realnum * ednstr
Definition: struc.h:25
void init_yield()
realnum ** molecules
Definition: struc.h:71
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:908
long int ncut[2][HS_NZ]
Definition: atmdat.h:337
const char * strstr_s(const char *haystack, const char *needle)
Definition: cddefines.h:1315
STATIC void init_struc()
iterator begin(void)
Definition: transition.h:339
double RefIndex(double EnergyWN)
realnum * volstr
Definition: struc.h:25
realnum g[210][4]
Definition: mewecoef.h:13
FILE * ipSDSFile
Definition: save.h:441
double * heatstr
Definition: struc.h:78
FILE * ioQQQ
Definition: cddefines.cpp:7
STATIC void validate_magic_number_1arg(const char *chFilename, FILE *ioFile, const long magicExp)
TransitionList TauLine2("TauLine2",&AnonStates)
#define MIN2(a, b)
Definition: cddefines.h:807
long int nSpecies
Definition: taulines.cpp:22
void HyperfineCreate(void)
double * coolstr
Definition: struc.h:78
t_dense dense
Definition: global.cpp:15
realnum * AccelGravity
Definition: struc.h:25
const int ipLI_LIKE
Definition: iso.h:66
static t_yield & Inst()
Definition: cddefines.h:209
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int ntemp[2][HS_NZ]
Definition: atmdat.h:335
realnum * DenMass
Definition: struc.h:25
#define NLINEHS
Definition: atmdat.h:266
t_trace trace
Definition: trace.cpp:5
#define MALLOC(exp)
Definition: cddefines.h:556
int nfl_nelem[MEWE_FLUOR]
Definition: yield.h:30
void resize(size_t newsize)
Definition: transition.h:319
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1198
bool lgChiantiHybrid
Definition: atmdat.h:376
const int ipIRON
Definition: cddefines.h:373
realnum & gf() const
Definition: emission.h:558
realnum * depth
Definition: struc.h:25
long int nDensity[2][HS_NZ]
Definition: atmdat.h:335
realnum & EnergyWN() const
Definition: transition.h:477
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
const ios_base::openmode mode_r
Definition: cpu.h:267
#define MEWE_FLUOR
Definition: yield.h:10
#define NUM_HS98_DATA_POINTS
Definition: atmdat.h:269
STATIC void validate_magic_number_3arg(const char *chFilename, FILE *ioFile, const long yearExp, const long monthExp, const long dayExp)
bool lgChiantiOn
Definition: atmdat.h:374
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
void push_back(const TransitionProxy &tr)
Definition: transition.h:347
realnum & dampXvel() const
Definition: emission.h:598
realnum * GasPressure
Definition: struc.h:25
EmissionList::reference Emis() const
Definition: transition.h:447
int nfl_nLine[MEWE_FLUOR]
Definition: yield.h:34
double **** HS_He1_Xsectn
Definition: atmdat.cpp:60
realnum * drad
Definition: struc.h:25
const int ipNA_LIKE
Definition: iso.h:74
float realnum
Definition: cddefines.h:124
realnum **** StatesElem
Definition: struc.h:67
#define EXIT_FAILURE
Definition: cddefines.h:168
double ***** OP_Helike_Xsectn
Definition: atmdat.cpp:62
long int nfl_lines
Definition: yield.h:41
qList::iterator Hi() const
Definition: transition.h:435
STATIC void read_SH98_He1_cross_sections(void)
#define cdEXIT(FAIL)
Definition: cddefines.h:484
STATIC void ReadBadnellAIData(const string &fnam, long nelem, long ion, TransitionList &UTA, bitset< IS_TOP > Skip)
double Density[2][HS_NZ][NHSDIM]
Definition: atmdat.h:329
double Emiss[2][HS_NZ][NHSDIM][NHSDIM][NLINEHS]
Definition: atmdat.h:329
bool lgSDSOn
Definition: save.h:440
STATIC void read_Hummer_Storey()
t_iterations iterations
Definition: iterations.cpp:6
long nWindLine
Definition: cdinit.cpp:19
const int ipMG_LIKE
Definition: iso.h:75
void pop_back(void)
Definition: transition.h:335
void atmdat_2phot_setSplineCoefs()
realnum * DenParticles
Definition: struc.h:25
realnum * H2_abund
Definition: struc.h:71
void atmdat_outer_shell(long int iz, long int in, long int *imax, long int *ig0, long int *ig1)
realnum * cs1_flag_lev2
Definition: taulines.cpp:37
bool lgUTAprint
Definition: atmdat.h:433
const TransitionProxy back(void)
Definition: transition.h:351
bool lgElmtOn[LIMELM]
Definition: dense.h:160
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
long int nzlim
Definition: struc.h:19
int nfl_ion_emit[MEWE_FLUOR]
Definition: yield.h:33
#define HS_NZ
Definition: atmdat.h:267
#define ASSERT(exp)
Definition: cddefines.h:617
realnum * drad_x_fillfac
Definition: struc.h:25
qList::iterator Lo() const
Definition: transition.h:431
void AddHiState() const
Definition: transition.cpp:651
const int ipCALCIUM
Definition: cddefines.h:367
bool lgInnerShell_Kisielius
Definition: atmdat.h:431
realnum * AccelTotalOutward
Definition: struc.h:25
double GetGF(double trans_prob, double enercm, double gup)
const int LIMELM
Definition: cddefines.h:307
long int n_elec_eject[30][30][7]
Definition: yield.h:26
double **** HS_He1_Energy
Definition: atmdat.cpp:61
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
const int ipHELIUM
Definition: cddefines.h:349
double ***** OP_Helike_Energy
Definition: atmdat.cpp:63
realnum * drad_last
Definition: struc.h:25
realnum ** gas_phase
Definition: struc.h:75
long **** OP_Helike_NumPts
Definition: atmdat.cpp:64
realnum frac_elec_eject[30][30][7][10]
Definition: yield.h:25
realnum * histr
Definition: struc.h:25
bool lgInnerShellLine_on
Definition: atmdat.h:429
void database_readin(void)
Definition: species.cpp:42
realnum * xLyman_depth
Definition: struc.h:25
#define MAX2(a, b)
Definition: cddefines.h:828
int nfl_ion[MEWE_FLUOR]
Definition: yield.h:31
void atmdat_readin(void)
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
STATIC void read_mewe_gbar()
const int ipPRD
Definition: cddefines.h:339
bool lgKillAuger
Definition: yield.h:44
#define S(I_, J_)
vector< TransitionList > AllTransitions
Definition: taulines.cpp:9
const realnum f_cutoff
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:72
long nint(double x)
Definition: cddefines.h:762
realnum & AutoIonizFrac() const
Definition: emission.h:628
realnum * testr
Definition: struc.h:25
exc_type
#define fixit(a)
Definition: cddefines.h:416
void InitTransition(const TransitionProxy &t)
bool lgLamdaOn
Definition: atmdat.h:391
t_save save
Definition: save.cpp:5
t_MeweCoef MeweCoef
Definition: mewecoef.cpp:5
const int ipHYDROGEN
Definition: cddefines.h:348
STATIC void read_Helike_cross_sections(void)
const int ipLITHIUM
Definition: cddefines.h:350
realnum * hiistr
Definition: struc.h:25
realnum & Aul() const
Definition: emission.h:668
void AddLoState() const
Definition: transition.cpp:640
realnum *** xIonDense
Definition: struc.h:64
realnum & WLAng() const
Definition: transition.h:468
realnum * o3str
Definition: struc.h:25
EmissionList & Emis()
Definition: transition.h:363
double ElecTemp[2][HS_NZ][NHSDIM]
Definition: atmdat.h:329
void AddLine2Stack() const
Definition: transition.cpp:628
realnum * pressure
Definition: struc.h:25
#define EXIT_SUCCESS
Definition: cddefines.h:166
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:393
#define UNUSED
Definition: cpu.h:14
realnum * pres_radiation_lines_curr
Definition: struc.h:25