cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_chianti.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 #include "cddefines.h"
4 #include "taulines.h"
5 #include "trace.h"
6 #include "thirdparty.h"
7 #include "atmdat.h"
8 #include "lines_service.h"
9 #include "parse_species.h"
10 
11 typedef vector< pair<double,long> > DoubleLongPairVector;
12 
20 static void AulTTError( const char chFilename[], const char chLine[], const char TT[] )
21 {
22  DEBUG_ENTRY( "AulTTError()" );
23 
24  fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chFilename);
25  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
26  fprintf( ioQQQ, " This transition already has an Aul value set for %s.\n",TT);
27  cdEXIT(EXIT_FAILURE);
28 
29 }
30 
32 template<class T>
33 long getSortedLevel( const T& input, long oldLev)
34 {
35  DEBUG_ENTRY( "getSortedLevel()" );
36 
37  long returnLev = -1;
38  /* Account for reordered energy levels */
39  typename T::const_iterator i;
40  for( i = input.begin(); i != input.end(); i++ )
41  {
42  if( oldLev == *i )
43  {
44  returnLev = i-input.begin();
45  break;
46  }
47  }
48 
49  return returnLev;
50 }
51 template<>
52 long getSortedLevel( const DoubleLongPairVector& dBaseStatesEnergy, long oldLev)
53 {
54  DEBUG_ENTRY( "getSortedLevel()" );
55 
56  long returnLev = -1;
57  /* Account for reordered energy levels */
58  for( DoubleLongPairVector::const_iterator i = dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
59  {
60  if( oldLev == i->second )
61  {
62  returnLev = i-dBaseStatesEnergy.begin();
63  break;
64  }
65  }
66 
67  return returnLev;
68 }
69 
70 static const bool DEBUGSTATE = false;
71 // minimum energy difference (wavenumbers) we will accept
72 const double ENERGY_MIN_WN = 1e-10;
73 
74 void atmdat_STOUT_readin( long intNS, char *chPrefix )
75 {
76  DEBUG_ENTRY( "atmdat_STOUT_readin()" );
77 
78  long int nMolLevs;
79  char chLine[FILENAME_PATH_LENGTH_2],
80  chNRGFilename[FILENAME_PATH_LENGTH_2],
81  chTPFilename[FILENAME_PATH_LENGTH_2],
82  chCOLLFilename[FILENAME_PATH_LENGTH_2];
83 
84  static const int MAX_NUM_LEVELS = 999;
85 
86  dBaseSpecies[intNS].lgMolecular = false;
87  dBaseSpecies[intNS].lgLTE = false;
88 
89  strcpy( chNRGFilename , chPrefix );
90  strcpy( chTPFilename , chPrefix );
91  strcpy( chCOLLFilename , chPrefix );
92 
93  /*Open the energy levels file*/
94  strcat( chNRGFilename , ".nrg");
95  uncaps( chNRGFilename );
96 
97  /*Open the files*/
98  if( trace.lgTrace )
99  fprintf( ioQQQ," atmdat_STOUT_readin opening %s:",chNRGFilename);
100 
101  FILE *ioDATA;
102  ioDATA = open_data( chNRGFilename, "r" );
103  bool lgEOL = false;
104  long index = 0;
105  double nrg = 0.0;
106  double stwt = 0.0;
107 
108  /* first line is a version number - now confirm that it is valid */
109  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
110  {
111  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chNRGFilename );
112  cdEXIT(EXIT_FAILURE);
113  }
114  long int ipFFmt = 1;
115  const long int iyr = 11, imo=10 , idy = 14;
116  long iyrread, imoread , idyread;
117  iyrread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
118  imoread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
119  idyread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
120 
121  if(( iyrread != iyr ) ||
122  ( imoread != imo ) ||
123  ( idyread != idy ) )
124  {
125  fprintf( ioQQQ,
126  " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chNRGFilename );
127  fprintf( ioQQQ,
128  " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
129  iyr, imo , idy ,iyrread, imoread , idyread );
130  cdEXIT(EXIT_FAILURE);
131  }
132 
133  nMolLevs = 0;
134  //Count number of levels
135  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
136  {
137  /* # - comment, *** ends data */
138  if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
139  {
140  ipFFmt = 1;
141  long n = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
142  if( n < 0 )
143  break;
144  nMolLevs++;
145  }
146  else if( (chLine[0] == '*' && chLine[1] == '*' ) )
147  {
148  /* stop reading when field of stars encountered.*/
149  break;
150  }
151  }
152 
153  /* now rewind the file so we can read it a second time*/
154  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
155  {
156  fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chNRGFilename );
157  cdEXIT(EXIT_FAILURE);
158  }
159  //Skip the magic numbers this time
160  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
161 
162  long HighestIndexInFile = nMolLevs;
163 
164  dBaseSpecies[intNS].numLevels_max = HighestIndexInFile;
165 
166  setProperties(dBaseSpecies[intNS]);
167 
168  if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
169  {
170  // Fe is special case with more levels
171  nMolLevs = MIN3(nMolLevs, atmdat.nStoutMaxLevelsFe, MAX_NUM_LEVELS );
172  }
173  else if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'o' && dBaseSpecies[intNS].chLabel[3] == '1')
174  {
175  // O 1 needs more levels for H Lyman Beta
176  nMolLevs = MIN3(nMolLevs, atmdat.nDefaultPhotoLevelsFe, MAX_NUM_LEVELS );
177  }
178  else
179  {
180  nMolLevs = MIN3(nMolLevs, atmdat.nStoutMaxLevels, MAX_NUM_LEVELS );
181  }
182 
183  //Consider the masterlist specified number of levels as the min. =1 if not specified
184  long numMasterlist = MIN2( dBaseSpecies[intNS].numLevels_masterlist , HighestIndexInFile );
185  nMolLevs = MAX2(nMolLevs,numMasterlist);
186 
187  if (dBaseSpecies[intNS].setLevels != -1)
188  {
189  if (dBaseSpecies[intNS].setLevels > HighestIndexInFile)
190  {
191  char chLabelChemical[CHARS_SPECIES] = "";
192  spectral_to_chemical( chLabelChemical, dBaseSpecies[intNS].chLabel ),
193  fprintf( ioQQQ,"Using STOUT spectrum %s (species: %s) with %li requested, only %li energy levels available.\n",
194  dBaseSpecies[intNS].chLabel, chLabelChemical, dBaseSpecies[intNS].setLevels, HighestIndexInFile );
195  nMolLevs = HighestIndexInFile;
196  }
197  else
198  {
199  nMolLevs = dBaseSpecies[intNS].setLevels;
200  }
201  }
202 
203  if( atmdat.lgStoutPrint == true)
204  {
205  char chLabelChemical[CHARS_SPECIES] = "";
206  spectral_to_chemical( chLabelChemical, dBaseSpecies[intNS].chLabel ),
207  fprintf( ioQQQ,"Using STOUT spectrum %s (species: %s) with %li levels of %li available.\n",
208  dBaseSpecies[intNS].chLabel, chLabelChemical, nMolLevs, HighestIndexInFile );
209  }
210 
211  dBaseSpecies[intNS].numLevels_max = nMolLevs;
212  dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
213 
214  /*Resize the States array*/
215  dBaseStates[intNS].init(dBaseSpecies[intNS].chLabel,nMolLevs);
216  /*Zero out the maximum wavenumber for each species */
217  dBaseSpecies[intNS].maxWN = 0.;
218 
219  /* allocate the Transition array*/
220  ipdBaseTrans[intNS].reserve(nMolLevs);
221  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
222  ipdBaseTrans[intNS].reserve(ipHi,ipHi);
223  ipdBaseTrans[intNS].alloc();
224  dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
225  dBaseTrans[intNS].states() = &dBaseStates[intNS];
226  dBaseTrans[intNS].chLabel() = dBaseSpecies[intNS].chLabel;
227  dBaseSpecies[intNS].database = "Stout";
228 
229  //This is creating transitions that we don't have collision data for
230  //Maybe use gbar or keep all of the Fe 2 even if it was assumed (1e-5)
231  int ndBase = 0;
232  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
233  {
234  for( long ipLo = 0; ipLo < ipHi; ipLo++)
235  {
236  ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
237  dBaseTrans[intNS][ndBase].Junk();
238  dBaseTrans[intNS][ndBase].setLo(ipLo);
239  dBaseTrans[intNS][ndBase].setHi(ipHi);
240  ++ndBase;
241  }
242  }
243 
244  /* Create arrays for holding energies and statistical weights so that
245  * we can put the energies in the correct order before moving them to
246  * dBaseStates */
247  DoubleLongPairVector dBaseStatesEnergy;
248  vector<double> dBaseStatesStwt(HighestIndexInFile,-1.0);
249  for( long ii = 0; ii < HighestIndexInFile; ii++ )
250  {
251  dBaseStatesEnergy.push_back(make_pair(-1.0,ii));
252  }
253 
254  if( DEBUGSTATE )
255  {
256  fprintf(ioQQQ,"\nStout Species: %s\n",dBaseSpecies[intNS].chLabel);
257  fprintf(ioQQQ,"Energy Level File: %s\n",chNRGFilename);
258  fprintf(ioQQQ,"Number of Energy Levels in File: %li\n",HighestIndexInFile);
259  fprintf(ioQQQ,"Number of Energy Levels Cloudy is Currently Using: %li\n",nMolLevs);
260  fprintf(ioQQQ,"Species|File Index|Energy(WN)|StatWT\n");
261  }
262 
263  /* Check for an end of file sentinel */
264  bool lgSentinelReached = false;
265 
266  //Read the first line of data
267  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
268  {
269  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the energy level file.\n");
270  cdEXIT(EXIT_FAILURE);
271  }
272  //Read the remaining lines of the energy level file
273  do
274  {
275  ipFFmt = 1;
276 
277  /* Stop on *** */
278  if( chLine[0] == '*' )
279  {
280  lgSentinelReached = true;
281  break;
282  }
283  //Comments start with #, skip them
284  if( chLine[0] != '#' )
285  {
286  /* Skip blank lines */
287  if( chLine[0] == '\n')
288  continue;
289 
290  index = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL) -1 ;
291  if( DEBUGSTATE )
292  {
293  fprintf(ioQQQ,"<%s>\t%li\t",dBaseSpecies[intNS].chLabel,index+1);
294  }
295 
296  if( index < 0 )
297  {
298  fprintf( ioQQQ, " PROBLEM An energy level index was less than 1 in file %s\n",chNRGFilename);
299  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
300  cdEXIT(EXIT_FAILURE);
301  }
302 
303  nrg = (double)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
304  stwt = (double)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
305  if( DEBUGSTATE )
306  {
307  fprintf(ioQQQ,"%.3f\t%.1f\n",nrg,stwt);
308  }
309 
310  if( lgEOL )
311  {
312  fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chNRGFilename);
313  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
314  cdEXIT(EXIT_FAILURE);
315  }
316 
317  /* Verify that energy levels are not overwritten */
318  if( dBaseStatesEnergy.at(index).first == -1. && dBaseStatesStwt.at(index) == -1 )
319  {
320  dBaseStatesEnergy.at(index) = make_pair(nrg,index);
321  dBaseStatesStwt.at(index) = stwt;
322  }
323  else
324  {
325  fprintf( ioQQQ, " PROBLEM Duplicate energy level index in file %s\n",chNRGFilename);
326  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
327  cdEXIT(EXIT_FAILURE);
328  }
329  }
330  }
331  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
332  if( !lgSentinelReached )
333  {
334  fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chNRGFilename);
335  fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column of that line.\n");
336  cdEXIT(EXIT_FAILURE);
337  }
338  fclose(ioDATA);
339 
340  //Sort levels by energy (then by the index in the file if necessary)
341  sort(dBaseStatesEnergy.begin(),dBaseStatesEnergy.end());
342 
343  if( DEBUGSTATE )
344  {
345  fprintf(ioQQQ,"\n\n***Energy levels have been sorted in order of increasing energy.***\n");
346  fprintf(ioQQQ,"Species|File Index|Sorted Index|Energy(WN)|StatWT\n");
347  }
348 
349  /* Store sorted energies in their permanent home */
350  for(DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin();i!=dBaseStatesEnergy.end();i++)
351  {
352  long oldindex = i->second;
353  long index = i - dBaseStatesEnergy.begin();
354  double nrg = i->first;
355  double stwt = dBaseStatesStwt.at(oldindex);
356 
357  if( index >= nMolLevs )
358  break;
359 
360  if( DEBUGSTATE )
361  {
362  fprintf(ioQQQ,"<%s>\t%li\t%li\t%.3f\t%.1f\n",dBaseSpecies[intNS].chLabel,oldindex+1,index+1,nrg,stwt);
363  }
364 
365  dBaseStates[intNS][index].energy().set(nrg,"cm^-1");
366  dBaseStates[intNS][index].g() = stwt;
367  }
368 
369 
370  /* fill in all transition energies, can later overwrite for specific radiative transitions */
371  double fenergyWN = 0.;
372  for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
373  tr!= dBaseTrans[intNS].end(); ++tr)
374  {
375  int ipHi = (*tr).ipHi();
376  int ipLo = (*tr).ipLo();
377  ASSERT(ipHi > ipLo );
378  fenergyWN = dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN();
379  if( fenergyWN <= 0.)
380  {
381  long ipLoFile = dBaseStatesEnergy[ipLo].second;
382  long ipHiFile = dBaseStatesEnergy[ipHi].second;
383 
384 
385  if( DEBUGSTATE )
386  {
387  if( fenergyWN == 0. && dBaseStates[intNS][ipHi].g() != dBaseStates[intNS][ipLo].g() && abs(ipHiFile - ipLoFile) == 1)
388  {
389  /* This likely means that the data source lists adjacent energy levels with the same energy.
390  * Will just use the minimum energy. */
391  if( atmdat.lgStoutPrint )
392  {
393  fprintf( ioQQQ, "\nCaution: A %s transition between adjacent energy levels has zero energy.\n",dBaseSpecies[intNS].chLabel);
394  fprintf( ioQQQ, "The levels may have been sorted by energy since being read from the data files.\n");
395  fprintf( ioQQQ, "The sorted lower and upper levels are %i and %i.\n",ipLo+1,ipHi+1);
396  fprintf( ioQQQ, "The level indices as they appear in the Stout energy level data file, %s, are %li and %li.\n",chNRGFilename,ipLoFile+1,ipHiFile+1);
397  fprintf( ioQQQ, "Verify with the data source that the correct energies are listed in the Stout data file.\n");
398  }
399 
400  }
401  else
402  {
403  fprintf( ioQQQ, "The levels may have been sorted by energy since being read from the data files.\n");
404  fprintf( ioQQQ, "The sorted lower and upper levels are %i and %i.\n",ipLo+1,ipHi+1);
405  fprintf( ioQQQ, "The level indices as they appear in the Stout energy level data file, %s, are %li and %li.\n",chNRGFilename,ipLoFile+1,ipHiFile+1);
406  fprintf( ioQQQ, "Check the data file for missing or duplicate energies.\n");
407  //cdEXIT(EXIT_FAILURE);
408  }
409  }
410  }
411  (*tr).EnergyWN() = (realnum)MAX2(ENERGY_MIN_WN,fenergyWN);
412  (*tr).WLAng() = (realnum)(1e+8/(*tr).EnergyWN()/RefIndex((*tr).EnergyWN()));
413  dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,(*tr).EnergyWN());
414  }
415 
416  /******************************************************
417  ************* Transition Probability File ************
418  ******************************************************/
419  strcat( chTPFilename , ".tp");
420  uncaps( chTPFilename );
421 
422  ioDATA = open_data( chTPFilename, "r" );
423 
424  /* first line is a version number - now confirm that it is valid */
425  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
426  {
427  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chTPFilename );
428  cdEXIT(EXIT_FAILURE);
429  }
430  ipFFmt = 1;
431  iyrread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
432  imoread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
433  idyread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
434 
435  if(( iyrread != iyr ) ||
436  ( imoread != imo ) ||
437  ( idyread != idy ) )
438  {
439  fprintf( ioQQQ,
440  " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chTPFilename );
441  fprintf( ioQQQ,
442  " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
443  iyr, imo , idy ,iyrread, imoread , idyread );
444  cdEXIT(EXIT_FAILURE);
445  }
446 
447  long numtrans = 0;
448  //Count number of transitions
449  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
450  {
451  /* # is comment, *** is end of data */
452  if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
453  {
454  ipFFmt = 1;
455  long n = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
456  if( n < 0 )
457  break;
458  numtrans++;
459  }
460  else if( (chLine[0] == '*' && chLine[1] == '*' ) )
461  {
462  /* stop reading when field of stars encountered.*/
463  break;
464  }
465  }
466  /* now rewind the file so we can read it a second time*/
467  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
468  {
469  fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chTPFilename );
470  cdEXIT(EXIT_FAILURE);
471  }
472  //Skip the magic numbers this time
473  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
474 
475 
476  //Read the first line of data
477  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
478  {
479  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the transition probability file.\n");
480  cdEXIT(EXIT_FAILURE);
481  }
482 
483  double tpData = 0.0;
484  lgSentinelReached = false;
485  static const int intNumCols = 6;
486 
487  /* lgLineStrengthTT functions as a checklist for line strengths
488  * from the various transition types (E1,M2, etc.)
489  * When a line strength is added to the total Aul from a specific
490  * transition type, the corresponding value (based on ipLo, ipHi, and transition type)
491  * is set to true.
492  * lgLineStrengthTT[ipLo][ipHi][k]:
493  * k = 0 => E1
494  * k = 1 => E2
495  * k = 2 => E3
496  * k = 3 => M1
497  * k = 4 => M2
498  * k = 5 => M3
499  */
500  bool ***lgLineStrengthTT;
501  lgLineStrengthTT = (bool ***)MALLOC(nMolLevs *sizeof(bool**));
502  for( int ii = 0; ii < nMolLevs; ii++ )
503  {
504  lgLineStrengthTT[ii] = (bool **)MALLOC(nMolLevs *sizeof(bool*));
505  for( int j = 0; j < nMolLevs; j++ )
506  {
507  lgLineStrengthTT[ii][j] = (bool *)MALLOC(intNumCols *sizeof(bool));
508  }
509  }
510 
511  /* Initialize lgLineStrengthTT values to false */
512  for( int ii = 0; ii < nMolLevs; ii++ )
513  {
514  for( int j = 0; j < nMolLevs; j++ )
515  {
516  for( int k = 0; k < intNumCols; k++ )
517  {
518  lgLineStrengthTT[ii][j][k] = false;
519  }
520  }
521  }
522 
523  if( DEBUGSTATE )
524  {
525  fprintf(ioQQQ,"\nStout Species: %s\n",dBaseSpecies[intNS].chLabel);
526  fprintf(ioQQQ,"Radiative Data File: %s\n",chTPFilename);
527  fprintf(ioQQQ,"Species|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Data Type (A,G,S)|Data\n");
528  }
529 
530  //Read the remaining lines of the transition probability file
531  do
532  {
533  if( chLine[0] == '*' )
534  {
535  lgSentinelReached = true;
536  break;
537  }
538 
539  //Comments start with #, skip them
540  if( chLine[0] != '#' )
541  {
542  /* skip null lines */
543  if( chLine[0] == '\n')
544  continue;
545 
546  if( nMatch("A",chLine) || nMatch("G",chLine) || nMatch("S",chLine) )
547  {
548  /* reset read pointer */
549  ipFFmt = 1;
550 
551  /* save original level for reference, array of levels will be sorted
552  * to become energy ordered */
553  long ipLoInFile = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
554  long ipHiInFile = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
555  long ipLo = ipLoInFile - 1;
556  long ipHi = ipHiInFile - 1;
557  tpData = (double)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
558 
559  if( tpData < atmdat.aulThreshold && nMatch("A",chLine))
560  {
561  // skip these lines
562  continue;
563  }
564 
565  ipLo = getSortedLevel(dBaseStatesEnergy, ipLo);
566  ipHi = getSortedLevel(dBaseStatesEnergy, ipHi);
567 
568  if( ipLo < 0 || ipLo >= nMolLevs || ipHi < 0 || ipHi >= nMolLevs )
569  {
570  // skip these lines
571  continue;
572  }
573 
574  ASSERT( ipLo != ipHi );
575 
576  // swap indices if energy levels were not correctly sorted
577  if( ipHi < ipLo )
578  {
579  long swap = ipHi;
580  ipHi = ipLo;
581  ipLo = swap;
582  }
583 
584  if( DEBUGSTATE )
585  {
586  fprintf(ioQQQ,"<%s>\t%li:%li\t%li:%li\t%c\t%.2e\n",
587  dBaseSpecies[intNS].chLabel,ipLoInFile,ipHiInFile,
588  ipLo+1,ipHi+1,chLine[0],tpData);
589  }
590 
591  if( lgEOL )
592  {
593  fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chTPFilename);
594  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
595  cdEXIT(EXIT_FAILURE);
596  }
597 
598  TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
599 
600  /* If we don't already have the transition in the stack, add it and
601  * zero out Aul */
602  if( !(*tr).hasEmis() )
603  {
604  (*tr).AddLine2Stack();
605  (*tr).Emis().Aul() = 0.;
606  (*tr).Emis().gf() = 0.;
607  }
608 
609  //This means last data column has Aul.
610  if( nMatch("A",chLine) )
611  {
612  if( (*tr).EnergyWN() > ENERGY_MIN_WN )
613  {
614  (*tr).Emis().Aul() += tpData;
615  // use updated total Aul to get gf
616  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
617  }
618  }
619  //This means last data column has gf.
620  else if( nMatch("G",chLine) )
621  {
622  if( (*tr).EnergyWN() > ENERGY_MIN_WN )
623  {
624  (*tr).Emis().gf() += tpData;
625  // use updated total gf to get Aul
626  (*tr).Emis().Aul() = (realnum)eina((*tr).Emis().gf(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
627  }
628  }
629  else if( nMatch("S",chLine) )
630  {
631  static const double BOHR_MAGNETON = ELEM_CHARGE_ESU*H_BAR/2/ELECTRON_MASS/SPEEDLIGHT;
634  /* Data column is composed of line strengths */
635  if( nMatch("E1",chLine) )
636  {
637  if( lgLineStrengthTT[ipLo][ipHi][0] )
638  {
639  AulTTError(chTPFilename,chLine,"E1\0");
640  }
641 
642  /*Convert line strength to Aul for E1 transitions
643  * Aul = 64*Pi^4*e^2*a0^2/3/h*S/gu/WLAng^3 */
644  static const double E1Coeff = 64*powi(PI,4)*pow(ELEM_CHARGE_ESU,2)*pow2(BOHR_RADIUS_CM)/3/HPLANCK/pow3(1e-8);
645  /* E1Coeff = 2.0261e18 */
646 
647  (*tr).Emis().Aul() += E1Coeff*tpData/(*(*tr).Hi()).g()/pow3((*tr).EnergyAng());
648  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
649  lgLineStrengthTT[ipLo][ipHi][0]= true;
650  }
651  else if( nMatch("E2",chLine) )
652  {
653  if( lgLineStrengthTT[ipLo][ipHi][1] )
654  {
655  AulTTError(chTPFilename,chLine,"E2\0");
656  }
657 
658  /*Convert line strength to Aul for E2 transitions
659  * Aul = 64*Pi^6*e^2*a0^4/15/h*S/gu/WLAng^5 */
660  static const double E2Coeff = 64*powi(PI,6)*pow2(ELEM_CHARGE_ESU)*pow4(BOHR_RADIUS_CM)/15/HPLANCK/powi(1e-8,5);
661  /* E2Coeff = 1.1199e18 */
662 
663  (*tr).Emis().Aul() += E2Coeff*tpData/(*(*tr).Hi()).g()/powi((*tr).EnergyAng(),5);
664  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
665  lgLineStrengthTT[ipLo][ipHi][1]= true;
666  }
667  else if( nMatch("E3",chLine) )
668  {
669  if( lgLineStrengthTT[ipLo][ipHi][2] )
670  {
671  AulTTError(chTPFilename,chLine,"E3\0");
672  }
673 
674  /*Convert line strength to Aul for E3 transitions
675  * Aul = 2048*Pi^8*e^2*a0^6/4725/h*S/gu/WLAng^7 */
676  static const double E3Coeff = 2048*powi(PI,8)*pow2(ELEM_CHARGE_ESU)*powi(BOHR_RADIUS_CM,6)/4725/HPLANCK/powi(1e-8,7);
677  /* E3Coeff = 3.1444165e17
678  * Atomic Transition Probabilites of Silicon. A Critical Compilation
679  * Kelleher & Podobedova 2006*/
680 
681  (*tr).Emis().Aul() += E3Coeff*tpData/(*(*tr).Hi()).g()/powi((*tr).EnergyAng(),7);
682  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
683  lgLineStrengthTT[ipLo][ipHi][2]= true;
684  }
685  else if( nMatch("M1",chLine) )
686  {
687  if( lgLineStrengthTT[ipLo][ipHi][3] )
688  {
689  AulTTError(chTPFilename,chLine,"M1\0");
690  }
691 
692  /*Convert line strength to Aul for M1 transitions
693  * Aul = 64*Pi^4*mu_B^2/3/h*S/gu/WLAng^3 */
694  static const double M1Coeff = 64*powi(PI,4)*pow2(BOHR_MAGNETON)/3/HPLANCK/pow3(1e-8);
695  /* M1Coeff = 2.697e13 */
696 
697  (*tr).Emis().Aul() += M1Coeff*tpData/(*(*tr).Hi()).g()/pow3((*tr).EnergyAng());
698  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
699  lgLineStrengthTT[ipLo][ipHi][3]= true;
700  }
701  else if( nMatch("M2",chLine) )
702  {
703  if( lgLineStrengthTT[ipLo][ipHi][4] )
704  {
705  AulTTError(chTPFilename,chLine,"M2\0");
706  }
707 
708  /*Convert line strength to Aul for M2 transitions
709  * Aul = 64*Pi^6*mu_B^2*a0^2/15/h*S/gu/WLAng^5 */
710  static const double M2Coeff = 64*powi(PI,6)*pow2(BOHR_MAGNETON)*pow2(BOHR_RADIUS_CM)/15/HPLANCK/powi(1e-8,5);
711  /* M2Coeff = 1.4909714e13
712  * Tables of Atomic Transition Probabilities for Be and B
713  * Fuhr & Wiese 2010*/
714 
715  (*tr).Emis().Aul() += M2Coeff*tpData/(*(*tr).Hi()).g()/powi((*tr).EnergyAng(),5);
716  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
717  lgLineStrengthTT[ipLo][ipHi][4]= true;
718  }
719  else if( nMatch("M3",chLine) )
720  {
721  if( lgLineStrengthTT[ipLo][ipHi][5] )
722  {
723  AulTTError(chTPFilename,chLine,"M3\0");
724  }
725 
726  /*Convert line strength to Aul for M3 transitions
727  * Aul = 2048*Pi^8*mu_B^2*a0^4/4725/h*S/gu/WLAng^7 */
728  static const double M3Coeff = 2048*powi(PI,8)*pow2(BOHR_MAGNETON)*pow4(BOHR_RADIUS_CM)/4725/HPLANCK/powi(1e-8,7);
729  /* M2Coeff = 4.18610e12
730  * Safronova & Safronova et al 2005*/
731 
732 
733  (*tr).Emis().Aul() += M3Coeff*tpData/(*(*tr).Hi()).g()/powi((*tr).EnergyAng(),7);
734  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
735  lgLineStrengthTT[ipLo][ipHi][5]= true;
736  }
737  else
738  {
739  fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chTPFilename);
740  fprintf( ioQQQ, " The line strength does not list a valid transition type.\n");
741  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
742  cdEXIT(EXIT_FAILURE);
743  }
744  }
745  else
746  {
747  /* The code has already passed a check looking for A, G or S.
748  * Now it can't find any of the three. Insanity. */
749  TotalInsanity();
750  }
751 
752  (*tr).setComment( db_comment_tran_levels( ipLoInFile, ipHiInFile ) );
753  }
754  else
755  {
756  fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chTPFilename);
757  fprintf( ioQQQ, " Data must either be in the form of Aul, gf, or S(line strength) and list"
758  " the data type in the first colum as A, G, or S respectively.");
759  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
760  cdEXIT(EXIT_FAILURE);
761  }
762  }
763  }
764  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
765  if( !lgSentinelReached )
766  {
767  fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chTPFilename);
768  fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column of that line.");
769  cdEXIT(EXIT_FAILURE);
770  }
771  fclose(ioDATA);
772 
773  /******************************************************
774  ************* Collision Data File ********************
775  ******************************************************/
776 
777  strcat( chCOLLFilename , ".coll");
778  uncaps( chCOLLFilename );
779 
780  ioDATA = open_data( chCOLLFilename, "r" );
781 
782  /* first line is a version number - now confirm that it is valid */
783  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
784  {
785  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chCOLLFilename );
786  cdEXIT(EXIT_FAILURE);
787  }
788  ipFFmt = 1;
789  iyrread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
790  imoread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
791  idyread = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
792 
793  if(( iyrread != iyr ) ||
794  ( imoread != imo ) ||
795  ( idyread != idy ) )
796  {
797  fprintf( ioQQQ,
798  " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chCOLLFilename );
799  fprintf( ioQQQ,
800  " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
801  iyr, imo , idy ,iyrread, imoread , idyread );
802  cdEXIT(EXIT_FAILURE);
803  }
804 
805  /****** Could add ability to count number of temperature changes, electron CS, and proton CS ****/
806 
807 
808  //Read the first line of data
809  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
810  {
811  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the collision data file.\n");
812  cdEXIT(EXIT_FAILURE);
813  }
814 
815  /* Malloc space for collision strengths */
816  StoutCollData[intNS].alloc(nMolLevs,nMolLevs,ipNCOLLIDER);
817  for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
818  {
819  for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
820  {
821  for( long k=0; k<ipNCOLLIDER; k++ )
822  {
823  /* initialize all spline variables */
824  StoutCollData[intNS].junk(ipHi,ipLo,k);
825  }
826  }
827  }
828 
829  int numpoints = 0;
830  vector<double> temps;
831  long ipCollider = -1;
832  lgSentinelReached = false;
833 
834  if( DEBUGSTATE )
835  {
836  fprintf(ioQQQ,"\nStout Species: %s\n",dBaseSpecies[intNS].chLabel);
837  fprintf(ioQQQ,"Collision Data File: %s\n",chCOLLFilename);
838  fprintf(ioQQQ,"Species|TEMP|Temperatures (K)\n");
839  fprintf(ioQQQ,"Species|Data Type (CS,RATE)|Collider|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Data\n");
840  }
841 
842  //Read the remaining lines of the collision data file
843  do
844  {
845  /* Stop on *** */
846  if( chLine[0] == '*' )
847  {
848  lgSentinelReached = true;
849  break;
850  }
851 
852  //Comments start with #, skip them
853  if( chLine[0] != '#' )
854  {
855  ipFFmt = 1;
856 
857  /* Skip blank lines */
858  if( chLine[0] == '\n')
859  continue;
860 
861  /* Make all letters of the line upper case so they are case insensitive */
862  caps( chLine );
863 
864  //This is a temperature line
865  if( nMatch("TEMP",chLine) )
866  {
867 
868  if( DEBUGSTATE )
869  {
870  fprintf(ioQQQ,"<%s>\tTEMP\t",dBaseSpecies[intNS].chLabel);
871  }
872 
873  // count number of temperature points
874  ipFFmt = 1;
875  numpoints = 0;
876  while( !lgEOL )
877  {
878  FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
879  if( lgEOL )
880  break;
881  numpoints++;
882  }
883  ASSERT( numpoints > 0 );
884 
885  temps.resize(numpoints);
886  for( int j = 0; j < numpoints; j++ )
887  {
888  temps[j] = 0.;
889  }
890 
891  ipFFmt = 1;
892  for( int j = 0; j < numpoints; j++ )
893  {
894  temps[j] = (double)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
895  if( DEBUGSTATE )
896  {
897  fprintf(ioQQQ,"%.2e\t",temps[j]);
898  }
899  ASSERT( temps[j] > 0 );
900  }
901  if( DEBUGSTATE )
902  {
903  fprintf(ioQQQ,"\n");
904  }
905  }
906  else if( nMatch("CS",chLine) || nMatch("RATE",chLine) )
907  {
908 
909  bool isRate = false;
910  ipFFmt = 1;
911  if( nMatch("RATE", chLine) )
912  isRate = true;
913 
914  if( nMatch( "ELECTRON",chLine ) )
915  {
916  ipCollider = ipELECTRON;
917  }
918  else if( nMatch( "PROTON",chLine ) || nMatch( "H+",chLine ) )
919  {
920  ipCollider = ipPROTON;
921  }
922  else if( nMatch( "HE+2", chLine ) )
923  {
924  ipCollider = ipALPHA;
925  //Move the cursor past the number in the species definition
926  (void)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
927  }
928  else if( nMatch( "HE+ ",chLine ) || nMatch( "HE+\t", chLine) )
929  {
930  ipCollider = ipHE_PLUS;
931  }
932  else if( nMatch( "H2 ",chLine ) || nMatch( "H2\t",chLine ) )
933  {
934  //Move the cursor past the number in the species definition
935  (void)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
936 
937  if( nMatch( "ORTHO",chLine ) )
938  {
939  ipCollider = ipH2_ORTHO;
940  }
941  else if( nMatch( "PARA",chLine ) )
942  {
943  ipCollider = ipH2_PARA;
944  }
945  else
946  {
947  ipCollider = ipH2;
948  }
949  }
950  else if( nMatch( "HE ",chLine ) || nMatch( "HE\t",chLine ) )
951  {
952  ipCollider = ipATOM_HE;
953  }
954  else if( nMatch( "H ",chLine ) || nMatch( "H\t",chLine ) )
955  {
956  ipCollider = ipATOM_H;
957  }
958  else
959  {
960  fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: Each line of the collision data"
961  "file must specify the collider.\n");
962  fprintf( ioQQQ, " Possible colliders are: ELECTRON, PROTON, HE, H,"
963  " HE+2, HE+, H2, H2 ORTHO, H2 PARA\n");
964  cdEXIT(EXIT_FAILURE);
965  }
966 
967  if( temps.empty() )
968  {
969  fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: The collision "
970  "file must specify temperatures before the collision strengths");
971  cdEXIT(EXIT_FAILURE);
972  }
973 
974  long ipLoInFile = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
975  long ipHiInFile = (long)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
976  long ipLo = ipLoInFile-1;
977  long ipHi = ipHiInFile-1;
978 
979  /* Account for reordered energy levels */
980  ipLo = getSortedLevel(dBaseStatesEnergy, ipLo);
981  ipHi = getSortedLevel(dBaseStatesEnergy, ipHi);
982 
983  if( ipLo < 0 || ipLo >= nMolLevs || ipHi < 0 || ipHi >= nMolLevs )
984  {
985  // skip these lines
986  continue;
987  }
988 
989  ASSERT( ipLo != ipHi );
990 
991  // swap indices if energy levels were not correctly sorted
992  if( ipHi < ipLo )
993  {
994  long swap = ipHi;
995  ipHi = ipLo;
996  ipLo = swap;
997  }
998 
999  if( DEBUGSTATE )
1000  {
1001  fprintf(ioQQQ,"<%s>\t%s\t%li\t%li:%li\t%li:%li",
1002  dBaseSpecies[intNS].chLabel,isRate?"RATE":"CS",ipCollider,
1003  ipLoInFile,ipHiInFile,ipLo+1,ipHi+1);
1004  }
1005 
1006  /* Set this as a collision strength not a collision rate coefficient*/
1007  StoutCollData[intNS].lgIsRate(ipHi,ipLo,ipCollider) = isRate;
1008 
1009  ASSERT( numpoints > 0 );
1010  StoutCollData[intNS].setpoints(ipHi,ipLo,ipCollider,numpoints);
1011 
1012  /* Loop over all but one CS value */
1013  for( int j = 0; j < numpoints; j++ )
1014  {
1015  StoutCollData[intNS].temps(ipHi,ipLo,ipCollider)[j] = temps[j];
1016  StoutCollData[intNS].collstrs(ipHi,ipLo,ipCollider)[j] = (double)FFmtRead(chLine,&ipFFmt,sizeof(chLine),&lgEOL);
1017  if( DEBUGSTATE )
1018  {
1019  fprintf(ioQQQ,"\t%.2e",StoutCollData[intNS].collstrs(ipHi,ipLo,ipCollider)[j]);
1020  }
1021  if( StoutCollData[intNS].temps(ipHi,ipLo,ipCollider)[j] <= 0 ||
1022  StoutCollData[intNS].collstrs(ipHi,ipLo,ipCollider)[j] <= 0 )
1023  {
1024  fprintf(ioQQQ,"PROBLEM: A Stout temperature or collision strength is less than or equal to zero.\n");
1025  fprintf(ioQQQ,"Species = %s\tipLo = %li\tipHi = %li\n",dBaseSpecies[intNS].chLabel,ipLo+1,ipHi+1);
1026  fprintf(ioQQQ,"numpoint = %i\tTemp = %e\tCS = %e\n",j,temps[j],
1027  StoutCollData[intNS].collstrs(ipHi,ipLo,ipCollider)[j]);
1028  cdEXIT(EXIT_FAILURE);
1029  }
1030  }
1031  if( DEBUGSTATE )
1032  fprintf(ioQQQ,"\n");
1033  }
1034  else
1035  {
1036  fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chCOLLFilename);
1037  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
1038  cdEXIT(EXIT_FAILURE);
1039  }
1040 
1041  if( lgEOL )
1042  {
1043  fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chCOLLFilename);
1044  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
1045  cdEXIT(EXIT_FAILURE);
1046  }
1047  }
1048  }
1049  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
1050  if( !lgSentinelReached )
1051  {
1052  fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chCOLLFilename);
1053  fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column.");
1054  cdEXIT(EXIT_FAILURE);
1055  }
1056  fclose(ioDATA);
1057 
1058  /* Free up memory from lgLineStrengthTT */
1059  for( int ii = 0; ii < nMolLevs; ii++ )
1060  {
1061  for( int j = 0; j < nMolLevs; j++ )
1062  {
1063  free( lgLineStrengthTT[ii][j] );
1064  }
1065  free( lgLineStrengthTT[ii] );
1066  }
1067  free( lgLineStrengthTT );
1068 
1069  return;
1070 }
1071 
1072 void atmdat_CHIANTI_readin( long intNS, char *chPrefix )
1073 {
1074  DEBUG_ENTRY( "atmdat_CHIANTI_readin()" );
1075 
1076  int intsplinepts,intTranType,intxs;
1077  long int nMolLevs,nMolExpLevs,nElvlcLines,nTheoLevs;// number of experimental and total levels
1078  FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
1079  realnum fstatwt,fenergyWN,fWLAng,fenergy,feinsteina;
1080  double fScalingParam,fEnergyDiff;
1081  const char chCommentChianti = '#';
1082 
1083  char chLine[FILENAME_PATH_LENGTH_2] ,
1084  chEnFilename[FILENAME_PATH_LENGTH_2],
1085  chTraFilename[FILENAME_PATH_LENGTH_2],
1086  chEleColFilename[FILENAME_PATH_LENGTH_2],
1087  chProColFilename[FILENAME_PATH_LENGTH_2];
1088 
1089  bool lgProtonData=false;
1090 
1091  // this is the largest number of levels allowed by the chianti format, I3
1092  static const int MAX_NUM_LEVELS = 999;
1093 
1094  dBaseSpecies[intNS].lgMolecular = false;
1095  dBaseSpecies[intNS].lgLTE = false;
1096 
1097  strcpy( chEnFilename , chPrefix );
1098  strcpy( chTraFilename , chPrefix );
1099  strcpy( chEleColFilename , chPrefix );
1100  strcpy( chProColFilename , chPrefix );
1101 
1102  /*For the CHIANTI DATABASE*/
1103  /*Open the energy levels file*/
1104  strcat( chEnFilename , ".elvlc");
1105  uncaps( chEnFilename );
1106 
1107  /*Open the files*/
1108  if( trace.lgTrace )
1109  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEnFilename);
1110 
1111  fstream elvlcstream;
1112  open_data(elvlcstream, chEnFilename,mode_r);
1113 
1114  /*Open the transition probabilities file*/
1115  strcat( chTraFilename , ".wgfa");
1116  uncaps( chTraFilename );
1117 
1118  /*Open the files*/
1119  if( trace.lgTrace )
1120  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chTraFilename);
1121 
1122  fstream wgfastream;
1123  open_data(wgfastream, chTraFilename,mode_r);
1124 
1125  /*Open the electron collision data*/
1126  strcat( chEleColFilename , ".splups");
1127  uncaps( chEleColFilename );
1128 
1129  /*Open the files*/
1130  if( trace.lgTrace )
1131  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEleColFilename);
1132 
1133  ioElecCollData = open_data( chEleColFilename, "r" );
1134 
1135  /*Open the proton collision data*/
1136  strcat( chProColFilename , ".psplups");
1137  uncaps( chProColFilename );
1138 
1139  /*Open the files*/
1140  if( trace.lgTrace )
1141  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chProColFilename);
1142 
1143  /*We will set a flag here to indicate if the proton collision strengths are available */
1144  if( ( ioProtCollData = open_data( chProColFilename, "r", AS_DATA_ONLY_TRY ) ) != NULL )
1145  {
1146  lgProtonData = true;
1147  }
1148  else
1149  {
1150  lgProtonData = false;
1151  }
1152 
1153  /*Loop finds how many theoretical and experimental levels are in the elvlc file */
1154  //eof_col is used get the first 4 charcters per line to find end of file
1155  const int eof_col = 5;
1156  //length (+1) of the nrg in the elvlc file
1157  const int lvl_nrg_col=16;
1158  //# of columns skipped from the left to get to nrg start
1159  const int lvl_skipto_nrg = 40;
1160  /* # of columns to skip from eof check to nrg start */
1161  const int lvl_eof_to_nrg = lvl_skipto_nrg - eof_col + 1;
1162  //# of columns to skip over the rydberg energy, we don't use it
1163  const int lvl_skip_ryd = 15;
1164  nElvlcLines = 0;
1165  nMolExpLevs = 1;
1166  nTheoLevs = 1;
1167  if (elvlcstream.is_open())
1168  {
1169  int nj = 0;
1170  char otemp[eof_col];
1171  char exptemp[lvl_nrg_col],theotemp[lvl_nrg_col];
1172  double tempexpenergy = 0.,theoenergy = 0.;
1173  /*This loop counts the number of valid rows within the elvlc file
1174  as well as the number of experimental energy levels.*/
1175  while(nj != -1)
1176  {
1177  elvlcstream.get(otemp,eof_col);
1178  nj = atoi(otemp);
1179  if( nj == -1)
1180  break;
1181  nElvlcLines++;
1182 
1183  elvlcstream.seekg(lvl_eof_to_nrg,ios::cur);
1184  elvlcstream.get(exptemp,lvl_nrg_col);
1185  tempexpenergy = (double) atof(exptemp);
1186  if( tempexpenergy != 0.)
1187  nMolExpLevs++;
1188 
1189  elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1190  elvlcstream.get(theotemp,lvl_nrg_col);
1191  theoenergy = (double) atof(theotemp);
1192  if( theoenergy != 0. )
1193  nTheoLevs++;
1194 
1195  elvlcstream.ignore(INT_MAX,'\n');
1196 
1197  }
1198  elvlcstream.seekg(0,ios::beg);
1199  }
1200 
1201  //Sometimes the theoretical chianti level data is incomplete.
1202  //If it is bad use experimental
1203  bool lgChiaBadTheo = false;
1204  if( !atmdat.lgChiantiExp && nTheoLevs < nElvlcLines )
1205  {
1206  lgChiaBadTheo = true;
1207  atmdat.lgChiantiExp = true;
1208  fprintf(ioQQQ,"Warning: The theoretical energy levels for %s are incomplete.",dBaseSpecies[intNS].chLabel);
1209  fprintf(ioQQQ,"Switching to the experimental levels for this species.");
1210  }
1211 
1212  long HighestIndexInFile = -1;
1213 
1214  /* The total number of levels depends on the experimental Chianti switch */
1215  if( atmdat.lgChiantiExp )
1216  {
1217  HighestIndexInFile = nMolExpLevs;
1218  }
1219  else
1220  {
1221  HighestIndexInFile = nElvlcLines;
1222  }
1223 
1224  dBaseSpecies[intNS].numLevels_max = HighestIndexInFile;
1225 
1226  setProperties(dBaseSpecies[intNS]);
1227 
1228  if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
1229  {
1230  // Fe is special case with more levels
1231  nMolLevs = MIN3(HighestIndexInFile, atmdat.nChiantiMaxLevelsFe,MAX_NUM_LEVELS );
1232  }
1233  else
1234  {
1235  nMolLevs = MIN3(HighestIndexInFile, atmdat.nChiantiMaxLevels,MAX_NUM_LEVELS );
1236  }
1237 
1238  if( nMolLevs <= 0 )
1239  {
1240  fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
1241  fprintf( ioQQQ, "The file must be corrupted.\n" );
1242  cdEXIT( EXIT_FAILURE );
1243  }
1244 
1245  //Consider the masterlist specified number of levels as the min. =1 if not specified
1246  long numMasterlist = MIN2( dBaseSpecies[intNS].numLevels_masterlist , HighestIndexInFile );
1247  nMolLevs = MAX2(nMolLevs,numMasterlist);
1248 
1249  if (dBaseSpecies[intNS].setLevels != -1)
1250  {
1251  if (dBaseSpecies[intNS].setLevels > HighestIndexInFile)
1252  {
1253  char chLabelChemical[CHARS_SPECIES] = "";
1254  spectral_to_chemical( chLabelChemical, dBaseSpecies[intNS].chLabel ),
1255  fprintf( ioQQQ,"Using CHIANTI spectrum %s (species: %s) with %li requested, only %li energy levels available.\n",
1256  dBaseSpecies[intNS].chLabel, chLabelChemical, dBaseSpecies[intNS].setLevels, HighestIndexInFile );
1257  nMolLevs = HighestIndexInFile;
1258  }
1259  else
1260  {
1261  nMolLevs = dBaseSpecies[intNS].setLevels;
1262  }
1263  }
1264 
1265  dBaseSpecies[intNS].numLevels_max = nMolLevs;
1266  dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
1267 
1268  if( atmdat.lgChiantiPrint == true)
1269  {
1270  if( atmdat.lgChiantiExp )
1271  {
1272  char chLabelChemical[CHARS_SPECIES] = "";
1273  spectral_to_chemical( chLabelChemical, dBaseSpecies[intNS].chLabel ),
1274  fprintf( ioQQQ,"Using CHIANTI spectrum %s (species: %s) with %li experimental energy levels of %li available.\n",
1275  dBaseSpecies[intNS].chLabel, chLabelChemical, nMolLevs , nMolExpLevs );
1276  }
1277  else
1278  {
1279  char chLabelChemical[CHARS_SPECIES] = "";
1280  spectral_to_chemical( chLabelChemical, dBaseSpecies[intNS].chLabel ),
1281  fprintf( ioQQQ,"Using CHIANTI spectrum %s (species: %s) with %li theoretical energy levels of %li available.\n",
1282  dBaseSpecies[intNS].chLabel, chLabelChemical, nMolLevs , nElvlcLines );
1283  }
1284  }
1285 
1286  /*malloc the States array*/
1287  dBaseStates[intNS].init(dBaseSpecies[intNS].chLabel,nMolLevs);
1288 
1289  /* allocate the Transition array*/
1290  ipdBaseTrans[intNS].reserve(nMolLevs);
1291  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
1292  ipdBaseTrans[intNS].reserve(ipHi,ipHi);
1293  ipdBaseTrans[intNS].alloc();
1294  dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
1295  dBaseTrans[intNS].states() = &dBaseStates[intNS];
1296  dBaseTrans[intNS].chLabel() = dBaseSpecies[intNS].chLabel;
1297  dBaseSpecies[intNS].database = "Chianti";
1298 
1299  int ndBase = 0;
1300  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
1301  {
1302  for( long ipLo = 0; ipLo < ipHi; ipLo++)
1303  {
1304  ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
1305  dBaseTrans[intNS][ndBase].Junk();
1306  dBaseTrans[intNS][ndBase].setLo(ipLo);
1307  dBaseTrans[intNS][ndBase].setHi(ipHi);
1308  ++ndBase;
1309  }
1310  }
1311 
1312  /*Keep track of which levels have experimental data and then create a vector
1313  which relates their indices to the default chianti energy indices.
1314  */
1315  long ncounter = 0;
1316 
1317  //Relate Chianti level indices to a set that only include experimental levels
1318  vector<long> intExperIndex(nElvlcLines,-1);
1319 
1320  DoubleLongPairVector dBaseStatesEnergy;
1321  vector<double> dBaseStatesStwt(HighestIndexInFile,-1.0);
1322  for( long ii = 0; ii < HighestIndexInFile; ii++ )
1323  {
1324  dBaseStatesEnergy.push_back(make_pair(-1.0,ii));
1325  }
1326 
1327  //lvl_skipto_statwt is the # of columns to skip to statwt from left
1328  const int lvl_skipto_statwt = 37;
1329  //lvl_statwt_col is the length (+1) of statwt
1330  const int lvl_statwt_col = 4;
1331  //Read in stat weight and energy
1332 
1333  //Read in nrg levels to see if they are in order
1334  for( long ipLev=0; ipLev<nElvlcLines; ipLev++ )
1335  {
1336  if(elvlcstream.is_open())
1337  {
1338  char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
1339  elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
1340  elvlcstream.get(gtemp,lvl_statwt_col);
1341  fstatwt = (realnum)atof(gtemp);
1342  elvlcstream.get(thtemp,lvl_nrg_col);
1343  fenergy = (double) atof(thtemp);
1344 
1345  if(fstatwt <= 0.)
1346  {
1347  fprintf( ioQQQ, " WARNING: A positive non zero value is expected for the "
1348  "statistical weight but was not found in %s"
1349  " level %li\n", chEnFilename,ipLev);
1350  cdEXIT(EXIT_FAILURE);
1351  }
1352 
1353  if( atmdat.lgChiantiExp )
1354  {
1355  /* Go through the entire level list selectively choosing only experimental level energies.
1356  * Store them, not zeroes, in order using ncounter to count the index.
1357  * Any row on the level list where there is no experimental energy, put a -1 in the relational vector.
1358  * If it is a valid experimental energy level store the new ncounter index.
1359  */
1360 
1361  if( fenergy != 0. || ipLev == 0 )
1362  {
1363  dBaseStatesEnergy.at(ncounter).first = fenergy;
1364  dBaseStatesEnergy.at(ncounter).second = ncounter;
1365  dBaseStatesStwt.at(ncounter) = fstatwt;
1366  intExperIndex.at(ipLev) = ncounter;
1367  ncounter++;
1368  }
1369  else
1370  {
1371  intExperIndex.at(ipLev) = -1;
1372  }
1373  }
1374  else
1375  {
1376  elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1377  elvlcstream.get(obtemp,lvl_nrg_col);
1378  fenergy = (double) atof(obtemp);
1379  if(fenergy != 0. || ipLev == 0)
1380  {
1381  dBaseStatesEnergy.at(ipLev).first = fenergy;
1382  dBaseStatesEnergy.at(ipLev).second = ipLev;
1383  dBaseStatesStwt.at(ipLev) = fstatwt;
1384  }
1385  else
1386  {
1387  dBaseStatesEnergy.at(ipLev).first = -1.;
1388  dBaseStatesEnergy.at(ipLev).second = ipLev;
1389  dBaseStatesStwt.at(ipLev) = -1.;
1390  }
1391  }
1392 
1393  elvlcstream.ignore(INT_MAX,'\n');
1394  }
1395  else
1396  {
1397  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
1398  fclose( ioProtCollData );
1399  cdEXIT(EXIT_FAILURE);
1400  }
1401  }
1402 
1403  elvlcstream.close();
1404 
1405  if( DEBUGSTATE )
1406  {
1407  fprintf(ioQQQ,"\nintExperIndex Vector:\n");
1408  fprintf(ioQQQ,"File Index|Exper Index\n");
1409  for( vector<long>::iterator i = intExperIndex.begin(); i != intExperIndex.end(); i++ )
1410  {
1411  // term on rhs is long in 64 bit, int in 32 bit, print with long format
1412  long iPrt = (i-intExperIndex.begin())+1;
1413  fprintf(ioQQQ,"%li\t%li\n",iPrt,(*i)+1);
1414  }
1415 
1416  for( DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
1417  {
1418  // term on rhs is long in 64 bit, int in 32 bit, print with long format
1419  long iPrt = (i-dBaseStatesEnergy.begin())+1;
1420  fprintf(ioQQQ,"PreSort:%li\t%li\t%f\t%f\n",iPrt,
1421  (i->second)+1,i->first,dBaseStatesStwt.at(i->second));
1422  }
1423  }
1424 
1425  //Sort energy levels
1426  sort(dBaseStatesEnergy.begin(),dBaseStatesEnergy.end());
1427 
1428  if( DEBUGSTATE )
1429  {
1430  for( DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
1431  {
1432  // term on rhs is long in 64 bit, int in 32 bit, print with long format
1433  long iPrt = (i-dBaseStatesEnergy.begin())+1;
1434  if( iPrt > nMolLevs )
1435  break;
1436  fprintf(ioQQQ,"PostSort:%li\t%li\t%f\t%f\n",iPrt,
1437  (i->second)+1,i->first,dBaseStatesStwt.at(i->second));
1438  }
1439 
1440  fprintf(ioQQQ,"\nChianti Species: %s\n",dBaseSpecies[intNS].chLabel);
1441  fprintf(ioQQQ,"Energy Level File: %s\n",chEnFilename);
1442  if( atmdat.lgChiantiExp )
1443  {
1444  fprintf(ioQQQ,"Number of Experimental Energy Levels in File: %li\n",nMolExpLevs);
1445  }
1446  else
1447  {
1448  fprintf(ioQQQ,"Number of Theoretical Energy Levels in File: %li\n",nElvlcLines);
1449  }
1450  fprintf(ioQQQ,"Number of Energy Levels Cloudy is Currently Using: %li\n",nMolLevs);
1451  fprintf(ioQQQ,"Species|File Index|Cloudy Index|StatWT|Energy(WN)\n");
1452  }
1453 
1454  for( DoubleLongPairVector::iterator i=dBaseStatesEnergy.begin(); i != dBaseStatesEnergy.end(); i++ )
1455  {
1456 
1457  long ipLevNew = i - dBaseStatesEnergy.begin();
1458  long ipLevFile = -1;
1459 
1460  if( ipLevNew >= nMolLevs )
1461  break;
1462 
1463  if( atmdat.lgChiantiExp )
1464  {
1465  ipLevFile = getSortedLevel(intExperIndex,ipLevNew);
1466  }
1467  else
1468  {
1469  ipLevFile = i->second;
1470  }
1471 
1472  if( DEBUGSTATE )
1473  {
1474  fprintf(ioQQQ,"<%s>\t%li\t%li\t",dBaseSpecies[intNS].chLabel,ipLevFile+1,ipLevNew+1);
1475  }
1476 
1477  dBaseStates[intNS][ipLevNew].g() = dBaseStatesStwt.at(i->second);
1478  dBaseStates[intNS][ipLevNew].energy().set(i->first,"cm^-1");
1479 
1480  if( DEBUGSTATE )
1481  {
1482  fprintf(ioQQQ,"%.1f\t",dBaseStatesStwt.at(i->second));
1483  fprintf(ioQQQ,"%.3f\n",i->first);
1484  }
1485  }
1486 
1487  // highest energy transition in chianti
1488  dBaseSpecies[intNS].maxWN = 0.;
1489  /* fill in all transition energies, can later overwrite for specific radiative transitions */
1490  for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
1491  tr!= dBaseTrans[intNS].end(); ++tr)
1492  {
1493  int ipHi = (*tr).ipHi();
1494  int ipLo = (*tr).ipLo();
1495  fenergyWN = (realnum)MAX2( ENERGY_MIN_WN , dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN() );
1496 
1497  (*tr).EnergyWN() = fenergyWN;
1498 
1499  (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1500  dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,fenergyWN);
1501  }
1502 
1503  /************************************************************************/
1504  /*** Read in the level numbers, Einstein A and transition wavelength ***/
1505  /************************************************************************/
1506 
1507  //Count the number of rows first
1508  long wgfarows = -1;
1509  if (wgfastream.is_open())
1510  {
1511  int nj = 0;
1512  char otemp[eof_col];
1513  while(nj != -1)
1514  {
1515  wgfastream.get(otemp,eof_col);
1516  wgfastream.ignore(INT_MAX,'\n');
1517  if( otemp[0] == chCommentChianti ) continue;
1518  nj = atoi(otemp);
1519  wgfarows++;
1520  }
1521  wgfastream.seekg(0,ios::beg);
1522  }
1523  else
1524  fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
1525 
1526 
1527  if( DEBUGSTATE )
1528  {
1529  fprintf(ioQQQ,"\n\nTransition Probability File: %s\n",chTraFilename);
1530  fprintf(ioQQQ,"Species|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Wavelength(A)|Ein A\n");
1531  }
1532 
1533 
1534  //line_index_col is the length(+1) of the level indexes in the WGFA file
1535  const int line_index_col = 6;
1536  //line_nrg_to_eina is the # of columns to skip from wavelength to eina in WGFA file
1537  const int line_nrg_to_eina =15;
1538  //line_eina_col is the length(+1) of einsteinA in WGFA
1539  const int line_eina_col = 16;
1540  char lvltemp[line_index_col];
1541  //Start reading WGFA file
1542  if (wgfastream.is_open())
1543  {
1544  for (long ii = 0;ii<wgfarows;ii++)
1545  {
1546  wgfastream.get(lvltemp,line_index_col);
1547  if( lvltemp[0] == chCommentChianti )
1548  {
1549  wgfastream.ignore(INT_MAX,'\n');
1550  continue;
1551  }
1552 
1553  long ipLoInFile = atoi(lvltemp);
1554  wgfastream.get(lvltemp,line_index_col);
1555  long ipHiInFile = atoi(lvltemp);
1556 
1557  // ipLo and ipHi will be manipulated below, want to retain original vals for prints
1558  long ipLo = ipLoInFile - 1;
1559  long ipHi = ipHiInFile - 1;
1560 
1561  if( atmdat.lgChiantiExp )
1562  {
1563  /* If either upper or lower index is -1 in the relational vector,
1564  * skip that line in the wgfa file.
1565  * Otherwise translate the level indices.*/
1566  if( intExperIndex[ipLo] == -1 || intExperIndex[ipHi] == -1 )
1567  {
1568  wgfastream.ignore(INT_MAX,'\n');
1569  continue;
1570  }
1571  else
1572  {
1573  ipLo = getSortedLevel(dBaseStatesEnergy,intExperIndex.at(ipLo));
1574  ipHi = getSortedLevel(dBaseStatesEnergy,intExperIndex.at(ipHi));
1575  }
1576  }
1577  else
1578  {
1579  long testlo = -1, testhi = -1;
1580 
1581  try
1582  {
1583  testlo = getSortedLevel(dBaseStatesEnergy,ipLo);
1584  testhi = getSortedLevel(dBaseStatesEnergy,ipHi);
1585  }
1586  catch ( out_of_range& /* e */ )
1587  {
1588  if( DEBUGSTATE )
1589  {
1590  fprintf(ioQQQ,"NOTE: An out of range exception has occurred"
1591  " reading in data from %s\n",chTraFilename);
1592  fprintf(ioQQQ," The line in the file containing the unidentifiable"
1593  " levels has been ignored.\n");
1594  fprintf(ioQQQ,"There is no reason for alarm."
1595  " This message is just for documentation.\n");
1596  }
1597 
1598  wgfastream.ignore(INT_MAX,'\n');
1599  continue;
1600  }
1601 
1602  if( testlo == -1 || testhi == -1 )
1603  {
1604  wgfastream.ignore(INT_MAX,'\n');
1605  continue;
1606  }
1607  else
1608  {
1609  ipLo = testlo;
1610  ipHi = testhi;
1611  }
1612  }
1613 
1614  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1615  {
1616  // skip these lines
1617  wgfastream.ignore(INT_MAX,'\n');
1618  continue;
1619  }
1620 
1621  if( ipHi == ipLo )
1622  {
1623  fprintf( ioQQQ," WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
1624  wgfastream.ignore(INT_MAX,'\n');
1625  continue;
1626  }
1627 
1628  if( DEBUGSTATE )
1629  {
1630  fprintf(ioQQQ,"<%s>\t%li:%li\t%li:%li\t",dBaseSpecies[intNS].chLabel,ipLoInFile,ipHiInFile,ipLo+1,ipHi+1);
1631  }
1632 
1633  ASSERT( ipHi != ipLo );
1634  ASSERT( ipHi >= 0 );
1635  ASSERT( ipLo >= 0 );
1636 
1637  // sometimes the CHIANTI datafiles list the highest index first as in the middle of these five lines in ne_10.wgfa:
1638  // ...
1639  // 8 10 187.5299 0.000e+00 4.127e+05 3d 2D1.5 - 4s 2S0.5 E2
1640  // 9 10 187.6573 0.000e+00 6.197e+05 3d 2D2.5 - 4s 2S0.5 E2
1641  // 11 10 4842624.0000 1.499e-05 9.423e-06 4p 2P0.5 - 4s 2S0.5 E1
1642  // 1 11 9.7085 1.892e-02 6.695e+11 1s 2S0.5 - 4p 2P0.5 E1
1643  // 2 11 48.5157 6.787e-02 9.618e+10 2s 2S0.5 - 4p 2P0.5 E1
1644  // ...
1645  // so, just set ipHi (ipLo) equal to the max (min) of the two indices.
1646  // NB NB NB it looks like this may depend upon whether one uses observed or theoretical energies.
1647 
1648  //Read in wavelengths
1649  char trantemp[lvl_nrg_col];
1650  wgfastream.get(trantemp,lvl_nrg_col);
1651  fWLAng = (realnum)atof(trantemp);
1652  if( DEBUGSTATE && atmdat.lgChiantiExp)
1653  {
1654  fprintf(ioQQQ,"%.4f\t",fWLAng);
1655  }
1656 
1657  /* \todo 2 CHIANTI labels the H 1 2-photon transition as z wavelength of zero.
1658  * Should we just ignore all of the wavelengths in this file and use the
1659  * difference of level energies instead. */
1660 
1661  if( ipHi < ipLo )
1662  {
1663  long swap = ipHi;
1664  ipHi = ipLo;
1665  ipLo = swap;
1666  }
1667 
1668  /* If the given wavelength is negative, then theoretical energies are being used.
1669  * Take the difference in stored theoretical energies.
1670  * It should equal the absolute value of the wavelength in the wgfa file. */
1671  if( fWLAng <= 0. ) // && !atmdat.lgChiantiExp )
1672  {
1673  //if( fWLAng < 0.)
1674  //fprintf( ioQQQ," WARNING: Negative wavelength for species %6s, indices %3li %3li \n", dBaseSpecies[intNS].chLabel, ipLo, ipHi);
1675  fWLAng = (realnum)(1e8/abs(dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN()));
1676  }
1677 
1678  if( DEBUGSTATE && !atmdat.lgChiantiExp)
1679  {
1680  fprintf(ioQQQ,"%.4f\t",fWLAng);
1681  }
1682  //Skip from end of Wavelength to Einstein A and read in
1683  wgfastream.seekg(line_nrg_to_eina,ios::cur);
1684  wgfastream.get(trantemp,line_eina_col);
1685  feinsteina = (realnum)atof(trantemp);
1686  if( feinsteina == 0. )
1687  {
1688  static bool notPrintedYet = true;
1689  if( notPrintedYet && atmdat.lgChiantiPrint)
1690  {
1691  fprintf( ioQQQ," WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
1692  notPrintedYet = false;
1693  }
1694  wgfastream.ignore(INT_MAX,'\n');
1695  continue;
1696  }
1697  if( DEBUGSTATE )
1698  {
1699  fprintf(ioQQQ,"%.3e\n",feinsteina);
1700  }
1701 
1702  fixit("may need to do something with these rates");
1703  //Read in the rest of the line and look for auto
1704  wgfastream.getline(chLine,INT_MAX);
1705  TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
1706  if( nMatch("auto", chLine) )
1707  {
1708  if( (*tr).hasEmis() )
1709  {
1710  (*tr).Emis().AutoIonizFrac() =
1711  feinsteina/((*tr).Emis().Aul() + feinsteina);
1712  ASSERT( (*tr).Emis().AutoIonizFrac() >= 0. );
1713  ASSERT( (*tr).Emis().AutoIonizFrac() <= 1. );
1714  }
1715  continue;
1716  }
1717 
1718  if( (*tr).hasEmis() )
1719  {
1720  fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_chianti in %s, "
1721  "wavelength=%f\n", chTraFilename,fWLAng);
1722  wgfastream.close();
1723  cdEXIT(EXIT_FAILURE);
1724  }
1725  (*tr).AddLine2Stack();
1726  (*tr).Emis().Aul() = feinsteina;
1727 
1728  fenergyWN = (realnum)(1e+8/fWLAng);
1729 
1730  // TODO::Check the wavelength in the file with the difference in energy levels
1731 
1732  (*tr).EnergyWN() = fenergyWN;
1733  (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1734  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
1735 
1736  (*tr).setComment( db_comment_tran_levels( ipLoInFile, ipHiInFile ) );
1737  }
1738  }
1739  else fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
1740  wgfastream.close();
1741 
1742  /* Malloc space for splines */
1743  AtmolCollSplines[intNS].reserve(nMolLevs);
1744  for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
1745  {
1746  AtmolCollSplines[intNS].reserve(ipHi,nMolLevs);
1747  for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
1748  {
1749  AtmolCollSplines[intNS].reserve(ipHi,ipLo,ipNCOLLIDER);
1750  }
1751  }
1752  AtmolCollSplines[intNS].alloc();
1753 
1754  for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
1755  {
1756  for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
1757  {
1758  for( long k=0; k<ipNCOLLIDER; k++ )
1759  {
1760  /* initialize all spline variables */
1761  AtmolCollSplines[intNS][ipHi][ipLo][k].collspline = NULL;
1762  AtmolCollSplines[intNS][ipHi][ipLo][k].SplineSecDer = NULL;
1763  AtmolCollSplines[intNS][ipHi][ipLo][k].nSplinePts = -1;
1764  AtmolCollSplines[intNS][ipHi][ipLo][k].intTranType = -1;
1765  AtmolCollSplines[intNS][ipHi][ipLo][k].EnergyDiff = BIGDOUBLE;
1766  AtmolCollSplines[intNS][ipHi][ipLo][k].ScalingParam = BIGDOUBLE;
1767  }
1768  }
1769  }
1770 
1771  /************************************/
1772  /*** Read in the collisional data ***/
1773  /************************************/
1774 
1775  // ipCollider 0 is electrons, 1 is protons
1776  for( long ipCollider=0; ipCollider<=1; ipCollider++ )
1777  {
1778  char chFilename[FILENAME_PATH_LENGTH_2];
1779 
1780  if( ipCollider == ipELECTRON )
1781  {
1782  strcpy( chFilename, chEleColFilename );
1783  }
1784  else if( ipCollider == ipPROTON )
1785  {
1786  if( !lgProtonData )
1787  break;
1788  strcpy( chFilename, chProColFilename );
1789  }
1790  else
1791  TotalInsanity();
1792 
1793  /*Dummy string used for convenience*/
1794  strcpy(chLine,"A");
1795 
1796  if( DEBUGSTATE )
1797  {
1798  fprintf(ioQQQ,"\n\nCollision Data File: %s\n",chTraFilename);
1799  fprintf(ioQQQ,"Species|File Index (Lo:Hi)|Cloudy Index (Lo:Hi)|Spline Points\n");
1800  }
1801 
1802  fstream splupsstream;
1803  open_data(splupsstream, chFilename,mode_r);
1804 
1805  //cs_eof_col is the length(+1) of the first column used for finding the end of file
1806  const int cs_eof_col = 4;
1807  //cs_index_col is the length(+1) of the indexes in the CS file
1808  const int cs_index_col = 4;
1809  //cs_trantype_col is the length(+1) of the transition type in the CS file
1810  const int cs_trantype_col = 4;
1811  //cs_values_col is the length(+1) of the other values in the CS file
1812  //including: GF, nrg diff, scaling parameter, and spline points
1813  const int cs_values_col = 11;
1814  //Determine the number of rows in the CS file
1815  if (splupsstream.is_open())
1816  {
1817  int nj = 0;
1818  //splupslines is -1 since the loop runs 1 extra time
1819  long splupslines = -1;
1820  char otemp[cs_eof_col];
1821  while(nj != -1)
1822  {
1823  splupsstream.get(otemp,cs_eof_col);
1824  splupsstream.ignore(INT_MAX,'\n');
1825  nj = atoi(otemp);
1826  splupslines++;
1827  }
1828  splupsstream.seekg(0,ios::beg);
1829 
1830  for (int m = 0;m<splupslines;m++)
1831  {
1832  if( ipCollider == ipELECTRON )
1833  {
1834  splupsstream.seekg(6,ios::cur);
1835  }
1836 
1837  /* level indices */
1838  splupsstream.get(otemp,cs_index_col);
1839  long ipLo = atoi(otemp)-1;
1840  splupsstream.get(otemp,cs_index_col);
1841  long ipHi = atoi(otemp)-1;
1842 
1843  long ipLoFile = ipLo;
1844  long ipHiFile = ipHi;
1845 
1846  /* If either upper or lower index is -1 in the relational vector,
1847  * skip that line in the splups file.
1848  * Otherwise translate the level indices.*/
1849  if( atmdat.lgChiantiExp )
1850  {
1851  if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1852  {
1853  splupsstream.ignore(INT_MAX,'\n');
1854  continue;
1855  }
1856  else
1857  {
1858  ipLo = getSortedLevel(dBaseStatesEnergy,intExperIndex.at(ipLo));
1859  ipHi = getSortedLevel(dBaseStatesEnergy,intExperIndex.at(ipHi));
1860  }
1861  }
1862  else
1863  {
1864  long testlo = -1, testhi = -1;
1865 
1866  /* With level trimming on it is possible that there can be rows that
1867  * have to be skipped when using theoretical
1868  * since the levels no longer exist */
1869  try
1870  {
1871  testlo = getSortedLevel(dBaseStatesEnergy,ipLo);
1872  testhi = getSortedLevel(dBaseStatesEnergy,ipHi);
1873  }
1874  catch ( out_of_range& /* e */ )
1875  {
1876  if( DEBUGSTATE )
1877  {
1878  fprintf(ioQQQ,"NOTE: An out of range exception has occurred"
1879  " reading in data from %s\n",chEleColFilename);
1880  fprintf(ioQQQ," The line in the file containing the unidentifiable"
1881  " levels has been ignored.\n");
1882  fprintf(ioQQQ,"There is no reason for alarm."
1883  " This message is for documentation.\n");
1884  }
1885  splupsstream.ignore(INT_MAX,'\n');
1886  continue;
1887  }
1888 
1889  if( testlo == -1 || testhi == -1 )
1890  {
1891  splupsstream.ignore(INT_MAX,'\n');
1892  continue;
1893  }
1894  else
1895  {
1896  ipLo = testlo;
1897  ipHi = testhi;
1898  }
1899  }
1900 
1901  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1902  {
1903  // skip these transitions
1904  splupsstream.ignore(INT_MAX,'\n');
1905  continue;
1906  }
1907 
1908  if( ipHi < ipLo )
1909  {
1910  long swap = ipHi;
1911  ipHi = ipLo;
1912  ipLo = swap;
1913  }
1914 
1915  if( DEBUGSTATE )
1916  {
1917  fprintf(ioQQQ,"<%s>\t%li:%li\t%li:%li",dBaseSpecies[intNS].chLabel,ipLoFile+1,ipHiFile+1,ipLo+1,ipHi+1);
1918  }
1919 
1920  /*Transition Type*/
1921  splupsstream.get(otemp,cs_trantype_col);
1922  intTranType = atoi(otemp);
1923  char qtemp[cs_values_col];
1924  splupsstream.get(qtemp,cs_values_col);
1925  /*Energy Difference*/
1926  splupsstream.get(qtemp,cs_values_col);
1927  fEnergyDiff = atof(qtemp);
1928  /*Scaling Parameter*/
1929  splupsstream.get(qtemp,cs_values_col);
1930  fScalingParam = atof(qtemp);
1931 
1932  ASSERT( ipLo != ipHi );
1933  ASSERT( ipLo >= 0 && ipLo < nMolLevs );
1934  ASSERT( ipHi >= 0 && ipHi < nMolLevs );
1935  ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline == NULL );
1936  ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer == NULL );
1937 
1938  const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
1939  STATIC_ASSERT(CHIANTI_SPLINE_MAX > CHIANTI_SPLINE_MIN);
1940 
1941  /*We malloc the space here*/
1942  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline =
1943  (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
1944  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer =
1945  (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
1946 
1947  /* always read at least CHIANTI_SPLINE_MIN */
1948  for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
1949  {
1950  //Look at the next character to see if it is the end of line.
1951  char p = splupsstream.peek();
1952  if( p == '\n' )
1953  {
1954  break;
1955  }
1956  else
1957  {
1958  if( intsplinepts >= CHIANTI_SPLINE_MAX )
1959  {
1960  fprintf( ioQQQ, " WARNING: More spline points than expected in %s, indices %3li %3li. Ignoring extras.\n", chFilename, ipHi, ipLo );
1961  break;
1962  }
1963  ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
1964  double temp;
1965  //Store a single spline point then look for more
1966  splupsstream.get(qtemp,cs_values_col);
1967  temp = atof(qtemp);
1968  if( DEBUGSTATE )
1969  {
1970  fprintf(ioQQQ,"\t%.3e",temp);
1971  }
1972  // intTranType == 6 means log10 of numbers have been fit => allow negative numbers
1973  // intTranType < 6 means linear numbers have been fit => negative numbers are unphysical
1974  if( intTranType < 6 )
1975  temp = max( temp, 0. );
1976  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
1977  }
1978  }
1979 
1980  if( DEBUGSTATE )
1981  {
1982  fprintf(ioQQQ,"\n");
1983  }
1984 
1985  ASSERT( intsplinepts > 2 );
1986 
1987  /*The zeroth element contains the number of spline points*/
1988  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].nSplinePts = intsplinepts;
1989  /*Transition type*/
1990  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].intTranType = intTranType;
1991  /*Energy difference between two levels in Rydbergs*/
1992  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].EnergyDiff = fEnergyDiff;
1993  /*Scaling parameter C*/
1994  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;
1995 
1996  /*Once the spline points have been filled,fill the second derivatives structure*/
1997  /*Creating spline points array*/
1998  vector<double> xs (intsplinepts),
1999  spl(intsplinepts),
2000  spl2(intsplinepts);
2001 
2002  for(intxs=0;intxs<intsplinepts;intxs++)
2003  {
2004  double coeff = (double)1/(intsplinepts-1);
2005  xs[intxs] = coeff*intxs;
2006  spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
2007  }
2008 
2009  spline(&xs[0], &spl[0],intsplinepts,2e31,2e31,&spl2[0]);
2010 
2011  /*Filling the second derivative structure*/
2012  for( long ii=0; ii<intsplinepts; ii++)
2013  {
2014  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer[ii] = spl2[ii];
2015  }
2016 
2017  splupsstream.ignore(INT_MAX,'\n');
2018  }
2019  splupsstream.close();
2020  }
2021  }
2022 
2023  // close open file handles
2024  fclose( ioElecCollData );
2025  if( lgProtonData )
2026  fclose( ioProtCollData );
2027 
2028  //Chianti had bad theo level data so we used experimental
2029  //Changing lgChiantiExp back to false so next speices will use theoretical
2030  if( lgChiaBadTheo )
2031  {
2032  atmdat.lgChiantiExp = false;
2033  }
2034 
2035  return;
2036 }
pow4
T pow4(T a)
Definition: cddefines.h:999
getSortedLevel
long getSortedLevel(const T &input, long oldLev)
Definition: atmdat_chianti.cpp:33
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:751
StoutCollData
vector< StoutCollArray > StoutCollData
Definition: taulines.cpp:21
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
DEBUGSTATE
static const bool DEBUGSTATE
Definition: atmdat_chianti.cpp:70
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:296
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
input
t_input input
Definition: input.cpp:12
BIGDOUBLE
const double BIGDOUBLE
Definition: cpu.h:249
parse_species.h
eina
double eina(double gf, double enercm, double gup)
Definition: lines_service.cpp:75
nMatch
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:537
t_atmdat::nStoutMaxLevelsFe
long nStoutMaxLevelsFe
Definition: atmdat.h:410
DoubleLongPairVector
vector< pair< double, long > > DoubleLongPairVector
Definition: atmdat_chianti.cpp:11
ipdBaseTrans
vector< multi_arr< int, 2 > > ipdBaseTrans
Definition: taulines.cpp:17
atmdat_STOUT_readin
void atmdat_STOUT_readin(long intNS, char *chFileName)
Definition: atmdat_chianti.cpp:74
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:167
t_atmdat::lgChiantiPrint
bool lgChiantiPrint
Definition: atmdat.h:378
pow3
T pow3(T a)
Definition: cddefines.h:992
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
MIN2
#define MIN2(a, b)
Definition: cddefines.h:807
spectral_to_chemical
STATIC void spectral_to_chemical(char *chLabelChemical, char *chLabel, long &nelem, long &IonStg)
Definition: species.cpp:952
uncaps
void uncaps(char *chCard)
Definition: service.cpp:287
trace
t_trace trace
Definition: trace.cpp:5
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:556
lines_service.h
spline
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
Definition: thirdparty.h:150
tolower
char tolower(char c)
Definition: cddefines.h:734
mode_r
const ios_base::openmode mode_r
Definition: cpu.h:267
t_atmdat::aulThreshold
static const double aulThreshold
Definition: atmdat.h:437
swap
void swap(count_ptr< T > &a, count_ptr< T > &b)
Definition: count_ptr.h:82
energy
double energy(const genericState &gs)
Definition: generic_state.cpp:49
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
AulTTError
static void AulTTError(const char chFilename[], const char chLine[], const char TT[])
Definition: atmdat_chianti.cpp:20
setProperties
void setProperties(species &sp)
Definition: parse_species.cpp:103
t_atmdat::nStoutMaxLevels
long nStoutMaxLevels
Definition: atmdat.h:413
realnum
float realnum
Definition: cddefines.h:124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
max
long max(int a, long b)
Definition: cddefines.h:821
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:484
powi
double powi(double, long int)
Definition: service.cpp:690
t_atmdat::nChiantiMaxLevels
long nChiantiMaxLevels
Definition: atmdat.h:386
AtmolCollSplines
vector< multi_arr< CollSplinesArray, 3 > > AtmolCollSplines
Definition: taulines.cpp:20
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:617
t_atmdat::nDefaultPhotoLevelsFe
const long nDefaultPhotoLevelsFe
Definition: atmdat.h:360
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:92
db_comment_tran_levels
string db_comment_tran_levels(long ipLoFile, long ipHiFile)
Definition: atmdat.h:515
pow2
T pow2(T a)
Definition: cddefines.h:985
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:16
ENERGY_MIN_WN
const double ENERGY_MIN_WN
Definition: atmdat_chianti.cpp:72
MIN3
#define MIN3(a, b, c)
Definition: cddefines.h:812
dBaseSpecies
vector< species > dBaseSpecies
Definition: taulines.cpp:15
MAX2
#define MAX2(a, b)
Definition: cddefines.h:828
t_atmdat::lgStoutPrint
bool lgStoutPrint
Definition: atmdat.h:406
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
atmdat_CHIANTI_readin
void atmdat_CHIANTI_readin(long intNS, char *chFileName)
Definition: atmdat_chianti.cpp:1072
pow
double pow(double x, int i)
Definition: cddefines.h:782
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:72
caps
void caps(char *chCard)
Definition: service.cpp:304
fixit
#define fixit(a)
Definition: cddefines.h:416
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
t_atmdat::nChiantiMaxLevelsFe
long nChiantiMaxLevelsFe
Definition: atmdat.h:384
t_atmdat::lgChiantiExp
bool lgChiantiExp
Definition: atmdat.h:380
STATIC_ASSERT
#define STATIC_ASSERT(x)
Definition: cddefines.h:140
ipH2
Definition: collision.h:19
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:393
Generated by   doxygen 1.8.5