cloudy  trunk
atmdat_chianti.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cddefines.h"
4 #include "lines_service.h"
5 #include "physconst.h"
6 #include "taulines.h"
7 #include "trace.h"
8 #include "string.h"
9 #include "thirdparty.h"
10 #include "rfield.h"
11 #include "atmdat.h"
12 
13 static const bool DEBUGSTATE = false;
14 // minimum energy difference (wavenumbers) we will accept
15 const double ENERGY_MIN_WN = 1e-10;
16 
17 void atmdat_STOUT_readin( long intNS, char *chPrefix )
18 {
19  DEBUG_ENTRY( "atmdat_STOUT_readin()" );
20 
21  long int nMolLevs;
22  char chLine[FILENAME_PATH_LENGTH_2],
23  chNRGFilename[FILENAME_PATH_LENGTH_2],
24  chTPFilename[FILENAME_PATH_LENGTH_2],
25  chCOLLFilename[FILENAME_PATH_LENGTH_2];
26 
27  static const int MAX_NUM_LEVELS = 999;
28 
29  dBaseSpecies[intNS].lgMolecular = false;
30  dBaseSpecies[intNS].lgLTE = false;
31  dBaseSpecies[intNS].lgLAMDA = false;
32 
33  strcpy( chNRGFilename , chPrefix );
34  strcpy( chTPFilename , chPrefix );
35  strcpy( chCOLLFilename , chPrefix );
36 
37  /*Open the energy levels file*/
38  strcat( chNRGFilename , ".nrg");
39  uncaps( chNRGFilename );
40  if(DEBUGSTATE)
41  printf("The data file is %s \n",chNRGFilename);
42 
43  /*Open the files*/
44  if( trace.lgTrace )
45  fprintf( ioQQQ," atmdat_STOUT_readin opening %s:",chNRGFilename);
46 
47  FILE *ioDATA;
48  ioDATA = open_data( chNRGFilename, "r" );
49  long int i = 0;
50  bool lgEOL = false;
51  long index = 0;
52  double nrg = 0.0;
53  double stwt = 0.0;
54 
55  /* first line is a version number - now confirm that it is valid */
56  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
57  {
58  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chNRGFilename );
59  cdEXIT(EXIT_FAILURE);
60  }
61  i = 1;
62  const long int iyr = 11, imo=10 , idy = 14;
63  long iyrread, imoread , idyread;
64  iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
65  imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
66  idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
67 
68  if(( iyrread != iyr ) ||
69  ( imoread != imo ) ||
70  ( idyread != idy ) )
71  {
72  fprintf( ioQQQ,
73  " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chNRGFilename );
74  fprintf( ioQQQ,
75  " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
76  iyr, imo , idy ,iyrread, imoread , idyread );
77  cdEXIT(EXIT_FAILURE);
78  }
79 
80  nMolLevs = 0;
81  //Count number of levels
82  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
83  {
84  /* # - comment, *** ends data */
85  if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
86  {
87  i = 1;
88  long n = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
89  if( n < 0 )
90  break;
91  nMolLevs++;
92  }
93  }
94 
95  /* now rewind the file so we can read it a second time*/
96  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
97  {
98  fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chNRGFilename );
99  cdEXIT(EXIT_FAILURE);
100  }
101  //Skip the magic numbers this time
102  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
103 
104  long HighestIndexInFile = nMolLevs;
105 
106  nMolLevs = MIN3(atmdat.nStoutMaxLevels,nMolLevs, MAX_NUM_LEVELS );
107 
108  if( atmdat.lgStoutPrint == true)
109  {
110  fprintf( ioQQQ," Using STOUT model %s with %li levels of %li available.\n",
111  dBaseSpecies[intNS].chLabel , nMolLevs , HighestIndexInFile );
112  }
113 
114  dBaseSpecies[intNS].numLevels_max = nMolLevs;
115  dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
116 
117  /*Resize the States array*/
118  dBaseStates[intNS].resize(nMolLevs);
119  /*Zero out the maximum wavenumber for each species */
120  dBaseSpecies[intNS].maxWN = 0.;
121 
122  /* allocate the Transition array*/
123  ipdBaseTrans[intNS].reserve(nMolLevs);
124  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
125  ipdBaseTrans[intNS].reserve(ipHi,ipHi);
126  ipdBaseTrans[intNS].alloc();
127  dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
128  dBaseTrans[intNS].states() = &dBaseStates[intNS];
129  dBaseTrans[intNS].chLabel() = "Stout";
130 
131  //This is creating transitions that we don't have collision data for
132  //Maybe use gbar or keep all of the Fe 2 even if it was assumed (1e-5)
133  int ndBase = 0;
134  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
135  {
136  for( long ipLo = 0; ipLo < ipHi; ipLo++)
137  {
138  ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
139  dBaseTrans[intNS][ndBase].Junk();
140  dBaseTrans[intNS][ndBase].setLo(ipLo);
141  dBaseTrans[intNS][ndBase].setHi(ipHi);
142  ++ndBase;
143  }
144  }
145 
146 
147  /* Check for an end of file sentinel */
148  bool lgSentinelReached = false;
149 
150  //Read the first line of data
151  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
152  {
153  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the energy level file.\n");
154  cdEXIT(EXIT_FAILURE);
155  }
156  //Read the remaining lines of the energy level file
157  do
158  {
159  i = 1;
160 
161  /* Stop on *** */
162  if( chLine[0] == '*' )
163  {
164  lgSentinelReached = true;
165  break;
166  }
167  //Comments start with #, skip them
168  if( chLine[0] != '#' )
169  {
170  /* Skip blank lines */
171  if( chLine[0] == '\n')
172  continue;
173 
174  index = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) -1 ;
175 
176  if( index < 0 )
177  {
178  fprintf( ioQQQ, " PROBLEM An energy level index was less than 1 in file %s\n",chNRGFilename);
179  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
180  cdEXIT(EXIT_FAILURE);
181  }
182 
183  nrg = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
184  stwt = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
185 
186  if( lgEOL )
187  {
188  fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chNRGFilename);
189  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
190  cdEXIT(EXIT_FAILURE);
191  }
192 
193  if( index >= nMolLevs )
194  {
195  //Skip unused levels
196  continue;
197  }
198 
199  dBaseStates[intNS][index].energy().set(nrg,"cm^-1");
200  dBaseStates[intNS][index].g() = stwt;
201 
202  /*information needed for label*/
203  strcpy(dBaseStates[intNS][index].chLabel(), " ");
204  strncpy(dBaseStates[intNS][index].chLabel(),dBaseSpecies[intNS].chLabel, 4);
205  dBaseStates[intNS][index].chLabel()[4] = '\0';
206  }
207  }
208  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
209  if( !lgSentinelReached )
210  {
211  fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chNRGFilename);
212  fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column of that line.\n");
213  cdEXIT(EXIT_FAILURE);
214  }
215  fclose(ioDATA);
216 
217  if( DEBUGSTATE )
218  {
219  /*Test whether energy levels are read in properly*/
220  printf("DEBUG STOUT ENERGY LEVELS:\n");
221  for( int i = 0; i< nMolLevs; i++ )
222  {
223  printf("DEBUG\t%i\t%11.4e\t%3.1f\n",i+1,dBaseStates[intNS][i].energy().WN(),dBaseStates[intNS][i].g());
224  }
225  }
226 
227  /* fill in all transition energies, can later overwrite for specific radiative transitions */
228  double fenergyWN = 0.;
229  for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
230  tr!= dBaseTrans[intNS].end(); ++tr)
231  {
232  int ipHi = (*tr).ipHi();
233  int ipLo = (*tr).ipLo();
234  ASSERT(ipHi > ipLo );
235  fenergyWN = dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN();
236  if( fenergyWN <= 0.)
237  {
238  printf("\nWARNING: The %s transition between level %i and %i has zero or negative energy.\n",
239  dBaseStates[intNS][ipHi].chLabel(),ipLo+1,ipHi+1);
240  printf("Check the Stout energy level data file (*.nrg) of this species for missing or duplicate energies.\n");
241  //cdEXIT(EXIT_FAILURE);
242  }
243  (*tr).EnergyWN() = (realnum)MAX2(ENERGY_MIN_WN,fenergyWN);
244  (*tr).WLAng() = (realnum)(1e+8/(*tr).EnergyWN()/RefIndex((*tr).EnergyWN()));
245  dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,(*tr).EnergyWN());
246  }
247 
248  /******************************************************
249  ************* Transition Probability File ************
250  ******************************************************/
251  strcat( chTPFilename , ".tp");
252  uncaps( chTPFilename );
253 
254  ioDATA = open_data( chTPFilename, "r" );
255 
256  /* first line is a version number - now confirm that it is valid */
257  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
258  {
259  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chTPFilename );
260  cdEXIT(EXIT_FAILURE);
261  }
262  i = 1;
263  iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
264  imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
265  idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
266 
267  if(( iyrread != iyr ) ||
268  ( imoread != imo ) ||
269  ( idyread != idy ) )
270  {
271  fprintf( ioQQQ,
272  " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chTPFilename );
273  fprintf( ioQQQ,
274  " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
275  iyr, imo , idy ,iyrread, imoread , idyread );
276  cdEXIT(EXIT_FAILURE);
277  }
278 
279  long numtrans = 0;
280  //Count number of transitions
281  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
282  {
283  /* # is comment, *** is end of data */
284  if( chLine[0] != '#' && chLine[0] != '\n' && chLine[0] != '*' )
285  {
286  i = 1;
287  long n = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
288  if( n < 0 )
289  break;
290  numtrans++;
291  }
292  }
293  /* now rewind the file so we can read it a second time*/
294  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
295  {
296  fprintf( ioQQQ, " atmdat_STOUT_readin could not rewind %s.\n", chTPFilename );
297  cdEXIT(EXIT_FAILURE);
298  }
299  //Skip the magic numbers this time
300  read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
301 
302 
303  //Read the first line of data
304  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
305  {
306  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the transition probability file.\n");
307  cdEXIT(EXIT_FAILURE);
308  }
309 
310  long ipLo = 0;
311  long ipHi = 0;
312  double Aul = 0.0;
313  lgSentinelReached = false;
314  //Read the remaining lines of the transition probability file
315  do
316  {
317  if( chLine[0] == '*' )
318  {
319  lgSentinelReached = true;
320  break;
321  }
322 
323  //Comments start with #, skip them
324  if( chLine[0] != '#' )
325  {
326  i = 1;
327 
328  /* skip null lines */
329  if( chLine[0] == '\n')
330  continue;
331 
332  //This means last data column has Aul.
333  if( nMatch("A",chLine) )
334  {
335  /* reset read pointer */
336  i = 1;
337 
338  ipLo= (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
339  ipHi = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
340  Aul = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
341  if( lgEOL )
342  {
343  fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chTPFilename);
344  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
345  cdEXIT(EXIT_FAILURE);
346  }
347 
348  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
349  {
350  // skip these lines
351  continue;
352  }
353 
354  TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
355  if( (*tr).hasEmis() )
356  {
357  fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_STOUT_readin in %s, "
358  "wavelength=%f\n", chTPFilename,dBaseStates[intNS][ipHi].energy().WN()
359  - dBaseStates[intNS][ipLo].energy().WN());
360  cdEXIT(EXIT_FAILURE);
361  }
362 
363  if( (*tr).EnergyWN() > ENERGY_MIN_WN )
364  {
365  (*tr).AddLine2Stack();
366  (*tr).Emis().Aul() = Aul;
367  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
368  }
369  }
370  else
371  {
372  fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chTPFilename);
373  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
374  cdEXIT(EXIT_FAILURE);
375  }
376  }
377  }
378  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
379  if( !lgSentinelReached )
380  {
381  fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chTPFilename);
382  fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column of that line.");
383  cdEXIT(EXIT_FAILURE);
384  }
385  fclose(ioDATA);
386 
387  if( DEBUGSTATE )
388  {
389  /*Test whether stout A's are read in properly */
390  printf("DEBUG STOUT A's:\n");
391  for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
392  tr!= dBaseTrans[intNS].end(); ++tr)
393  {
394  printf("DEBUG.STOUT.A:\t%i\t%i\t%8.2e\n",(*tr).ipLo()+1,(*tr).ipHi()+1,(*tr).Emis().Aul());
395  }
396  }
397 
398 
399 
400  /******************************************************
401  ************* Collision Data File ********************
402  ******************************************************/
403 
404  strcat( chCOLLFilename , ".coll");
405  uncaps( chCOLLFilename );
406 
407  ioDATA = open_data( chCOLLFilename, "r" );
408 
409  /* first line is a version number - now confirm that it is valid */
410  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
411  {
412  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of %s.\n", chCOLLFilename );
413  cdEXIT(EXIT_FAILURE);
414  }
415  i = 1;
416  iyrread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
417  imoread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
418  idyread = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
419 
420  if(( iyrread != iyr ) ||
421  ( imoread != imo ) ||
422  ( idyread != idy ) )
423  {
424  fprintf( ioQQQ,
425  " PROBLEM atmdat_STOUT_readin: the version of %s is not the current version.\n", chCOLLFilename );
426  fprintf( ioQQQ,
427  " atmdat_STOUT_readin: I expected the magic numbers %li %li %li but found %li %li %li.\n",
428  iyr, imo , idy ,iyrread, imoread , idyread );
429  cdEXIT(EXIT_FAILURE);
430  }
431 
432  /****** Could add ability to count number of temperature changes, electron CS, and proton CS ****/
433 
434 
435  //Read the first line of data
436  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
437  {
438  fprintf( ioQQQ, " atmdat_STOUT_readin could not read first line of the collision data file.\n");
439  cdEXIT(EXIT_FAILURE);
440  }
441 
442  /* Malloc space for collision strengths */
443  StoutCollData[intNS] = (StoutColls***)MALLOC(nMolLevs *sizeof(StoutColls**));
444  for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
445  {
446  StoutCollData[intNS][ipHi] = (StoutColls **)MALLOC((unsigned long)(nMolLevs)*sizeof(StoutColls *));
447  for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
448  {
449  StoutCollData[intNS][ipHi][ipLo] =
450  (StoutColls *)MALLOC((unsigned long)(ipNCOLLIDER)*sizeof(StoutColls ));
451 
452  for( long k=0; k<ipNCOLLIDER; k++ )
453  {
454  /* initialize all spline variables */
455  StoutCollData[intNS][ipHi][ipLo][k].ntemps = -1;
456  StoutCollData[intNS][ipHi][ipLo][k].temps = NULL;
457  StoutCollData[intNS][ipHi][ipLo][k].collstrs = NULL;
458  StoutCollData[intNS][ipHi][ipLo][k].lgIsRate = false;
459  }
460  }
461  }
462 
463  ipLo = 0;
464  ipHi = 0;
465  int numpoints = 0;
466  double *temps = NULL;
467  long ipCollider = -1;
468  lgSentinelReached = false;
469 
470  //Read the remaining lines of the collision data file
471  do
472  {
473  /* Stop on *** */
474  if( chLine[0] == '*' )
475  {
476  lgSentinelReached = true;
477  break;
478  }
479 
480  //Comments start with #, skip them
481  if( chLine[0] != '#' )
482  {
483  i = 1;
484 
485  /* Skip blank lines */
486  if( chLine[0] == '\n')
487  continue;
488 
489  //This is a temperature line
490  if( nMatch("TEMP",chLine) )
491  {
492  if( temps != NULL)
493  free( temps );
494 
495  i = 1;
496  numpoints = (int)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
497  ASSERT( numpoints > 0 );
498 
499  temps = (double*)MALLOC(numpoints*sizeof(double));
500  memset( temps, 0, (unsigned long)(numpoints)*sizeof(double) );
501  for( int j = 0; j < numpoints; j++ )
502  {
503  temps[j] = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
504  ASSERT( temps[j] > 0 );
505  }
506  }
507  else if( nMatch("CS",chLine) || nMatch("RATE",chLine) )
508  {
509 
510  bool isRate = false;
511  if( nMatch("RATE", chLine) )
512  isRate = true;
513 
514  if( nMatch( "ELECTRON",chLine ) )
515  {
516  ipCollider = ipELECTRON;
517  }
518  else if( nMatch( "PROTON",chLine ) )
519  {
520  ipCollider = ipPROTON;
521  }
522  else
523  {
524  fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: Each line of the collision data"
525  "file must specify the collider.\n");
526  fprintf( ioQQQ, " Possible colliders are: ELECTRON, PROTON\n");
527  cdEXIT(EXIT_FAILURE);
528  }
529 
530  if( temps == NULL )
531  {
532  fprintf( ioQQQ, " PROBLEM atmdat_STOUT_readin: The collision "
533  "file must specify temperatures before the collision strengths");
534  cdEXIT(EXIT_FAILURE);
535  }
536 
537  i = 1;
538  ipLo= (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
539  ipHi = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL) - 1;
540 
541  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
542  {
543  // skip these lines
544  continue;
545  }
546 
547  /* Set this as a collision strength not a collision rate */
548  StoutCollData[intNS][ipHi][ipLo][ipCollider].lgIsRate = isRate;
549 
550  StoutCollData[intNS][ipHi][ipLo][ipCollider].ntemps = numpoints;
551  ASSERT( numpoints > 0 );
552 
553  /*malloc the space here for the temps and collision strengths*/
554  StoutCollData[intNS][ipHi][ipLo][ipCollider].temps =
555  (double *)MALLOC((unsigned long)(numpoints)*sizeof(double));
556  StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs =
557  (double *)MALLOC((unsigned long)(numpoints)*sizeof(double));
558  /* Loop over all but one CS value */
559  for( int j = 0; j < numpoints; j++ )
560  {
561  StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] = temps[j];
562  ASSERT( StoutCollData[intNS][ipHi][ipLo][ipCollider].temps[j] > 0 );
563  StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] = (double)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
564  ASSERT( StoutCollData[intNS][ipHi][ipLo][ipCollider].collstrs[j] > 0 );
565  }
566  }
567  else
568  {
569  fprintf( ioQQQ, " PROBLEM File %s contains an invalid line.\n",chCOLLFilename);
570  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
571  cdEXIT(EXIT_FAILURE);
572  }
573 
574  if( lgEOL )
575  {
576  fprintf( ioQQQ, " PROBLEM End of line reached prematurely in file %s\n",chCOLLFilename);
577  fprintf( ioQQQ, " The line being read is between the braces {%.*s}\n",int(strlen(chLine)-1),chLine);
578  cdEXIT(EXIT_FAILURE);
579  }
580  }
581  }
582  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL );
583  if( !lgSentinelReached )
584  {
585  fprintf( ioQQQ, " PROBLEM End of data sentinel was not reached in file %s\n",chCOLLFilename);
586  fprintf( ioQQQ, " Check that stars (*****) appear after the last line of data and start in the first column.");
587  cdEXIT(EXIT_FAILURE);
588  }
589  free(temps);
590  fclose(ioDATA);
591 
592  return;
593 }
594 
595 void atmdat_CHIANTI_readin( long intNS, char *chPrefix )
596 {
597  DEBUG_ENTRY( "atmdat_CHIANTI_readin()" );
598 
599  int intsplinepts,intTranType,intxs;
600  long int nMolLevs,nMolExpLevs,nElvlcLines;// number of experimental and total levels
601  FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
602  realnum fstatwt,fenergyWN,fWLAng,fenergy,feinsteina;
603  double fScalingParam,fEnergyDiff,*xs,*spl,*spl2;
604 
605  char chLine[FILENAME_PATH_LENGTH_2] ,
606  chEnFilename[FILENAME_PATH_LENGTH_2],
607  chTraFilename[FILENAME_PATH_LENGTH_2],
608  chEleColFilename[FILENAME_PATH_LENGTH_2],
609  chProColFilename[FILENAME_PATH_LENGTH_2];
610 
611  realnum *dBaseStatesEnergy;
612  long int *intNewIndex=NULL,*intOldIndex=NULL, *intExperIndex=NULL;
613  int interror;
614  bool lgProtonData=false,lgEneLevOrder;
615 
616  // this is the largest number of levels allowed by the chianti format, I3
617  static const int MAX_NUM_LEVELS = 999;
618 
619  dBaseSpecies[intNS].lgMolecular = false;
620  dBaseSpecies[intNS].lgLAMDA = false;
621  dBaseSpecies[intNS].lgLTE = false;
622 
623  strcpy( chEnFilename , chPrefix );
624  strcpy( chTraFilename , chPrefix );
625  strcpy( chEleColFilename , chPrefix );
626  strcpy( chProColFilename , chPrefix );
627 
628  /*For the CHIANTI DATABASE*/
629  /*Open the energy levels file*/
630  strcat( chEnFilename , ".elvlc");
631  uncaps( chEnFilename );
632  if(DEBUGSTATE)
633  printf("The data file is %s \n",chEnFilename);
634 
635  /*Open the files*/
636  if( trace.lgTrace )
637  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEnFilename);
638 
639  fstream elvlcstream;
640  open_data(elvlcstream, chEnFilename,mode_r);
641 
642  /*Open the transition probabilities file*/
643  strcat( chTraFilename , ".wgfa");
644  uncaps( chTraFilename );
645  if(DEBUGSTATE)
646  printf("The data file is %s \n",chTraFilename);
647 
648  /*Open the files*/
649  if( trace.lgTrace )
650  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chTraFilename);
651 
652  fstream wgfastream;
653  open_data(wgfastream, chTraFilename,mode_r);
654 
655  /*Open the electron collision data*/
656  strcat( chEleColFilename , ".splups");
657  uncaps( chEleColFilename );
658  if(DEBUGSTATE)
659  printf("The data file is %s \n",chEleColFilename);
660  /*Open the files*/
661  if( trace.lgTrace )
662  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEleColFilename);
663 
664  ioElecCollData = open_data( chEleColFilename, "r" );
665 
666  /*Open the proton collision data*/
667  strcat( chProColFilename , ".psplups");
668  uncaps( chProColFilename );
669  if(DEBUGSTATE)
670  printf("The data file is %s \n",chProColFilename);
671  /*Open the files*/
672  if( trace.lgTrace )
673  fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chProColFilename);
674 
675  /*We will set a flag here to indicate if the proton collision strengths are available */
676  if( ( ioProtCollData = fopen( chProColFilename , "r" ) ) != NULL )
677  {
678  lgProtonData = true;
679  //fclose( ioProtCollData );
680  //ioProtCollData = NULL;
681  }
682  else
683  {
684  lgProtonData = false;
685  }
686 
687  /*Loop finds how many theoretical and experimental levels are in the elvlc file */
688  //eof_col is used get the first 4 charcters per line to find end of file
689  const int eof_col = 5;
690  //length (+1) of the nrg in the elvlc file
691  const int lvl_nrg_col=16;
692  //# of columns skipped from the left to get to nrg start
693  const int lvl_skipto_nrg = 40;
694  /* # of columns to skip from eof check to nrg start */
695  const int lvl_eof_to_nrg = lvl_skipto_nrg - eof_col + 1;
696  nElvlcLines = 0;
697  nMolExpLevs = 1;
698  lgEneLevOrder = true;
699  if (elvlcstream.is_open())
700  {
701  int nj = 0;
702  char otemp[eof_col];
703  char exptemp[lvl_nrg_col];
704  double tempexpenergy = 0.;
705  /*This loop counts the number of valid rows within the elvlc file
706  as well as the number of experimental energy levels.*/
707  while(nj != -1)
708  {
709  elvlcstream.get(otemp,eof_col);
710  nj = atoi(otemp);
711  if( nj == -1)
712  break;
713  nElvlcLines++;
714 
715  elvlcstream.seekg(lvl_eof_to_nrg,ios::cur);
716  elvlcstream.get(exptemp,lvl_nrg_col);
717  tempexpenergy = (realnum) atof(exptemp);
718  if( tempexpenergy != 0.)
719  nMolExpLevs++;
720 
721  elvlcstream.ignore(INT_MAX,'\n');
722 
723  }
724  elvlcstream.seekg(0,ios::beg);
725  }
726 
727  if(DEBUGSTATE)
728  {
729  printf("DEBUG: The file %s contains %li experimental energy levels of the %li total levels. \n",chEnFilename,nMolExpLevs,nElvlcLines);
730  }
731 
732  /* The total number of levels depends on the experimental Chianti switch */
733  if( atmdat.lgChiantiExp )
734  {
735  if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
736  {
737  // Fe is special case with more levels
738  nMolLevs = MIN3(nMolExpLevs, atmdat.nChiantiMaxLevelsFe, MAX_NUM_LEVELS );
739  }
740  else
741  {
742  nMolLevs = MIN3(nMolExpLevs, atmdat.nChiantiMaxLevels, MAX_NUM_LEVELS );
743  }
744  }
745  else
746  {
747  if( tolower(dBaseSpecies[intNS].chLabel[0]) == 'f' && tolower(dBaseSpecies[intNS].chLabel[1]) == 'e')
748  {
749  // Fe is special case with more levels
750  nMolLevs = MIN3(nElvlcLines, atmdat.nChiantiMaxLevelsFe,MAX_NUM_LEVELS );
751  }
752  else
753  {
754  nMolLevs = MIN3(nElvlcLines, atmdat.nChiantiMaxLevels,MAX_NUM_LEVELS );
755  }
756  }
757 
758  if( nMolLevs <= 0 )
759  {
760  fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
761  fprintf( ioQQQ, "The file must be corrupted.\n" );
762  fclose( ioProtCollData );
763  cdEXIT( EXIT_FAILURE );
764  }
765 
766  dBaseSpecies[intNS].numLevels_max = nMolLevs;
767  dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
768 
769  if( atmdat.lgChiantiPrint == true)
770  {
771  if( atmdat.lgChiantiExp )
772  {
773  fprintf( ioQQQ," Using CHIANTI model %s with %li experimental energy levels of %li available.\n",
774  dBaseSpecies[intNS].chLabel , nMolLevs , nMolExpLevs );
775  }
776  else
777  {
778  fprintf( ioQQQ," Using CHIANTI model %s with %li theoretical energy levels of %li available.\n",
779  dBaseSpecies[intNS].chLabel , nMolLevs , nElvlcLines );
780  }
781  }
782 
783  /*Malloc the energy array to check if the energy levels are in order*/
784  dBaseStatesEnergy = (realnum*)MALLOC((unsigned long)(nMolLevs)*sizeof(realnum));
785 
786  /*malloc the States array*/
787  dBaseStates[intNS].resize(nMolLevs);
788 
789  /* allocate the Transition array*/
790  ipdBaseTrans[intNS].reserve(nMolLevs);
791  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
792  ipdBaseTrans[intNS].reserve(ipHi,ipHi);
793  ipdBaseTrans[intNS].alloc();
794  dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
795  dBaseTrans[intNS].states() = &dBaseStates[intNS];
796  dBaseTrans[intNS].chLabel() = "Chianti";
797 
798  int ndBase = 0;
799  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
800  {
801  for( long ipLo = 0; ipLo < ipHi; ipLo++)
802  {
803  ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
804  dBaseTrans[intNS][ndBase].Junk();
805  dBaseTrans[intNS][ndBase].setLo(ipLo);
806  dBaseTrans[intNS][ndBase].setHi(ipHi);
807  ++ndBase;
808  }
809  }
810 
811  /*Reading in the energy and checking that they are in order*/
812  // C++ io works in terms of cursor movement rather than position on line
813  //# of columns to skip over the rydberg energy, we don't use it
814  const int lvl_skip_ryd = 15;
815 
816  /*Keep track of which levels have experimental data and then create a vector
817  which relates their indices to the default chianti energy indices.
818  */
819  long ncounter = 0;
820 
821  if( atmdat.lgChiantiExp )
822  {
823  //Relate Chianti level indices to a set that only include experimental levels
824  intExperIndex = (long int*)MALLOC((unsigned long)(nElvlcLines)* sizeof(long int));
825 
826  //Fill relational vector with bad values
827  for(int i = 0; i < nElvlcLines; i++)
828  {
829  intExperIndex[i] = -1;
830  }
831  }
832 
833  //Read in nrg levels to see if they are in order
834  for( long ipLev=0; ipLev<nElvlcLines; ipLev++ )
835  {
836  if(elvlcstream.is_open())
837  {
838  char thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
839  elvlcstream.seekg(lvl_skipto_nrg,ios::cur);
840  elvlcstream.get(thtemp,lvl_nrg_col);
841  fenergy = (realnum) atof(thtemp);
842 
843  if( atmdat.lgChiantiExp )
844  {
845  /* Go through the entire level list selectively choosing only experimental level energies.
846  * Store them, not zeroes, in order using ncounter to count the index.
847  * Any row on the level list where there is no experimental energy, put a -1 in the relational vector.
848  * If it is a valid experimental energy level store the new ncounter index.
849  */
850 
851  //Once we collect enough experimental levels stop looking for more.
852  if( ncounter == nMolLevs)
853  break;
854  ASSERT( ncounter < nMolLevs );
855 
856  if(fenergy != 0. || ipLev == 0 )
857  {
858  dBaseStatesEnergy[ncounter] = fenergy;
859  intExperIndex[ipLev] = ncounter;
860  ncounter++;
861  }
862  else
863  {
864  intExperIndex[ipLev] = -1;
865  }
866 
867  }
868  else
869  {
870  //Stop looking for levels when the array is full
871  if( ipLev == nMolLevs )
872  break;
873 
874  elvlcstream.seekg(lvl_skip_ryd,ios::cur);
875  elvlcstream.get(obtemp,lvl_nrg_col);
876  if(atof(obtemp) != 0.)
877  {
878  fenergy = (realnum) atof(obtemp);
879  }
880  //else
881  //fprintf( ioQQQ," WARNING: Switched to theoretical energies for species %s, level %3li\n", dBaseSpecies[intNS].chLabel, ipLev );
882  dBaseStatesEnergy[ipLev] = fenergy;
883  }
884 
885  elvlcstream.ignore(INT_MAX,'\n');
886 
887  }
888  else
889  {
890  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
891  fclose( ioProtCollData );
892  cdEXIT(EXIT_FAILURE);
893  }
894  }
895 
896  for( long ipLev=1; ipLev<nMolLevs; ipLev++ )
897  {
898  /*To check if energy levels are in order*/
899  /*If the energy levels are not in order a flag is set*/
900  if (dBaseStatesEnergy[ipLev] < dBaseStatesEnergy[ipLev-1])
901  {
902  lgEneLevOrder = false;
903  }
904  }
905 
906  // malloc vector for new energy-sorted indices
907  intNewIndex = (long int*)MALLOC((unsigned long)(nElvlcLines)* sizeof(long int));
908 
909  if(lgEneLevOrder == false)
910  {
911  /*If the energy levels are not in order and the flag is set(lgEneLevOrder=FALSE)
912  then we sort the energy levels to find the relation matrix between the old and new indices*/
913  intOldIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
914  /*We now do a modified quick sort*/
915  spsort(dBaseStatesEnergy,nMolLevs,intOldIndex,2,&interror);
916  /*intNewIndex has the key to map*/
917  for( long i=0; i<nMolLevs; i++ )
918  {
919  ASSERT( intOldIndex[i] < nMolLevs );
920  intNewIndex[intOldIndex[i]] = i;
921  }
922 
923  if( nMolLevs < nElvlcLines )
924  {
925  /* If any chianti energy levels do not have experimental values
926  * the vector that relates experimental levels to the default will not be the same size
927  * as the stored experimental energy levels.
928  * Therefore this code is required to pad the rest of the sorted vector.
929  * Any index larger than nMolLevs set to -1 */
930  for( long i=nMolLevs; i<nElvlcLines; i++)
931  {
932  intNewIndex[i] = -1;
933  }
934  }
935 
936  if(DEBUGSTATE)
937  {
938  for( long i=0; i<nMolLevs; i++ )
939  {
940  printf("The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
941  }
942  for( long i=0; i<nMolLevs; i++ )
943  {
944  printf("The %ld value of the sorted energy array is %f \n",i,dBaseStatesEnergy[i]);
945  }
946  }
947  free( intOldIndex );
948  }
949  else
950  {
951  /* if energies already in correct order, new index is same as original. */
952  for( long i=0; i<nMolLevs; i++ )
953  {
954  intNewIndex[i] = i;
955  }
956 
957  if( nMolLevs < nElvlcLines )
958  {
959  /*Pad the sorted vector so that it has the correct number of elements
960  * Any index larger than nMolLevs set to -1.*/
961  for( long i=nMolLevs; i<nElvlcLines; i++)
962  {
963  intNewIndex[i] = -1;
964  }
965  }
966  }
967 
968  /* insist that there the intNewIndex values are unique */
969  for( long i=0; i<nMolLevs; i++ )
970  {
971  for( long j=i+1; j<nMolLevs; j++ )
972  {
973  ASSERT( intNewIndex[i] != intNewIndex[j] );
974  }
975  }
976 
977  /*This will print out the relational matrix
978  * which includes the original index,the new experimental index,
979  * and the sorted experimental index.
980  * The numbers NOT be on the C scale, so 1 is the lowest level. */
981  if( DEBUGSTATE && atmdat.lgChiantiExp )
982  {
983  printf("\n\n%s Energy Level Matrix\n",dBaseSpecies[intNS].chLabel);
984  printf("(Original, Experimental, Sorted)\n");
985  for( long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
986  {
987  if( intExperIndex[ipLevOld] != -1)
988  printf("%li\t%li\t%li\n",ipLevOld+1,intExperIndex[ipLevOld]+1,intNewIndex[intExperIndex[ipLevOld]]+1);
989  }
990  printf("\n");
991  }
992 
993  //After we read in the energies we rewind the energy levels file
994  elvlcstream.seekg(0,ios::beg);
995 
996  //lvl_skipto_statwt is the # of columns to skip to statwt from left
997  const int lvl_skipto_statwt = 37;
998  //lvl_statwt_col is the length (+1) of statwt
999  const int lvl_statwt_col = 4;
1000  //Read in stat weight and energy
1001  for( long ipLevOld=0; ipLevOld<nElvlcLines; ipLevOld++ )
1002  {
1003  long ipLevNew = 0;
1004  if( atmdat.lgChiantiExp )
1005  {
1006  /* We know that non-experimental levels are stored as -1
1007  * in the observed/experimental vector.
1008  * Use that to skip over those values. */
1009 
1010  if( intExperIndex[ipLevOld] == -1 )
1011  {
1012  elvlcstream.ignore(INT_MAX,'\n');
1013  continue;
1014  }
1015  else
1016  {
1017  /* Store values based on the sorted experimental indices */
1018  ipLevNew = intNewIndex[intExperIndex[ipLevOld]];
1019  }
1020  }
1021  else
1022  {
1023  /* With level trimming on it is possible that there can be rows that
1024  * have to be skipped when using theoretical
1025  * since the levels no longer exist */
1026  if( intNewIndex[ipLevOld] == -1 )
1027  {
1028  elvlcstream.ignore(INT_MAX,'\n');
1029  continue;
1030  }
1031  else
1032  {
1033  ipLevNew = intNewIndex[ipLevOld];
1034  }
1035  }
1036 
1037  char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
1038 
1039  /*information needed for label*/
1040  strcpy(dBaseStates[intNS][ipLevNew].chLabel(), " ");
1041  strncpy(dBaseStates[intNS][ipLevNew].chLabel(),dBaseSpecies[intNS].chLabel, 4);
1042  dBaseStates[intNS][ipLevNew].chLabel()[4] = '\0';
1043 
1044  //Move cursor to and read statwt
1045  elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
1046  elvlcstream.get(gtemp,lvl_statwt_col);
1047  fstatwt = (realnum)atof(gtemp);
1048 
1049  if(fstatwt <= 0.)
1050  {
1051  fprintf( ioQQQ, " WARNING: A positive non zero value is expected for the "
1052  "statistical weight but was not found in %s"
1053  " level %li\n", chEnFilename,ipLevNew);
1054  fstatwt = 1.f;
1055  //cdEXIT(EXIT_FAILURE);
1056  }
1057  dBaseStates[intNS][ipLevNew].g() = fstatwt;
1058 
1059  //Read nrg
1060  elvlcstream.get(thtemp,lvl_nrg_col);
1061 
1062  fenergy = (realnum) atof(thtemp);
1063 
1064  /* If we are looking for theoretical energies
1065  * move over a couple columns before reading in the energy*/
1066  if( !atmdat.lgChiantiExp )
1067  {
1068  elvlcstream.seekg(lvl_skip_ryd,ios::cur);
1069  elvlcstream.get(obtemp,lvl_nrg_col);
1070  if(atof(obtemp) != 0.)
1071  {
1072  fenergy = (realnum) atof(obtemp);
1073  }
1074  }
1075  elvlcstream.ignore(INT_MAX,'\n');
1076  dBaseStates[intNS][ipLevNew].energy().set(fenergy,"cm^-1");
1077  }
1078  //Close the level stream
1079  elvlcstream.close();
1080 
1081  // highest energy transition in chianti
1082  dBaseSpecies[intNS].maxWN = 0.;
1083  /* fill in all transition energies, can later overwrite for specific radiative transitions */
1084  for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
1085  tr!= dBaseTrans[intNS].end(); ++tr)
1086  {
1087  int ipHi = (*tr).ipHi();
1088  int ipLo = (*tr).ipLo();
1089  fenergyWN = (realnum)MAX2( ENERGY_MIN_WN , dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN() );
1090 
1091  (*tr).EnergyWN() = fenergyWN;
1092  (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1093  dBaseSpecies[intNS].maxWN = MAX2(dBaseSpecies[intNS].maxWN,fenergyWN);
1094  }
1095 
1096  if(DEBUGSTATE)
1097  {
1098  /* Print out the stored data for each level.
1099  * Level index is NOT on C scale. */
1100  printf("\n%s Stored Energy Levels\n",dBaseSpecies[intNS].chLabel);
1101  printf("(Index,Energy in WN, Statistical Weight)\n");
1102  for( long ipLo = 0; ipLo < nMolLevs; ipLo++ )
1103  {
1104  printf("%li\t%e\t%.1f\n",ipLo+1,dBaseStates[intNS][ipLo].energy().WN(),dBaseStates[intNS][ipLo].g());
1105  }
1106  }
1107 
1108  /************************************************************************/
1109  /*** Read in the level numbers, Einstein A and transition wavelength ***/
1110  /************************************************************************/
1111 
1112  //Count the number of rows first
1113  long wgfarows = -1;
1114  if (wgfastream.is_open())
1115  {
1116  int nj = 0;
1117  char otemp[eof_col];
1118  while(nj != -1)
1119  {
1120  wgfastream.get(otemp,eof_col);
1121  wgfastream.ignore(INT_MAX,'\n');
1122  nj = atoi(otemp);
1123  wgfarows++;
1124  }
1125  wgfastream.seekg(0,ios::beg);
1126  }
1127  else
1128  fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
1129 
1130  //line_index_col is the length(+1) of the level indexes in the WGFA file
1131  const int line_index_col = 6;
1132  //line_nrg_to_eina is the # of columns to skip from wavelength to eina in WGFA file
1133  const int line_nrg_to_eina =15;
1134  //line_eina_col is the length(+1) of einsteinA in WGFA
1135  const int line_eina_col = 16;
1136  char lvltemp[line_index_col];
1137  //Start reading WGFA file
1138  if (wgfastream.is_open())
1139  {
1140  for (long i = 0;i<wgfarows;i++)
1141  {
1142  wgfastream.get(lvltemp,line_index_col);
1143  long ipLo = atoi(lvltemp)-1;
1144  wgfastream.get(lvltemp,line_index_col);
1145  long ipHi = atoi(lvltemp)-1;
1146 
1147  if( atmdat.lgChiantiExp )
1148  {
1149  /* If either upper or lower index is -1 in the relational vector,
1150  * skip that line in the wgfa file.
1151  * Otherwise translate the level indices.*/
1152  if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1153  {
1154  wgfastream.ignore(INT_MAX,'\n');
1155  continue;
1156  }
1157  else
1158  {
1159  ipLo = intNewIndex[intExperIndex[ipLo]];
1160  ipHi = intNewIndex[intExperIndex[ipHi]];
1161  }
1162  }
1163  else
1164  {
1165  /* With level trimming on it is possible that there can be rows that
1166  * have to be skipped when using theoretical
1167  * since the levels no longer exist */
1168  if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
1169  {
1170  wgfastream.ignore(INT_MAX,'\n');
1171  continue;
1172  }
1173  else
1174  {
1175  ipLo = intNewIndex[ipLo];
1176  ipHi = intNewIndex[ipHi];
1177  }
1178  }
1179 
1180  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1181  {
1182  // skip these lines
1183  wgfastream.ignore(INT_MAX,'\n');
1184  continue;
1185  }
1186 
1187  if( ipHi == ipLo )
1188  {
1189  fprintf( ioQQQ," WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
1190  wgfastream.ignore(INT_MAX,'\n');
1191  continue;
1192  }
1193 
1194 
1195  ASSERT( ipHi != ipLo );
1196  ASSERT(ipHi >= 0);
1197  ASSERT(ipLo >= 0);
1198 
1199  // sometimes the CHIANTI datafiles list the highest index first as in the middle of these five lines in ne_10.wgfa:
1200  // ...
1201  // 8 10 187.5299 0.000e+00 4.127e+05 3d 2D1.5 - 4s 2S0.5 E2
1202  // 9 10 187.6573 0.000e+00 6.197e+05 3d 2D2.5 - 4s 2S0.5 E2
1203  // 11 10 4842624.0000 1.499e-05 9.423e-06 4p 2P0.5 - 4s 2S0.5 E1
1204  // 1 11 9.7085 1.892e-02 6.695e+11 1s 2S0.5 - 4p 2P0.5 E1
1205  // 2 11 48.5157 6.787e-02 9.618e+10 2s 2S0.5 - 4p 2P0.5 E1
1206  // ...
1207  // so, just set ipHi (ipLo) equal to the max (min) of the two indices.
1208  // NB NB NB it looks like this may depend upon whether one uses observed or theoretical energies.
1209 
1210  //Read in wavelengths
1211  char trantemp[lvl_nrg_col];
1212  wgfastream.get(trantemp,lvl_nrg_col);
1213  fWLAng = (realnum)atof(trantemp);
1214 
1215  /* \todo 2 CHIANTI labels the H 1 2-photon transition as z wavelength of zero.
1216  * Should we just ignore all of the wavelengths in this file and use the
1217  * difference of level energies instead. */
1218 
1219  if( atmdat.lgChiantiExp )
1220  {
1221  /* Sometimes the wgfa file lists the levels to where ipLo > ipHi.
1222  * This seems to be related to the theoretical energies being out of order for those transitions.
1223  * When this is true it seems that the associated wavelength is negative which means they are using theoretical energies,
1224  * even though they have experimental energies.
1225  * For now, print that the levels are switched and ignore the lines if the wavelength is negative. */
1226 
1227  if( ipHi < ipLo )
1228  {
1229  if( strcmp(dBaseSpecies[intNS].chLabel,"Fe 3") == 0)
1230  {
1231  long ipLoTemp = ipLo;
1232  long ipHiTemp = ipHi;
1233  ipHi = max( ipHiTemp, ipLoTemp );
1234  ipLo = min( ipHiTemp, ipLoTemp );
1235  if( atmdat.lgChiantiPrint )
1236  {
1237  fprintf( ioQQQ," WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
1238  dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
1239  }
1240  }
1241  else if( atmdat.lgChiantiPrint )
1242  {
1243  fprintf( ioQQQ," WARNING: Upper and Lower Indices are reversed.Species: %6s, Indices: %3li %3li Wavelength: %e \n",
1244  dBaseSpecies[intNS].chLabel, ipLo+1, ipHi+1,fWLAng );
1245  }
1246  }
1247 
1248  if( fWLAng < 0.)
1249  {
1250  // skip these lines
1251  wgfastream.ignore(INT_MAX,'\n');
1252  if( atmdat.lgChiantiPrint )
1253  {
1254  fprintf(ioQQQ,"WARNING: Skipping Transition %li to %li of %s.\n",ipLo+1,ipHi+1,dBaseSpecies[intNS].chLabel);
1255  }
1256  continue;
1257  }
1258 
1259  }
1260  else
1261  {
1262  /* If Cloudy is using theoretical energies, then just flip the levels. */
1263  if( ipHi < ipLo )
1264  {
1265  long ipLoTemp = ipLo;
1266  long ipHiTemp = ipHi;
1267  ipHi = max( ipHiTemp, ipLoTemp );
1268  ipLo = min( ipHiTemp, ipLoTemp );
1269  fprintf( ioQQQ," WARNING: Swapped level indices for species %6s, indices %3li %3li with wavelength %e \n",
1270  dBaseSpecies[intNS].chLabel, ipLoTemp, ipHiTemp,fWLAng );
1271  }
1272  }
1273 
1274  /* If the given wavelength is negative, then theoretical enegies are being used.
1275  * Take the difference in stored theoretical energies.
1276  * It should equal the absolute value of the wavelength in the wgfa file. */
1277  if( !atmdat.lgChiantiExp && fWLAng <= 0. )
1278  {
1279  //if( fWLAng < 0.)
1280  //fprintf( ioQQQ," WARNING: Negative wavelength for species %6s, indices %3li %3li \n", dBaseSpecies[intNS].chLabel, ipLo, ipHi);
1281  fWLAng = (realnum)(1e8/
1282  (dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN()));
1283  }
1284  //Skip from end of Wavelength to Einstein A and read in
1285  wgfastream.seekg(line_nrg_to_eina,ios::cur);
1286  wgfastream.get(trantemp,line_eina_col);
1287  feinsteina = (realnum)atof(trantemp);
1288  if( feinsteina < 1e-20 )
1289  {
1290  static bool notPrintedYet = true;
1291  if( notPrintedYet && atmdat.lgChiantiPrint)
1292  {
1293  fprintf( ioQQQ," WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
1294  notPrintedYet = false;
1295  }
1296  wgfastream.ignore(INT_MAX,'\n');
1297  continue;
1298  }
1299 
1300  fixit(); // may need to do something with these rates
1301  //Read in the rest of the line and look for auto
1302  wgfastream.getline(chLine,INT_MAX);
1303  TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
1304  if( nMatch("auto", chLine) )
1305  {
1306  if( (*tr).hasEmis() )
1307  {
1308  (*tr).Emis().AutoIonizFrac() =
1309  feinsteina/((*tr).Emis().Aul() + feinsteina);
1310  ASSERT( (*tr).Emis().AutoIonizFrac() >= 0. );
1311  ASSERT( (*tr).Emis().AutoIonizFrac() <= 1. );
1312  }
1313  continue;
1314  }
1315 
1316  if( (*tr).hasEmis() )
1317  {
1318  fprintf(ioQQQ," PROBLEM duplicate transition found by atmdat_chianti in %s, "
1319  "wavelength=%f\n", chTraFilename,fWLAng);
1320  fclose( ioProtCollData );
1321  cdEXIT(EXIT_FAILURE);
1322  }
1323  (*tr).AddLine2Stack();
1324  (*tr).Emis().Aul() = feinsteina;
1325  (*tr).WLAng() = fWLAng;
1326 
1327  fenergyWN = (realnum)(1e+8/fWLAng);
1328 
1329  // TODO::Check the wavelength in the file with the difference in energy levels
1330 
1331  (*tr).EnergyWN() = fenergyWN;
1332  (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
1333  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(), (*tr).EnergyWN(), (*(*tr).Hi()).g());
1334 
1335  if(DEBUGSTATE)
1336  {
1337  fprintf( ioQQQ, "The lower level is %ld \n",ipLo);
1338  fprintf( ioQQQ, "The upper level is %ld \n",ipHi);
1339  fprintf( ioQQQ, "The Einstein A is %f \n",(*tr).Emis().Aul());
1340  fprintf( ioQQQ, "The wavelength of the transition is %f \n",(*tr).WLAng() );
1341  }
1342  }
1343  }
1344  else fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
1345  wgfastream.close();
1346 
1347  /* Malloc space for splines */
1348  AtmolCollSplines[intNS] = (CollSplinesArray***)MALLOC(nMolLevs *sizeof(CollSplinesArray**));
1349  for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
1350  {
1351  AtmolCollSplines[intNS][ipHi] = (CollSplinesArray **)MALLOC((unsigned long)(nMolLevs)*sizeof(CollSplinesArray *));
1352  for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
1353  {
1354  AtmolCollSplines[intNS][ipHi][ipLo] =
1355  (CollSplinesArray *)MALLOC((unsigned long)(ipNCOLLIDER)*sizeof(CollSplinesArray ));
1356 
1357  for( long k=0; k<ipNCOLLIDER; k++ )
1358  {
1359  /* initialize all spline variables */
1360  AtmolCollSplines[intNS][ipHi][ipLo][k].collspline = NULL;
1361  AtmolCollSplines[intNS][ipHi][ipLo][k].SplineSecDer = NULL;
1362  AtmolCollSplines[intNS][ipHi][ipLo][k].nSplinePts = -1;
1363  AtmolCollSplines[intNS][ipHi][ipLo][k].intTranType = -1;
1364  AtmolCollSplines[intNS][ipHi][ipLo][k].EnergyDiff = BIGDOUBLE;
1365  AtmolCollSplines[intNS][ipHi][ipLo][k].ScalingParam = BIGDOUBLE;
1366  }
1367  }
1368  }
1369 
1370  /************************************/
1371  /*** Read in the collisional data ***/
1372  /************************************/
1373 
1374  // ipCollider 0 is electrons
1375  for( long ipCollider=0; ipCollider < 1; ipCollider++ )
1376  {
1377  char chFilename[FILENAME_PATH_LENGTH_2];
1378 
1379  if( ipCollider == ipELECTRON )
1380  {
1381  strcpy( chFilename, chEleColFilename );
1382  }
1383  else if( ipCollider == ipPROTON )
1384  {
1385  if( !lgProtonData )
1386  break;
1387  fprintf( ioQQQ," WARNING: Chianti proton collision data have different format than electron data files!\n" );
1388  strcpy( chFilename, chProColFilename );
1389  }
1390  else
1391  TotalInsanity();
1392 
1393  /*Dummy string used for convenience*/
1394  strcpy(chLine,"A");
1395 
1396  fstream splupsstream;
1397  open_data(splupsstream, chFilename,mode_r);
1398 
1399  //cs_eof_col is the length(+1) of the first column used for finding the end of file
1400  const int cs_eof_col = 4;
1401  //cs_index_col is the length(+1) of the indexes in the CS file
1402  const int cs_index_col = 4;
1403  //cs_trantype_col is the length(+1) of the transition type in the CS file
1404  const int cs_trantype_col = 4;
1405  //cs_values_col is the length(+1) of the other values in the CS file
1406  //including: GF, nrg diff, scaling parameter, and spline points
1407  const int cs_values_col = 11;
1408  //Determine the number of rows in the CS file
1409  if (splupsstream.is_open())
1410  {
1411  int nj = 0;
1412  //splupslines is -1 since the loop runs 1 extra time
1413  long splupslines = -1;
1414  char otemp[cs_eof_col];
1415  while(nj != -1)
1416  {
1417  splupsstream.get(otemp,cs_eof_col);
1418  splupsstream.ignore(INT_MAX,'\n');
1419  nj = atoi(otemp);
1420  splupslines++;
1421  }
1422  splupsstream.seekg(0,ios::beg);
1423 
1424  for (int m = 0;m<splupslines;m++)
1425  {
1426  if( ipCollider == ipELECTRON )
1427  {
1428  splupsstream.seekg(6,ios::cur);
1429  }
1430 
1431  /* level indices */
1432  splupsstream.get(otemp,cs_index_col);
1433  long ipLo = atoi(otemp)-1;
1434  splupsstream.get(otemp,cs_index_col);
1435  long ipHi = atoi(otemp)-1;
1436 
1437  /* If either upper or lower index is -1 in the relational vector,
1438  * skip that line in the splups file.
1439  * Otherwise translate the level indices.*/
1440  if( atmdat.lgChiantiExp )
1441  {
1442  if( intExperIndex[ipLo] == - 1 || intExperIndex[ipHi] == -1 )
1443  {
1444  splupsstream.ignore(INT_MAX,'\n');
1445  continue;
1446  }
1447  else
1448  {
1449  ipLo = intNewIndex[intExperIndex[ipLo]];
1450  ipHi = intNewIndex[intExperIndex[ipHi]];
1451  }
1452 
1453  }
1454  else
1455  {
1456  /* With level trimming on it is possible that there can be rows that
1457  * have to be skipped when using theoretical
1458  * since the levels no longer exist */
1459  if( intNewIndex[ipLo] == - 1 || intNewIndex[ipHi] == -1 )
1460  {
1461  splupsstream.ignore(INT_MAX,'\n');
1462  continue;
1463  }
1464  else
1465  {
1466  ipLo = intNewIndex[ipLo];
1467  ipHi = intNewIndex[ipHi];
1468  }
1469  }
1470 
1471  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
1472  {
1473  // skip these transitions
1474  splupsstream.ignore(INT_MAX,'\n');
1475  continue;
1476  }
1477 
1478  ASSERT( ipLo < nMolLevs );
1479  ASSERT( ipHi < nMolLevs );
1480  /*Transition Type*/
1481  splupsstream.get(otemp,cs_trantype_col);
1482  intTranType = atoi(otemp);
1483  char qtemp[cs_values_col];
1484  splupsstream.get(qtemp,cs_values_col);
1485  /*Energy Difference*/
1486  splupsstream.get(qtemp,cs_values_col);
1487  fEnergyDiff = atof(qtemp);
1488  /*Scaling Parameter*/
1489  splupsstream.get(qtemp,cs_values_col);
1490  fScalingParam = atof(qtemp);
1491 
1492  ASSERT( ipLo < nMolLevs );
1493  ASSERT( ipHi < nMolLevs );
1494  ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline == NULL );
1495  ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer == NULL );
1496 
1497  const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
1498  STATIC_ASSERT(CHIANTI_SPLINE_MAX > CHIANTI_SPLINE_MIN);
1499 
1500  /*We malloc the space here*/
1501  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline =
1502  (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
1503  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer =
1504  (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
1505 
1506  /* always read at least CHIANTI_SPLINE_MIN */
1507  for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
1508  {
1509  //Look at the next character to see if it is the end of line.
1510  char p = splupsstream.peek();
1511  if( p == '\n' )
1512  {
1513  /* there should be 5 or 9 spline points. If we got EOL,
1514  * insist that we were trying to read the 6th or 10th. */
1515  if( (intsplinepts != CHIANTI_SPLINE_MAX) && (intsplinepts != CHIANTI_SPLINE_MIN) )
1516  {
1517  static bool notPrintedYet = true;
1518  if( notPrintedYet )
1519  {
1520  fprintf( ioQQQ, " WARNING: Wrong number of spline points in %s.\n", chFilename);
1521  notPrintedYet = false;
1522  }
1523  for( long j=intsplinepts-1; j<CHIANTI_SPLINE_MAX; j++ )
1524  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[j] = 0.;
1525  }
1526  else
1527  break;
1528  }
1529  else
1530  {
1531  if( intsplinepts >= CHIANTI_SPLINE_MAX )
1532  {
1533  fprintf( ioQQQ, " WARNING: More spline points than expected in %s, indices %3li %3li. Ignoring extras.\n", chFilename, ipHi, ipLo );
1534  break;
1535  }
1536  ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
1537  double temp;
1538  //Store a single spline point then look for more
1539  splupsstream.get(qtemp,cs_values_col);
1540  temp = atof(qtemp);
1541  if(temp < 0)
1542  {
1543  temp = 0.;
1544  }
1545  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
1546  }
1547 
1548  }
1549 
1550  ASSERT( (intsplinepts == CHIANTI_SPLINE_MAX) || (intsplinepts == CHIANTI_SPLINE_MIN) );
1551 
1552  /*The zeroth element contains the number of spline points*/
1553  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].nSplinePts = intsplinepts;
1554  /*Transition type*/
1555  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].intTranType = intTranType;
1556  /*Energy difference between two levels in Rydbergs*/
1557  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].EnergyDiff = fEnergyDiff;
1558  /*Scaling parameter C*/
1559  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;
1560 
1561  /*Once the spline points have been filled,fill the second derivatives structure*/
1562  /*Creating spline points array*/
1563  xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
1564  spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
1565  spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
1566 
1567  // should be able to just loop to intsplinepts -- ASSERT above protects
1568  if(intsplinepts == CHIANTI_SPLINE_MIN)
1569  {
1570  for(intxs=0;intxs<CHIANTI_SPLINE_MIN;intxs++)
1571  {
1572  xs[intxs] = 0.25*intxs;
1573  spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
1574  }
1575  }
1576  else if(intsplinepts == CHIANTI_SPLINE_MAX)
1577  {
1578  for(intxs=0;intxs<CHIANTI_SPLINE_MAX;intxs++)
1579  {
1580  xs[intxs] = 0.125*intxs;
1581  spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
1582  }
1583  }
1584  else
1585  TotalInsanity();
1586 
1587  spline(xs, spl,intsplinepts,2e31,2e31,spl2);
1588 
1589  /*Filling the second derivative structure*/
1590  for( long i=0; i<intsplinepts; i++)
1591  {
1592  AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer[i] = spl2[i];
1593  }
1594 
1595  free( xs );
1596  free( spl );
1597  free( spl2 );
1598  splupsstream.ignore(INT_MAX,'\n');
1599  }
1600  splupsstream.close();
1601  }
1602  }
1603 
1604  // free some memory
1605  free( dBaseStatesEnergy );
1606  free( intNewIndex );
1607  if( atmdat.lgChiantiExp)
1608  {
1609  free( intExperIndex );
1610  }
1611 
1612  // close open file handles
1613  fclose( ioElecCollData );
1614  if( lgProtonData )
1615  fclose( ioProtCollData );
1616 
1617  return;
1618 }
AtmolCollSplines
CollSplinesArray **** AtmolCollSplines
Definition: taulines.cpp:19
t_atmdat::lgChiantiExp
bool lgChiantiExp
Definition: atmdat.h:237
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
t_CollSplinesArray::intTranType
long intTranType
Definition: cddefines.h:1286
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:17
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
realnum
float realnum
Definition: cddefines.h:103
rfield.h
t_atmdat::nChiantiMaxLevelsFe
long nChiantiMaxLevelsFe
Definition: atmdat.h:212
spline
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
Definition: thirdparty.h:117
MIN3
#define MIN3(a, b, c)
Definition: cddefines.h:766
ENERGY_MIN_WN
const double ENERGY_MIN_WN
Definition: atmdat_chianti.cpp:15
spsort
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1100
thirdparty.h
t_StoutColls::lgIsRate
bool lgIsRate
Definition: cddefines.h:1303
trace.h
GetGF
double GetGF(double trans_prob, double enercm, double gup)
Definition: lines_service.cpp:101
DEBUGSTATE
static const bool DEBUGSTATE
Definition: atmdat_chianti.cpp:13
ProxyIterator
Definition: proxy_iterator.h:58
lines_service.h
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
atmdat.h
t_StoutColls::temps
double * temps
Definition: cddefines.h:1299
ipPROTON
@ ipPROTON
Definition: collision.h:10
t_StoutColls::collstrs
double * collstrs
Definition: cddefines.h:1301
mode_r
const ios_base::openmode mode_r
Definition: cpu.h:212
t_species::maxWN
double maxWN
Definition: cddefines.h:1257
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_CollSplinesArray
Definition: cddefines.h:1275
t_species::numLevels_local
long numLevels_local
Definition: cddefines.h:1241
t_CollSplinesArray::collspline
double * collspline
Definition: cddefines.h:1282
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
uncaps
void uncaps(char *chCard)
Definition: service.cpp:263
t_atmdat::lgChiantiPrint
bool lgChiantiPrint
Definition: atmdat.h:235
nMatch
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:451
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
t_CollSplinesArray::ScalingParam
double ScalingParam
Definition: cddefines.h:1288
RefIndex
double RefIndex(double EnergyWN)
Definition: lines_service.cpp:141
atmdat_CHIANTI_readin
void atmdat_CHIANTI_readin(long intNS, char *chPrefix)
Definition: atmdat_chianti.cpp:595
ipNCOLLIDER
@ ipNCOLLIDER
Definition: collision.h:18
MAX2
#define MAX2
Definition: cddefines.h:782
t_species::lgLTE
bool lgLTE
Definition: cddefines.h:1259
ipELECTRON
@ ipELECTRON
Definition: collision.h:9
t_StoutColls
Definition: cddefines.h:1294
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
StoutCollData
StoutColls **** StoutCollData
Definition: taulines.cpp:20
t_CollSplinesArray::SplineSecDer
double * SplineSecDer
Definition: cddefines.h:1283
t_species::lgLAMDA
bool lgLAMDA
Definition: cddefines.h:1247
BIGDOUBLE
const double BIGDOUBLE
Definition: cpu.h:194
tolower
char tolower(char c)
Definition: cddefines.h:691
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
atmdat_STOUT_readin
void atmdat_STOUT_readin(long intNS, char *chPrefix)
Definition: atmdat_chianti.cpp:17
min
long min(int a, long b)
Definition: cddefines.h:723
t_atmdat::nChiantiMaxLevels
long nChiantiMaxLevels
Definition: atmdat.h:214
physconst.h
t_CollSplinesArray::nSplinePts
long nSplinePts
Definition: cddefines.h:1285
ipdBaseTrans
vector< multi_arr< int, 2 > > ipdBaseTrans
Definition: taulines.cpp:16
t_atmdat::lgStoutPrint
bool lgStoutPrint
Definition: atmdat.h:245
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
fixit
void fixit(void)
Definition: service.cpp:991
t_atmdat::nStoutMaxLevels
long nStoutMaxLevels
Definition: atmdat.h:218
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
taulines.h
STATIC_ASSERT
#define STATIC_ASSERT(x)
Definition: cddefines.h:113
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:15
t_CollSplinesArray::EnergyDiff
double EnergyDiff
Definition: cddefines.h:1287
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
t_species::numLevels_max
long numLevels_max
Definition: cddefines.h:1239
t_species::lgMolecular
bool lgMolecular
Definition: cddefines.h:1245
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
t_StoutColls::ntemps
long ntemps
Definition: cddefines.h:1297
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17