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