cloudy  trunk
 All Data Structures 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-2008 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 "yield.h"
13 #include "trace.h"
14 #include "lines.h"
15 #include "lines_service.h"
16 #include "ionbal.h"
17 #include "struc.h"
18 #include "geometry.h"
19 #include "dynamics.h"
20 #include "elementnames.h"
21 #include "hyperfine.h"
22 #include "atmdat.h"
23 /* */
24 /* this was needed to get array to crash out of bounds if not set.
25  * std limits on limits.h did not work with visual studio! */
26 #define INTBIG 2000000000
27 
28 /* these are the individual pointers to the level 1 lines, they are set to
29  * very large negative.
30  * NB NB NB!!
31  * these occur two times in the code!!
32  * They are declared in taulines.h, and defined here */
33 long
34  ipT1656=INTBIG, ipT9830=INTBIG, /*ipH21cm=INTBIG, ipHe3cm=INTBIG,*/
99  ipFe106375=INTBIG,ipT353=INTBIG, /*ipFe1310=INTBIG, ipFe1311=INTBIG,*/
101  ipT191=INTBIG, /*ipTFe07=INTBIG, ipTFe61=INTBIG,*/ ipFe18975=INTBIG, ipTFe23=INTBIG,
113 
114 /* above are the individual pointers to the level 1 lines, they are set to
115  * very large negative.
116  * NB NB NB!!
117  * these occur two times in the code!!
118  * They are declared in taulines.h, and defined here */
119 
120 /* definition for whether level 2 lines are enabled, will be set to -1
121  * with no level2 command */
122 /*long nWindLine = NWINDDIM;*/
123 /*realnum TauLine2[NWINDDIM][NTA];*/
124 /*realnum **TauLine2;*/
125 #include "atomfeii.h"
126 
127 void atmdat_readin(void)
128 {
129  long int i,
130  iCase ,
131  ion,
132  ipDens ,
133  ipISO ,
134  ipTemp ,
135  j,
136  J,
137  ig0,
138  ig1,
139  imax,
140  nelem,
141  nelec,
142  magic1,
143  magic2,
144  mol ,
145 
146  nUTA_Romas,
147  nFeIonRomas[ipIRON],
148 
149  nUTA_Behar ,
150  nFeIonBehar[ipIRON],
151 
152  nUTA_Gu06,
153  nFeIonGu[ipIRON];
154  double
155  WLloRomas[ipIRON],
156  WLhiRomas[ipIRON],
157  WLloBehar[ipIRON],
158  WLhiBehar[ipIRON],
159  WLloGu[ipIRON],
160  WLhiGu[ipIRON];
161 
162  char cha;
163  char chS2[3];
164 
165  long ipZ;
166  bool lgEOL;
167 
168  FILE *ioDATA, *ioROMAS , *ioBEHAR=NULL;
169  FILE *ioGU06=NULL;
170  char chLine[FILENAME_PATH_LENGTH_2] ,
171  chFilename[FILENAME_PATH_LENGTH_2];
172 
173  static bool lgFirstCall = true;
174 
175  DEBUG_ENTRY( "atmdat_readin()" );
176 
177  /* do nothing if not first call */
178  if( !lgFirstCall )
179  {
180  /* do not do anything, but make sure that number of zones has not increased */
181  bool lgTooBig = false;
182  for( j=0; j < iterations.iter_malloc; j++ )
183  {
184  if( geometry.nend[j]>=struc.nzlim )
185  lgTooBig = true;
186  }
187  if( lgTooBig )
188  {
189  fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
190  fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
191  fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
192  fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
193  fprintf(ioQQQ," Sorry.\n");
194  cdEXIT(EXIT_FAILURE);
195  }
196  return;
197  }
198 
199  lgFirstCall = false; /* do not reevaluate again */
200 
201  /* make sure that molecules have been initialized - this will fail
202  * if this routine is called before size of molecular network is known */
203  if( !mole.num_comole_calc )
204  {
205  /* mole.num_comole_calc can't be zero */
206  TotalInsanity();
207  }
208 
209  /* create space for the structure variables */
210  /* nzlim will be limit, and number allocated */
211  /* >>chng 01 jul 28, define this var, do all following mallocs */
212  for( j=0; j < iterations.iter_malloc; j++ )
213  {
214  struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
215  }
216 
217  /* sloppy, but add one extra for safely */
218  ++struc.nzlim;
219 
220  struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
221 
222  struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
223 
224  struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
225 
226  struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
227 
228  struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
229 
230  struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
231 
232  struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
233 
234  struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
235 
236  struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
237 
238  struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
239 
240  struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
241 
242  struc.AccelTot = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
243 
244  struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
245 
246  struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
247 
248  struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
249 
250  struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
251 
252  struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
253 
254  struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
255 
256  struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
257 
258  struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
259 
260  struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
261 
262  struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
263 
264  struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
265 
266  struc.H2_molec = ((realnum**)MALLOC( (size_t)(N_H_MOLEC)*sizeof(realnum* )));
267 
268  struc.CO_molec = ((realnum**)MALLOC( (size_t)(mole.num_comole_calc)*sizeof(realnum* )));
269 
270  for(mol=0;mol<N_H_MOLEC;mol++)
271  {
272  struc.H2_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
273  }
274 
275  for(mol=0;mol<mole.num_comole_calc;mol++)
276  {
277  struc.CO_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
278  }
279 
280  /* create space for gas phase abundances array, first create space for the elements */
281  struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)LIMELM );
282  /* last the zones */
283  for( ipZ=0; ipZ< LIMELM;++ipZ )
284  {
285  struc.gas_phase[ipZ] =
286  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
287  }
288 
289  /* create space for struc.xIonDense array, first create space for the zones */
290  struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM+3) );
291 
292  for( ipZ=0; ipZ< (LIMELM+3);++ipZ )
293  {
294  /* space for the ionization stage */
295  struc.xIonDense[ipZ] =
296  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+1) );
297  /* now create diagonal of space for zones */
298  for( ion=0; ion < (LIMELM+1); ++ion )
299  {
300  struc.xIonDense[ipZ][ion] =
301  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
302  }
303  }
304 
305  /* some structure variables */
306  for( i=0; i < struc.nzlim; i++ )
307  {
308  struc.testr[i] = 0.;
309  struc.volstr[i] = 0.;
310  struc.drad_x_fillfac[i] = 0.;
311  struc.histr[i] = 0.;
312  struc.hiistr[i] = 0.;
313  struc.ednstr[i] = 0.;
314  struc.o3str[i] = 0.;
315  struc.heatstr[i] = 0.;
316  struc.coolstr[i] = 0.;
317  struc.pressure[i] = 0.;
319  struc.GasPressure[i] = 0.;
320  struc.DenParticles[i] = 0.;
321  struc.depth[i] = 0.;
322  for(mol=0;mol<N_H_MOLEC;mol++)
323  {
324  struc.H2_molec[mol][i] = 0.;
325  }
326  for(mol=0;mol<mole.num_comole_calc;mol++)
327  {
328  struc.CO_molec[mol][i] = 0.;
329  }
330  }
331 
332  /* allocate space for some arrays used by dynamics routines, and zero out vars */
333  DynaCreateArrays( );
334 
335  /*************************************************************
336  * *
337  * get the level 1 lines, with real collision data set *
338  * *
339  *************************************************************/
340 
341  if( trace.lgTrace )
342  fprintf( ioQQQ," atmdat_readin opening level1.dat:");
343 
344  ioDATA = open_data( "level1.dat", "r" );
345 
346  /* first line is a version number and does not count */
347  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
348  {
349  fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
350  cdEXIT(EXIT_FAILURE);
351  }
352 
353  /* count how many lines are in the file, ignoring all lines
354  * starting with '#' */
355  nLevel1 = 0;
356  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
357  {
358  /* we want to count the lines that do not start with #
359  * since these contain data */
360  if( chLine[0] != '#')
361  ++nLevel1;
362  }
363 
364  /* now malloc the TauLines array */
365  TauLines = (transition *)MALLOC( (size_t)(nLevel1+1)*sizeof(transition ) );
366 
367  for( i=0; i<nLevel1+1; i++ )
368  {
369  TransitionJunk( &TauLines[i] );
370  TauLines[i].Hi = AddState2Stack();
371  TauLines[i].Lo = AddState2Stack();
372  TauLines[i].Emis = AddLine2Stack( true );
373  }
374 
375  /* now rewind the file so we can read it a second time*/
376  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
377  {
378  fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
379  cdEXIT(EXIT_FAILURE);
380  }
381 
382  /* check that magic number is ok */
383  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
384  {
385  fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
386  cdEXIT(EXIT_FAILURE);
387  }
388  i = 1;
389  /* level 1 magic number */
390  nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
391  nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
392  ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
393 
394  /* magic number
395  * the following is the set of numbers that appear at the start of level1.dat */
396  if( ( nelem != 5 ) || ( nelec != 12 ) || ( ion != 19 ) )
397  {
398  fprintf( ioQQQ,
399  " atmdat_readin: the version of level1.dat is not the current version.\n" );
400  fprintf( ioQQQ,
401  " Please obtain the current version from the Cloudy web site.\n" );
402  fprintf( ioQQQ,
403  " I expected to find the number 04 11 5 and got %li %li %li instead.\n" ,
404  nelem , nelec , ion );
405  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
406  cdEXIT(EXIT_FAILURE);
407  }
408 
409  /* this starts at 1 not 0 since zero is reserved for the
410  * dummy line - we counted number of lines above */
411  i = 1;
412  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
413  {
414  if( i >= (nLevel1+1) )
415  {
416  fprintf(ioQQQ," version conflict.\n");
417  }
418  /* only look at lines without '#' in first col */
419  if( chLine[0] != '#')
420  {
421  /* >>chng 00 apr 24, prevents crash in PresTotCurrent on call from ContSetIntensity, by PvH */
422 
423 
424  /*sscanf( chLine ,
425  "%2s%2li%9li %11f %2li %2li %10g %10g %2li\n",
426  chS2,&i2,&i3,&f1,&i4,&i5,&f2,&f3,&i6 );*/
427 
428  /* first two cols of line has chem element symbol,
429  * try to match it with the list of names
430  * there are no H lines in the level 1 set so zero
431  * is an invalid result */
432 
433  /* copy first two char into chS2 and null terminate */
434  strncpy( chS2 , chLine , 2);
435  chS2[2] = 0;
436 
437  ipZ = 0;
438  for( j=0; j<LIMELM; ++j)
439  {
440  if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
441  {
442  /* ipZ is the actual atomic number starting from 1 */
443  ipZ = j + 1;
444  break;
445  }
446  }
447 
448  /* this happens if no valid chemical element in first two cols on this line */
449  if( ipZ == 0 )
450  {
451  fprintf( ioQQQ,
452  " atmdat_readin could not identify chem symbol on this level 1line:\n");
453  fprintf( ioQQQ, "%s\n", chLine );
454  fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
455  cdEXIT(EXIT_FAILURE);
456  }
457 
458  /* now stuff them into the array, will convert over to proper units later */
459  TauLines[i].Hi->nelem = (int)ipZ;
460  /* start the scan early on the line -- the element label will not be
461  * picked up, first number will be the ion stage after the label, as in c 4 */
462  j = 1;
463  TauLines[i].Hi->IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
464  if( lgEOL )
465  {
466  fprintf( ioQQQ, " There should have been a number on this level1 line 1. Sorry.\n" );
467  fprintf( ioQQQ, "string==%s==\n" ,chLine );
468  cdEXIT(EXIT_FAILURE);
469  }
470 
471  TauLines[i].Lo->nelem = TauLines[i].Hi->nelem;
472  TauLines[i].Lo->IonStg = TauLines[i].Hi->IonStg;
473 
474  TauLines[i].WLAng = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
475  if( lgEOL )
476  {
477  fprintf( ioQQQ, " There should have been a number on this level1 line 2. Sorry.\n" );
478  fprintf( ioQQQ, "string==%s==\n" ,chLine );
479  cdEXIT(EXIT_FAILURE);
480  }
481  /*TauLines[i].WLAng = (realnum)i3;*/
482  TauLines[i].EnergyWN = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
483  if( lgEOL )
484  {
485  fprintf( ioQQQ, " There should have been a number on this level1 line 3. Sorry.\n" );
486  fprintf( ioQQQ, "string==%s==\n" ,chLine );
487  cdEXIT(EXIT_FAILURE);
488  }
489  /*TauLines[i].EnergyWN = (realnum)f1;*/
490  TauLines[i].Lo->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
491  if( lgEOL )
492  {
493  fprintf( ioQQQ, " There should have been a number on this level1 line 4. Sorry.\n" );
494  fprintf( ioQQQ, "string==%s==\n" ,chLine );
495  cdEXIT(EXIT_FAILURE);
496  }
497  /*TauLines[i].Lo->g = (realnum)i4;*/
498  TauLines[i].Hi->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
499  if( lgEOL )
500  {
501  fprintf( ioQQQ, " There should have been a number on this level1 line 5. Sorry.\n" );
502  fprintf( ioQQQ, "string==%s==\n" ,chLine );
503  cdEXIT(EXIT_FAILURE);
504  }
505  /*TauLines[i].Hi->g = (realnum)i5;*/
506 
507  /* gf is log if negative */
508  TauLines[i].Emis->gf = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
509  if( lgEOL )
510  {
511  fprintf( ioQQQ, " There should have been a number on this level1 line 6. Sorry.\n" );
512  fprintf( ioQQQ, "string==%s==\n" ,chLine );
513  cdEXIT(EXIT_FAILURE);
514  }
515  /*TauLines[i].Emis->gf = (realnum)f2;*/
516  if( TauLines[i].Emis->gf < 0. )
517  TauLines[i].Emis->gf = (realnum)pow((realnum)10.f,TauLines[i].Emis->gf);
518 
519  /* Emis.Aul is log if negative */
520  TauLines[i].Emis->Aul = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
521  if( lgEOL )
522  {
523  fprintf( ioQQQ, " There should have been a number on this level1 line 7. Sorry.\n" );
524  fprintf( ioQQQ, "string==%s==\n" ,chLine );
525  cdEXIT(EXIT_FAILURE);
526  }
527  if( TauLines[i].Emis->Aul < 0. )
528  TauLines[i].Emis->Aul = (realnum)pow((realnum)10.f,TauLines[i].Emis->Aul);
529  /*TauLines[i].Emis->Aul = f3;*/
530 
531  TauLines[i].Emis->iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
532  if( lgEOL )
533  {
534  fprintf( ioQQQ, " There should have been a number on this level1 line 8. Sorry.\n" );
535  fprintf( ioQQQ, "string==%s==\n" ,chLine );
536  cdEXIT(EXIT_FAILURE);
537  }
538  /*TauLines[i].Emis->iRedisFun = i6;*/
539 
540  /* this is new in c94.01 and returns nothing (0.) for most lines */
541  TauLines[i].Coll.cs = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
542 
543  /* finally increment i */
544  ++i;
545  }
546  }
547 
548  /* now close the file */
549  fclose(ioDATA);
550 
551  /*************************************************************
552  * *
553  * get the level 2 line, opacity project, data set *
554  * *
555  *************************************************************/
556 
557  /* nWindLine is initialized to the dimension of the vector when it is
558  * initialized in the definition at the start of this file.
559  * it is set to -1 with the "no level2" command, which
560  * stops us from trying to establish this vector */
561  if( nWindLine > 0 )
562  {
563  /* begin level2 level 2 wind block */
564  /* open file with level 2 line data */
565 
566  /* create the TauLine2 emline array */
567  TauLine2 = ((transition *)MALLOC( (size_t)nWindLine*sizeof(transition )));
568  cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum )));
569 
570  /* first initialize entire array to dangerously large negative numbers */
571  for( i=0; i< nWindLine; ++i )
572  {
573  /* >>chng 99 jul 16, from setting all t[] to flt_max, to call for
574  * following, each member of structure set to own type of impossible value */
575  TransitionJunk( &TauLine2[i] );
576 
577  TauLine2[i].Hi = AddState2Stack();
578  TauLine2[i].Lo = AddState2Stack();
579  TauLine2[i].Emis = AddLine2Stack( true );
580  }
581 
582  if( trace.lgTrace )
583  fprintf( ioQQQ," atmdat_readin opening level2.dat:");
584 
585  ioDATA = open_data( "level2.dat", "r" );
586 
587  /* get magic number off first line */
588  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
589  {
590  fprintf( ioQQQ, " level2.dat error getting magic number\n" );
591  cdEXIT(EXIT_FAILURE);
592  }
593 
594  /* scan off three numbers, should be the yr mn dy of creation date */
595  i = 1;
596  nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
597  nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
598  ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
599 
600  /* level2.dat magic number
601  * the following is the set of numbers that appear at the start of level2.dat 01 05 29
602  * >>chng 06 aug 10, chng to 06, 08, 10 */
603  if( lgEOL || ( nelem != 6 ) || ( nelec != 8 ) || ( ion != 10 ) )
604  {
605  fprintf( ioQQQ,
606  " atmdat_readin: the version of level2.dat is not the current version.\n" );
607  fprintf( ioQQQ,
608  " I expected to find the number 06 08 10 and got %li %li %li instead.\n" ,
609  nelem , nelec , ion );
610  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
611  fprintf( ioQQQ, "Please obtain the correct version.\n");
612  cdEXIT(EXIT_FAILURE);
613  }
614 
615  /* now get the actual data */
616  i = 0;
617  while( i < nWindLine )
618  {
619  /* this must be double for sscanf to work below */
620  double tt[7];
621  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
622  {
623  fprintf( ioQQQ, " level2.dat error getting line %li\n", i );
624  cdEXIT(EXIT_FAILURE);
625  }
626  /* option to skip any line starting with # */
627  if( chLine[0]!='#' )
628  {
629  /*printf("string is %s\n",chLine );*/
630  sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " ,
631  &tt[0] ,
632  &tt[1] ,
633  &tt[2] ,
634  &tt[3] ,
635  &tt[4] ,
636  &tt[5] ,
637  &tt[6] );
638  /* these are readjusted into their final form in the structure
639  * in routine lines_setup*/
640  TauLine2[i].Hi->nelem = (int)tt[0];
641  TauLine2[i].Hi->IonStg = (int)tt[1];
642  TauLine2[i].Lo->g = (realnum)tt[2];
643  TauLine2[i].Hi->g = (realnum)tt[3];
644  TauLine2[i].Emis->gf = (realnum)tt[4];
645  TauLine2[i].EnergyWN = (realnum)tt[5];
646  cs1_flag_lev2[i] = (realnum)tt[6];
647  ++i;
648  }
649  }
650  /* get magic number off last line */
651  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
652  {
653  fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
654  cdEXIT(EXIT_FAILURE);
655  }
656  sscanf( chLine , "%ld" , &magic2 );
657  if( 999 != magic2 )
658  {
659  fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n",
660  magic2 );
661  cdEXIT(EXIT_FAILURE);
662  }
663  fclose( ioDATA );
664  if( trace.lgTrace )
665  fprintf( ioQQQ," OK \n");
666 
667  /* end of level2 block*/
668  }
669 
670  /*************************
671  *
672  * get the UTA line sets
673  * >>chng 06 nov 26, use Gu et al. 2006 UTA data by default
674  * with option to use older Behar
675  *
676  *************************/
677 
678  /* these will count number of UTA lines of each type */
679  nUTA = 0;
680  nUTA_Gu06 = 0;
681  nUTA_Behar = 0;
682  nUTA_Romas = 0;
683  /* this is option to use older Behar data rather than new Gu data
684  * variable lgInnerShell_Gu06 is true by default, set false with
685  * SET UTA BEHAR command */
687  {
688  if( trace.lgTrace )
689  fprintf( ioQQQ," atmdat_readin opening UTA_Gu06.dat:");
690 
691  ioGU06 = open_data( "UTA_Gu06.dat", "r" );
692 
693  /* count how many non-comment lines there are in Gu06 file */
694  while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
695  {
696  /* we want to count the lines that do not start with #
697  * since these contain data */
698  if( chLine[0] != '#')
699  ++nUTA_Gu06;
700  }
701  /* we counted a magic number as a real line, back off by one */
702  ASSERT( nUTA_Gu06 > 0 );
703  --nUTA_Gu06;
704 
705  /* now rewind the file so we can read it a second time*/
706  if( fseek( ioGU06 , 0 , SEEK_SET ) != 0 )
707  {
708  fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Gu06.dat.\n");
709  cdEXIT(EXIT_FAILURE);
710  }
711 
712  /* get up to magic number in Gu 06 file - there is a large header of
713  * comments with the first data the magic number */
714  /* count how many non-comment lines there are in Gu06 file */
715  while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
716  {
717  /* we want to break on first line that does not start with #
718  * since that contains data */
719  if( chLine[0] != '#')
720  break;
721  }
722  /* now get Gu magic number */
723  j = 1;
724  nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
725  nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
726  ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
727 
728  /* is magic number correct?
729  * the following is the set of numbers that appear at the start of UTA_Gu06.dat 2006 11 17 */
730  if( ( nelem != 2007 ) || ( nelec != 1 ) || ( ion != 23 ) )
731  {
732  fprintf( ioQQQ,
733  " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
734  fprintf( ioQQQ,
735  " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
736  nelem , nelec , ion );
737  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
738  cdEXIT(EXIT_FAILURE);
739  }
740  }
741  else
742  {
743  /* this branch, get the Behar 01 data */
744  /* >>refer Fe inner shell UTA Behar, E., Sako, M, Kahn, S.M.,
745  * >>refercon 2001, ApJ, 563, 497-504
746  * the fits were based on the paper
747  * >>refer Fe inner shell UTA Behar, E., & Netzer, H., 2002, ApJ, 570, 165-170 */
748 
749  if( trace.lgTrace )
750  fprintf( ioQQQ," atmdat_readin opening UTA_Behar.dat:");
751 
752  ioBEHAR = open_data( "UTA_Behar.dat", "r" );
753 
754  /* count number of lines in Behar file
755  * first line is a version number and does not count */
756  if( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
757  {
758  fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n");
759  cdEXIT(EXIT_FAILURE);
760  }
761 
762  j = 1;
763  /* UTA_Behar.dat magic number */
764  nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
765  nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
766  ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
767 
768  /* magic number
769  * the following is the set of numbers that appear at the start of UTA_Behar.dat 2002 8 19 */
770  /* >>chng 02 oct 22, had incorrect stat weights for all stages of ionization except atom.
771  * no longer used stored values when stat weight is zero */
772  if( ( nelem != 2002 ) || ( nelec != 10 ) || ( ion != 22 ) )
773  {
774  fprintf( ioQQQ,
775  " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" );
776  fprintf( ioQQQ,
777  " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
778  nelem , nelec , ion );
779  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
780  cdEXIT(EXIT_FAILURE);
781  }
782  }
783 
784  if( trace.lgTrace )
785  fprintf( ioQQQ," atmdat_readin UTA data files open OK\n");
786 
787  /* count how many lines are in the file, ignoring all lines
788  * starting with '#' */
789  nUTA_Behar = 0;
790 
791  /* first fill in Behar & Netzer interpolation - first pass just count number of lines */
792  for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
793  {
794  for( nelem=ipISO; nelem<LIMELM; ++nelem )
795  {
796  ++nUTA_Behar;
797  }
798  }
799 
801  {
802  /* count number of lines in Behar file if we are going to use it
803  * wee always use above inner shell UTA from B&N */
804  while( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) != NULL )
805  {
806  /* we want to count the lines that do not start with #
807  * since these contain data */
808  if( chLine[0] != '#')
809  ++nUTA_Behar;
810  }
811  /* now rewind the file so we can read it a second time*/
812  if( fseek( ioBEHAR , 0 , SEEK_SET ) != 0 )
813  {
814  fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Behar.dat.\n");
815  cdEXIT(EXIT_FAILURE);
816  }
817  /* first line is a version number and does not count */
818  if( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
819  {
820  fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n");
821  cdEXIT(EXIT_FAILURE);
822  }
823  }
824 
825  /* now read Romas Kisielius data, taken from
826  * >>refer Fe UTA Kisielius, R., Hibbert, A., Ferland, G.J., Foord, M.E., Rose, S.J.,
827  * >>refercon van Hoof, P.A.M. & Keenan, F.P. 2003, MNRAS, 344, 696 */
828 
829  if( trace.lgTrace )
830  fprintf( ioQQQ," atmdat_readin opening UTA_Kisielius.dat:");
831 
832  ioROMAS = open_data( "UTA_Kisielius.dat", "r" );
833 
834  nUTA_Romas = 0;
835  while( read_whole_line( chLine , (int)sizeof(chLine) , ioROMAS ) != NULL )
836  {
837  /* we want to count the lines that do not start with #
838  * since these contain data */
839  if( chLine[0] != '#')
840  {
841  ++nUTA_Romas;
842  }
843  }
844 
845  /* now malloc the UTA lines array with enough space for both
846  * data sets */
847  UTALines = (transition *)MALLOC( (size_t)(nUTA_Behar+nUTA_Romas+nUTA_Gu06)*sizeof(transition) );
848 
849  for( i=0; i<nUTA_Behar+nUTA_Romas+nUTA_Gu06; i++ )
850  {
851  TransitionJunk( &UTALines[i] );
852  UTALines[i].Hi = AddState2Stack();
853  UTALines[i].Lo = AddState2Stack();
854  UTALines[i].Emis = AddLine2Stack( true );
855  }
856 
857  /* these will be used to keep track of what is in each data set */
858  for( ion=0; ion<ipIRON; ++ion )
859  {
860  nFeIonRomas[ion] = 0;
861  WLloRomas[ion] = 1e10;
862  WLhiRomas[ion] = 0.;
863 
864  nFeIonBehar[ion] = 0;
865  WLloBehar[ion] = 1e10;
866  WLhiBehar[ion] = 0.;
867 
868  nFeIonGu[ion] = 0;
869  WLloGu[ion] = 1e10;
870  WLhiGu[ion] = 0.;
871  }
872 
873  /* this will count number of lines, must equal sum of UTAs at end */
874  i = 0;
875  /* now rewind the file so we can read it a second time*/
876  if( fseek( ioROMAS , 0 , SEEK_SET ) != 0 )
877  {
878  fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Kisielius.dat.\n");
879  cdEXIT(EXIT_FAILURE);
880  }
881 
882  i = 0;
883  /* first fill in Behar & Netzer interpolation - now do real thing,
884  * these are Nick Abel's fits to the BN data */
885  for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
886  {
887  for( nelem=ipISO; nelem<LIMELM; ++nelem )
888  {
889  double f1;
890  double sum_spon_auto;
891 # define N 7
892 
893  double alam[N]={ 1.859035 , 4.692138 , 10.324543 , 19.294364 ,
894  40.401292 , 44.442754 , 72.959719 };
895  double blam[N]={-9.7855968,-11.159739, -12.914931 , -14.987272,
896  -18.853446 , -19.271781 ,-23.828388 };
897  double clam[N]={ 8.2874628, 8.3043002, 8.3164038 , 8.3269937,
898  8.4312895 , 8.3730893 , 8.493802 };
899  /* >>chng 05 mar 07, by NA, update fits to better form */
900 
901  double aA[N] = {1.9067 , 1.8715 , 2.1033, 2.2815 ,
902  1.9511 , 1.9966 , 2.0852};
903  double bA[N] = {8.4538 , 8.5528, 7.616, 6.9247 ,
904  7.9479, 7.9777 , 7.8349 };
905 
906  double aDep[N] = {29.54 , 12.07 , 24.23 , 7.35 , 7.37 , 11.14 , 11.99 };
907  double bDep[N] = {-14.853 , -7.606 , -7.844 , 9.987 , 10.503 , 8.541 , 8.865 };
908 
909  /* NB - these are plus one is because 1 will be subtracted when used */
910  /* derive element and ion */
911  UTALines[i].Hi->nelem = nelem+1;
912  ion = nelem + 1 - ipISO;
913  UTALines[i].Hi->IonStg = ion;
914 
915  /* these were the statistical weights */
917  UTALines[i].Hi->g = 4.f;
918  UTALines[i].Lo->g = 2.f;
919 
920  f1 = alam[ipISO-2]*1e-4 + blam[ipISO-2]*1e-4*(nelem+1) +
921  clam[ipISO-2]*1e-4*POW2(nelem+1.);
922  f1 = (1./f1);
923  UTALines[i].WLAng = (realnum)f1;
924  UTALines[i].EnergyWN = (realnum)(1e8/f1);
926 
927  f1 = aA[ipISO-2]*log((double)(nelem+1)) + bA[ipISO-2];
928 
929  UTALines[i].Emis->Aul = (realnum)pow(10., f1 );
930  /* generate Emis.gf from A */
931  UTALines[i].Emis->gf =
932  (realnum)(UTALines[i].Emis->Aul*UTALines[i].Hi->g*
933  1.4992e-8*POW2(1e4/UTALines[i].EnergyWN));
934 
935  UTALines[i].Emis->iRedisFun = 1;
936 
937  /* find sum of spontaneous relax and autoionization As */
938  f1 = log((double)(nelem+1));
939  f1 = POW2( f1 );
940  sum_spon_auto = aDep[ipISO-2]*1e-2*f1 + bDep[ipISO-2]*1e-1;
941 
942  if( ipISO == ipBE_LIKE )
943  {
944  sum_spon_auto = exp( sum_spon_auto );
945  sum_spon_auto = pow(10., sum_spon_auto )*1e13;
946  }
947  else if( ipISO <= ipF_LIKE )
948  {
949  sum_spon_auto = pow(10., sum_spon_auto )*1e13;
950  }
951  else
952  TotalInsanity();
953 
954  /* last number is negative branching ratio for autoionization,
955  * which we will store as negative of cs so not misinterpreted as cs */
956  /* >>chng 03 feb 18, change to negative of br */
957 
958  /*ASSERT( UTALines[i].Emis->Aul / sum_spon_auto <= 1. && UTALines[i].Aul / sum_spon_auto >= 0. );*/
959  UTALines[i].Coll.cs = (realnum)MIN2(0., UTALines[i].Emis->Aul / sum_spon_auto - 1. );
960 
961  /* option to print UTAs */
962 # if 0
963  fprintf(ioQQQ," %2s%2li\t%.4f\t%.3e\t%.3e\t%.3e\n",
964  elementnames.chElementSym[nelem],
965  ion,
966  UTALines[i].WLAng ,
967  UTALines[i].Emis->Aul ,
968  sum_spon_auto ,
969  -UTALines[i].cs);
970 # endif
971 
972  /* finally increment i */
973  ++i;
974 # undef N
975  }
976  }
977 
978  /* next read in the Gu file */
980  {
981  ioDATA = ioGU06;
982  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
983  {
984  /* only look at lines without '#' in first col */
985  if( chLine[0] != '#')
986  {
987  long int i2;
988  double WLAng, Aul, oscill , Aauto;
989 
990  sscanf( chLine ,"%4li%5li%8lf%13le%12le",
991  &ion,&i2,&WLAng,&Aul,&Aauto);
992  sscanf( &chLine[54] ,"%13le", &oscill);
993 
994  /* all these are iron, first number was ion stage with 0 the atom */
995  UTALines[i].Hi->nelem = ipIRON+1;
996 
997  /* now do stage of ionization */
998  ASSERT( ion >= 0 && ion < ipIRON );
999  /* the ion stage - 1 is atom - this data file has different
1000  * format from other two - others are 0 for atom */
1001  UTALines[i].Hi->IonStg = ion;
1002 
1003  /* these were the statistical weights
1004  * not included in the original data files, and not needed,
1005  * so use unity */
1006  UTALines[i].Hi->g = 1.f;
1007  UTALines[i].Lo->g = 1.f;
1008 
1009  /* wavelength in Angstroms */
1010  UTALines[i].WLAng = (realnum)WLAng;
1011 
1012  /* keep track of what we have */
1013  ++nFeIonGu[ion-1];
1014  WLloGu[ion-1] = MIN2( WLAng , WLloGu[ion-1] );
1015  WLhiGu[ion-1] = MAX2( WLAng , WLloGu[ion-1] );
1016 
1017  /* energy in wavenumbers */
1018  UTALines[i].EnergyWN = (realnum)(1e8/WLAng);
1020 
1021  /* this is true spontaneous rate for doubly excited state to inner
1022  * shell UTA, and is used for pumping, and also relaxing back to inner shell */
1023  UTALines[i].Emis->Aul = (realnum)Aul;
1024 
1025  /* cs is branching ratio for autoionization,
1026  * which we will store as negative cs so not misinterpreted as real cs
1027  * this multiplies the direct pump rate to get ionization rate
1028  * rate evaluated in conv_base.cpp*/
1029  double frac_ioniz = Aauto/(Aul + Aauto);
1030  ASSERT( frac_ioniz >=0. && frac_ioniz <=1. );
1031  UTALines[i].Coll.cs = -(realnum)frac_ioniz;
1032 
1033  /* save gf scanned from line */
1034  UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill;
1035 
1036  /* statistical weights are not given in the data file -
1037  * assume unity and convert absorption oscillator strength
1038  * into an effective Aul, which is correct for unit statistical
1039  * weight */
1040  UTALines[i].Emis->Aul = (realnum)eina(UTALines[i].Emis->gf,
1041  UTALines[i].EnergyWN,UTALines[i].Hi->g);
1042 
1043  UTALines[i].Emis->iRedisFun = 1;
1044  {
1045  /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */
1046  /*@-redef@*/
1047  enum {DEBUG_LOC=false};
1048  /*@+redef@*/
1049  if( DEBUG_LOC && UTALines[i].Hi->nelem==ipIRON+1 && UTALines[i].Hi->IonStg==9)
1050  {
1051  fprintf(ioQQQ,"DEBUG Gu UTA\t%li\t%.3f\t%.3e\t%.3e\t%.3e\n",
1052  ion , WLAng,Aul ,
1053  eina(UTALines[i].Emis->gf,
1054  UTALines[i].EnergyWN,UTALines[i].Hi->g),
1055  Aauto );
1056  }
1057  }
1058 
1059  /* finally increment i */
1060  ++i;
1061  }
1062  }
1063  }
1064  else
1065  {
1066  /* ionbal.lgInnerShell_Gu06 is false, do not use this data */
1067  nUTA_Gu06 = 0;
1068 
1069  /* the Behar 01 data */
1070  /* read in the Behar data */
1071  ioDATA = ioBEHAR;
1072  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1073  {
1074  /* only look at lines without '#' in first col */
1075  if( chLine[0] != '#')
1076  {
1077  long int i1, i2, i3;
1078  double f1, f2, oscill;
1079  double frac_relax;
1080 
1081  sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
1082  &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
1083 
1084  /* all these are iron, first number was ion stage with 0 the atom */
1085  UTALines[i].Hi->nelem = ipIRON+1;
1086 
1087  /* now do stage of ionization */
1088  ASSERT( i1 >= 0 && i1 < ipIRON );
1089  /* NB - the plus one is because 1 will be subtracted when used,
1090  * in the original data file i1 is 0 for the atom */
1091  UTALines[i].Hi->IonStg = i1 + 1;
1092 
1093  /* these were the statistical weights */
1094  if( i2>0 && i3>0 )
1095  {
1096  UTALines[i].Hi->g = (realnum)i2;
1097  UTALines[i].Lo->g = (realnum)i3;
1098  }
1099  else
1100  {
1101  /* nearly all stages of ionization do not include the stat weights,
1102  * so use unity */
1103  UTALines[i].Hi->g = 1.f;
1104  UTALines[i].Lo->g = 1.f;
1105  }
1106 
1107  /* wavelength in Angstroms */
1108  UTALines[i].WLAng = (realnum)f1;
1109 
1110  /* keep track of what we have */
1111  ++nFeIonBehar[i1];
1112  WLloBehar[i1] = MIN2( f1 , WLloBehar[i1] );
1113  WLhiBehar[i1] = MAX2( f1 , WLhiBehar[i1] );
1114 
1115  /* energy in wavenumbers */
1116  UTALines[i].EnergyWN = (realnum)(1e8/f1);
1118 
1119  /* this is true spontaneous rate for doubly excited state to inner shell UTA,
1120  * and is used for pumping, and also relaxing back to inner shell */
1121  UTALines[i].Emis->Aul = (realnum)f2;
1122 
1123  /* >>chng 02 oct 22, use absorption oscillator strength in data file
1124  * since don't have stat weights for most lines */
1125  UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill;
1126 
1127  UTALines[i].Emis->iRedisFun = 1;
1128 
1129  /* last number is branching ratio for autoionizing,
1130  * which we will store in negative in cs so not misinterpreted as real cs*/
1131  /* >>chng 03 feb 18, cs to negative of branching ratio */
1132  ASSERT( frac_relax >=0. && frac_relax <=1. );
1133  UTALines[i].Coll.cs = (realnum)(frac_relax-1.);
1134 # if 0
1135  fprintf(ioQQQ,"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n",
1136  nDataSet,
1137  i,
1139  UTALines[i].Hi->IonStg,
1140  UTALines[i].WLAng ,
1141  UTALines[i].Emis->Aul ,
1142  -UTALines[i].cs);
1143 # endif
1144 
1145  /* finally increment the counter i */
1146  ++i;
1147  }
1148  }
1149  }
1150 
1151  /* last read in the Romas Kisielius data last since will use it in preference to other
1152  * two sets
1153  *>>refer Fe UTA Kisielius, R.; Hibbert, A.; Ferland, G. J.; Foord, M. E.; Rose, S. J.;
1154  *>>refercon van Hoof, P. A. M.; Keenan, F. P. 2003, MNRAS, 344, 696
1155  */
1156  ioDATA = ioROMAS;
1157  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
1158  {
1159  /* only look at lines without '#' in first col */
1160  if( chLine[0] != '#')
1161  {
1162  long int i1, i2, i3;
1163  double f1, f2, oscill;
1164  double frac_relax;
1165 
1166  sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
1167  &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
1168 
1169  /* all these are iron, first number was ion stage with 0 the atom */
1170  UTALines[i].Hi->nelem = ipIRON+1;
1171 
1172  /* now do stage of ionization */
1173  ASSERT( i1 >= 0 && i1 < ipIRON );
1174  /* NB - the plus one is because 1 will be subtracted when used,
1175  * in the original data file i1 is 0 for the atom */
1176  UTALines[i].Hi->IonStg = i1 + 1;
1177 
1178  /* these were the statistical weights */
1179  if( i2>0 && i3>0 )
1180  {
1181  UTALines[i].Hi->g = (realnum)i2;
1182  UTALines[i].Lo->g = (realnum)i3;
1183  }
1184  else
1185  {
1186  /* nearly all stages of ionization do not include the stat weights,
1187  * so use unity */
1188  UTALines[i].Hi->g = 1.f;
1189  UTALines[i].Lo->g = 1.f;
1190  }
1191 
1192  /*TauLines[i].Hi->IonStg = i2;*/
1193  UTALines[i].WLAng = (realnum)f1;
1194 
1195  /* keep track of what we have */
1196  ++nFeIonRomas[i1];
1197  WLloRomas[i1] = MIN2( f1 , WLloRomas[i1] );
1198  WLhiRomas[i1] = MAX2( f1 , WLhiRomas[i1] );
1199 
1200  UTALines[i].EnergyWN = (realnum)(1e8/f1);
1202 
1203  /* this is true spontaneous rate for doubly excited state to inner shell,
1204  * and is used for pumping, and also relaxing back to inner shell */
1205  UTALines[i].Emis->Aul = (realnum)f2;
1206  /* generate gf from A */
1207  /* >>chng 02 oct 22, use absorption oscillator strength in data file
1208  * since don't have stat weights for most lines */
1209  UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill;
1210 
1211  UTALines[i].Emis->iRedisFun = 1;
1212 
1213  /* last number is branching ratio for autoionizing,
1214  * which we will store in negative in cs so not misinterpreted as real cs*/
1215  /* >>chng 03 feb 18, cs to negative of br */
1216  ASSERT( frac_relax >=0. && frac_relax <=1. );
1217  UTALines[i].Coll.cs = (realnum)(frac_relax-1.);
1218 
1219  /* finally increment i */
1220  ++i;
1221  }
1222  }
1223 
1224  /* these have to agree */
1225  ASSERT( i == (nUTA_Behar+nUTA_Romas+nUTA_Gu06) );
1226 
1227  if( trace.lgTrace )
1228  {
1229  fprintf(ioQQQ , "summary of UTA data sets read in\nion numb WLlo WLhi\n");
1230  fprintf(ioQQQ , "Behar 01 total lines %li\n" , nUTA_Behar );
1231  for( ion=0; ion<ipIRON;++ion )
1232  {
1233  if( nFeIonBehar[ion] )
1234  {
1235  fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n",
1236  ion,
1237  nFeIonBehar[ion],
1238  WLloBehar[ion],
1239  WLhiBehar[ion]);
1240  }
1241  }
1242 
1243  fprintf(ioQQQ , "Gu 06 total lines %li\n" , nUTA_Gu06 );
1244  for( ion=0; ion<ipIRON;++ion )
1245  {
1246  if( nFeIonGu[ion] )
1247  {
1248  fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n",
1249  ion,
1250  nFeIonGu[ion],
1251  WLloGu[ion],
1252  WLhiGu[ion]);
1253  }
1254  }
1255 
1256  fprintf(ioQQQ , "Romas Kisielius 03 total lines %li\n", nUTA_Romas );
1257  for( ion=0; ion<ipIRON;++ion )
1258  {
1259  if( nFeIonRomas[ion] )
1260  {
1261  fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n",
1262  ion,
1263  nFeIonRomas[ion],
1264  WLloRomas[ion],
1265  WLhiRomas[ion] );
1266  }
1267  }
1268  }
1269 
1270  /* this is default option to use the data Romas Kisielius computed
1271  * to replace any existing data . */
1272 
1273  /* here there are two options - we can add the Romas Kisielius data to the
1274  * first block of UTA data, but zero out any overlapping data,
1275  * or simply ignore the Romas Kisielius data */
1277  {
1278  for( i=0; i<nUTA_Behar+nUTA_Gu06; ++i )
1279  {
1280  if( UTALines[i].Hi->nelem-1 == ipIRON )
1281  {
1282  /* nFeIonRomas counts how many of the Romas Kisielius lines
1283  * fall in this stage of ionization - if zero then
1284  * none, so keep Behar data, if positive, then
1285  * there are new data, and zero out Behar data */
1286  ASSERT( UTALines[i].Hi->IonStg-1< ipIRON );
1287  if( nFeIonRomas[UTALines[i].Hi->IonStg-1] )
1288  {
1289  UTALines[i].Emis->Aul = -1.;
1290  }
1291  }
1292  }
1293  /* the total number of UTAs, a global variable */
1294  nUTA = nUTA_Behar + nUTA_Gu06 + nUTA_Romas;
1295  }
1296  else
1297  {
1298  /* the Kisielius data come last, and nUTA is the global variable
1299  * that counts the number of UTA lines. Setting nUTA to only the
1300  * first two amounts to ignoring the Kisielius data */
1301  nUTA = nUTA_Behar + nUTA_Gu06;
1302  }
1303 
1304  /* now close the files */
1305  fclose(ioROMAS);
1307  fclose(ioGU06);
1308  else
1309  fclose(ioBEHAR);
1310 
1311  /* end UTA */
1312 
1313  /* allocate space for the carbon monoxide CO rotation levels */
1314  ASSERT( nCORotate > 0);
1315  C12O16Rotate = (transition *)MALLOC( (size_t)nCORotate*sizeof(transition) );
1316  C13O16Rotate = (transition *)MALLOC( (size_t)nCORotate*sizeof(transition) );
1317 
1318  /* this says we cannot change the number of levels again, since space is allocated */
1319  lgCORotateMalloc = true;
1320 
1321  /* next initialize entire array to dangerously large negative numbers */
1322  for( J=0; J< nCORotate; ++J )
1323  {
1327  C12O16Rotate[J].Emis = AddLine2Stack( true );
1328 
1332  C13O16Rotate[J].Emis = AddLine2Stack( true );
1333  }
1334 
1335  /* read in data for the set of hyperfine structure lines, and allocate
1336  * space for the transition HFLines[nHFLines] structure */
1337  HyperfineCreate();
1338 
1339  /* this is Humeshkar's work on generic atoms and molecules. */
1340  //Nemala_Start();
1341 
1342  /* initialize the large block of level 1 real lines, and OP level 2 lines */
1343  lines_setup();
1344 
1345  /* mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat ========*/
1346  /* read in g-bar data taken from
1347  *>>refer all cs-gbar Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985,
1348  *>>refercon A&AS, 62, 197 */
1349  /* open file with Mewe coefficients */
1350 
1351  if( trace.lgTrace )
1352  fprintf( ioQQQ," atmdat_readin opening mewe_gbar.dat:");
1353 
1354  ioDATA = open_data( "mewe_gbar.dat", "r" );
1355 
1356  /* get magic number off first line */
1357  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1358  {
1359  fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" );
1360  cdEXIT(EXIT_FAILURE);
1361  }
1362  /* check the number */
1363  sscanf( chLine , "%ld" , &magic1 );
1364  if( magic1 != 9101 )
1365  {
1366  fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n",
1367  magic1 );
1368  cdEXIT(EXIT_FAILURE);
1369  }
1370 
1371  /* now get the actual data, indices are correct for c, in Fort went from 2 to 210 */
1372  for( i=1; i < 210; i++ )
1373  {
1374  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1375  {
1376  fprintf( ioQQQ, " mewe_gbar.dat error getting line %li\n", i );
1377  cdEXIT(EXIT_FAILURE);
1378  }
1379  /*printf("%s\n",chLine);*/
1380  double help[4];
1381  sscanf( chLine, "%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] );
1382  for( int l=0; l < 4; ++l )
1383  MeweCoef.g[i][l] = (realnum)help[l];
1384  }
1385 
1386  /* get magic number off last line */
1387  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1388  {
1389  fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
1390  cdEXIT(EXIT_FAILURE);
1391  }
1392 
1393  sscanf( chLine , "%ld" , &magic2 );
1394 
1395  if( magic1 != magic2 )
1396  {
1397  fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n",
1398  magic2 );
1399  cdEXIT(EXIT_FAILURE);
1400  }
1401 
1402  fclose( ioDATA );
1403 
1404  if( trace.lgTrace )
1405  fprintf( ioQQQ," OK \n");
1406 
1407  /* This is what remains of the t_yield initialization
1408  * this should not be in the constructor of t_yield
1409  * since it initializes a different struct */
1410  for( nelem=0; nelem < LIMELM; nelem++ )
1411  for( ion=0; ion < LIMELM; ion++ )
1412  Heavy.nsShells[nelem][ion] = LONG_MAX;
1413 
1414  /* now read in all auger yields
1415  * will do elements from li on up,
1416  * skip any line starting with #
1417  * this loop goes from lithium to Zn */
1418  for( nelem=2; nelem < LIMELM; nelem++ )
1419  {
1420  /* nelem is on the shifted C scale, so 2 is Li */
1421  for( ion=0; ion <= nelem; ion++ )
1422  {
1423  /* number of bound electrons, = atomic number for neutral */
1424  nelec = nelem - ion + 1;
1425  /* one of dima's routines to determine the number of electrons
1426  * for this species, nelem +1 to shift to physical number */
1427  /* subroutine atmdat_outer_shell(iz,in,imax,ig0,ig1)
1428  * iz - atomic number from 1 to 30 (integer)
1429  * in - number of electrons from 1 to iz (integer)
1430  * Output: imax - number of the outer shell
1431  */
1432  atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);
1433 
1434  ASSERT( imax > 0 && imax <= 10 );
1435 
1436  /* nsShells[nelem][ion] is outer shell number for ion with nelec electrons
1437  * on physics scale, with K shell being 1 */
1438  Heavy.nsShells[nelem][ion] = imax;
1439  }
1440  }
1441 
1442  /*************************************************************
1443  * *
1444  * get the Auger electron yield data set *
1445  * *
1446  *************************************************************/
1447 
1449 
1450  /****************************************************************
1451  * *
1452  * get the Hummer and Storey model case A, B results, these are *
1453  * the two data files e1b.dat and e2b.dat, for H and He *
1454  * *
1455  ****************************************************************/
1456 
1457  /* now get Hummer and Storey case b data, loop over H, He */
1458  /* >>chng 01 aug 08, generalized this to both case A and B, and all elements
1459  * up to HS_NZ, now 8 */
1460  for( ipZ=0; ipZ<HS_NZ; ++ipZ )
1461  {
1462  /* don't do the minor elements, Li-B */
1463  if( ipZ>1 && ipZ<5 ) continue;
1464 
1465  for( iCase=0; iCase<2; ++iCase )
1466  {
1467  /* open file with case B data */
1468  /* >>chng 01 aug 08, add HS_ to start of file names to indicate Hummer Storey origin */
1469  /* create file name for this charge
1470  * first character after e is charge number for this element,
1471  * then follows whether this is case A or case B data */
1472  sprintf( chFilename, "HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ? 'a' : 'b' );
1473 
1474  if( trace.lgTrace )
1475  fprintf( ioQQQ," atmdat_readin opening Hummer Storey emission file %s:",chFilename);
1476 
1477  /* open the file */
1478  ioDATA = open_data( chFilename, "r" );
1479 
1480  if( trace.lgTrace )
1481  {
1482  fprintf( ioQQQ," OK\n");
1483  }
1484 
1485  /* read in the number of temperatures and densities*/
1486  i = fscanf( ioDATA, "%li %li ",
1487  &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );
1488 
1489  /* check that ntemp and nDensity are below NHSDIM,
1490  * set to 15 in atmdat_HS_caseB.h */
1491  assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
1492  assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM);
1493 
1494  /* loop reading in line emissivities for all temperatures*/
1495  for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
1496  {
1497  for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
1498  {
1499  long int junk, junk2 , ne;
1500  fscanf( ioDATA, " %lf %li %lf %c %li %ld ",
1501  &atmdat.Density[iCase][ipZ][ipDens], &junk ,
1502  &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 ,
1503  &atmdat.ncut[iCase][ipZ] );
1504 
1505  ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
1506  ASSERT( ne<=NLINEHS );
1507  for( j=0; j < ne; j++ )
1508  {
1509  fscanf( ioDATA, "%lf ", &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
1510  }
1511  }
1512  }
1513 
1514  /*this is end of read-in loop */
1515  fclose(ioDATA);
1516 # if 0
1517  /* print list of densities and temperatures */
1518  for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
1519  {
1520  fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]);
1521  }
1522  fprintf(ioQQQ,"\n");
1523  for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
1524  {
1525  fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]);
1526  }
1527  fprintf(ioQQQ,"\n");
1528 # endif
1529  }
1530  }
1531 
1532  /* Verner's model FeII atom
1533  * do not call if Netzer model used, or it Fe(2) is zero
1534  * exception is when code is searching for first soln */
1535  /* read the atomic data from files */
1536  /* convert line arrays to internal form needed by code */
1537  FeIICreate();
1538  return;
1539 }
1540 
1542 {
1543  /* preset two arrays that will hold auger data
1544  * set to very large values so that code will blow
1545  * if not set properly below */
1546  for( int nelem=0; nelem < LIMELM; nelem++ )
1547  {
1548  for( int ion=0; ion < LIMELM; ion++ )
1549  {
1550  for( int ns=0; ns < 7; ns++ )
1551  {
1552  n_elec_eject[nelem][ion][ns] = LONG_MAX;
1553  for( int nelec=0; nelec < 10; nelec++ )
1554  {
1555  frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
1556  }
1557  }
1558  }
1559  }
1560 
1561  lgKillAuger = false;
1562 }
1563 
1565 {
1566  char chLine[FILENAME_PATH_LENGTH_2];
1567  const char* chFilename;
1568 
1569  /* following is double for sscanf to work */
1570  double temp[15];
1571 
1572  DEBUG_ENTRY( "init_yield()" );
1573 
1574  /*************************************************************
1575  * *
1576  * get the Auger electron yield data set *
1577  * *
1578  *************************************************************/
1579 
1580  /* NB NB -- This test of Heavy.nsShells remains here to assure
1581  * that it contains meaningful values since needed below, once
1582  * t_Heavy is a Singleton, it can be removed !!! */
1583  ASSERT( Heavy.nsShells[2][0] > 0 );
1584 
1585  /* hydrogen and helium will not be done below, so set yields here*/
1586  n_elec_eject[ipHYDROGEN][0][0] = 1;
1587  n_elec_eject[ipHELIUM][0][0] = 1;
1588  n_elec_eject[ipHELIUM][1][0] = 1;
1589 
1590  frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
1591  frac_elec_eject[ipHELIUM][0][0][0] = 1;
1592  frac_elec_eject[ipHELIUM][1][0][0] = 1;
1593 
1594  /* open file auger.dat, yield data file that came from
1595  * >>refer all auger Kaastra, J.S., & Mewe, R. 1993, A&AS, 97, 443-482 */
1596  chFilename = "mewe_nelectron.dat";
1597 
1598  if( trace.lgTrace )
1599  fprintf( ioQQQ, " init_yield opening %s:", chFilename );
1600 
1601  FILE *ioDATA;
1602 
1603  /* open the file */
1604  ioDATA = open_data( chFilename, "r" );
1605 
1606  /* now read in all auger yields
1607  * will do elements from li on up,
1608  * skip any line starting with #
1609  * this loop goes from lithium to Zn */
1610  for( int nelem=2; nelem < LIMELM; nelem++ )
1611  {
1612  /* nelem is on the shifted C scale, so 2 is Li */
1613  for( int ion=0; ion <= nelem; ion++ )
1614  {
1615  for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1616  {
1617  char ch1 = '#';
1618  /* the * is the old comment char, accept it, but really want # */
1619  while( ch1 == '#' || ch1 == '*' )
1620  {
1621  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1622  {
1623  fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns );
1624  cdEXIT(EXIT_FAILURE);
1625  }
1626  ch1 = chLine[0];
1627  }
1628  /*printf("string is %s\n",chLine );*/
1629  sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
1630  &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1631  &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
1632  &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] );
1633  n_elec_eject[nelem][ion][ns] = (long int)temp[3];
1634 
1635  ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 );
1636 
1637  /* can pop off up to 10 auger electrons, these are the probabilities*/
1638  for( int j=0; j < 10; j++ )
1639  {
1640  frac_elec_eject[nelem][ion][ns][j] = (realnum)temp[j+4];
1641  ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. );
1642  }
1643  }
1644  }
1645  /* activate this print statement to get yield for k-shell of all atoms */
1646  /*fprintf(ioQQQ,"yyyield\t%li\t%.4e\n", nelem+1, frac_elec_eject[nelem][0][0][0] );*/
1647  }
1648 
1649  fclose( ioDATA );
1650 
1651  if( trace.lgTrace )
1652  {
1653  /* this is set with no auger command */
1654  if( lgKillAuger )
1655  fprintf( ioQQQ, " Auger yields will be killed.\n");
1656  fprintf( ioQQQ, " OK\n");
1657  }
1658 
1659  /* open file mewe_fluor.dat, yield data file that came from
1660  * >>refer all auger Kaastra, J.S., & Mewe, R. 1993, 97, 443-482 */
1661  chFilename = "mewe_fluor.dat";
1662 
1663  if( trace.lgTrace )
1664  fprintf( ioQQQ, " init_yield opening %s:", chFilename );
1665 
1666  /* open the file */
1667  ioDATA = open_data( chFilename, "r" );
1668 
1669  /* now read in all auger yields
1670  * will do elements from li on up,
1671  * skip any line starting with # */
1672  do
1673  {
1674  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1675  {
1676  fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 );
1677  cdEXIT(EXIT_FAILURE);
1678  }
1679  }
1680  while( chLine[0] == '#' );
1681 
1682  bool lgEOL = false;
1683 
1684  nfl_lines = 0;
1685  do
1686  {
1687  const int NKM = 10;
1688  int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
1689  int nAuger;
1690 
1691  if( nfl_lines >= MEWE_FLUOR )
1692  TotalInsanity();
1693 
1694  /*printf("string is %s\n",chLine );*/
1695  sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf",
1696  &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
1697  &temp[5], &temp[6] );
1698 
1699  /* the atomic number, C is 5 */
1700  nfl_nelem[nfl_lines] = (int)temp[0]-1;
1701  ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM );
1702 
1703  /* the ion stage for target, atom is 0 */
1704  nfl_ion[nfl_lines] = (int)temp[1]-1;
1706 
1707  /* the target's shell */
1708  nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1];
1709  ASSERT( nfl_nshell[nfl_lines] >= 0 &&
1710  /* nsShells is shell number, where K is 1 */
1712 
1713  /* this is the number of Auger electrons ejected */
1714  nAuger = (int)temp[3];
1715  /* so this is the spectrum of the photons that are emitted */
1716  nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1;
1717  /* must be gt 0 since at least photoelectron comes off */
1719 
1720  /* this is the type of line as defined in their paper */
1721  nfl_nLine[nfl_lines] = (int)temp[4];
1722 
1723  /* energy in Ryd */
1724  fl_energy[nfl_lines] = (realnum)temp[5] / (realnum)EVRYD;
1725  ASSERT( fl_energy[nfl_lines] > 0. );
1726 
1727  /* fluor yield */
1728  fl_yield[nfl_lines] = (realnum)temp[6];
1729  /* NB cannot assert <=1 since data file has yields around 1.3 - 1.4 */
1730  ASSERT( fl_yield[nfl_lines] >= 0 );
1731 
1732  ++nfl_lines;
1733 
1734  do
1735  {
1736  if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
1737  lgEOL = true;
1738  }
1739  while( chLine[0]=='#' && !lgEOL );
1740  }
1741  while( !lgEOL );
1742 
1743  fclose( ioDATA );
1744 }

Generated for cloudy by doxygen 1.8.1.1