cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_lamda.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 "atmdat.h"
5 #include "rfield.h"
6 #include "taulines.h"
7 #include "trace.h"
8 #include "lines_service.h"
9 #include "parse_species.h"
10 
11 static const bool DEBUGSTATE=false;
12 
13 static void check_LAMDA_comment(const char* chLine, char* chEFilename)
14 {
15  const char* s= chLine;
16  while (isspace(*s))
17  ++s;
18  if (*s != '!')
19  {
20  fprintf (ioQQQ," PROBLEM: The LAMDA data file '%s' appears to be corrupt\n",
21  chEFilename);
22  fprintf (ioQQQ," Expected line starting '!', found\n%s\n",chLine);
24  }
25 }
26 
27 /*Separate Routine to read in the molecules*/
28 void atmdat_LAMDA_readin( long intNS, char *chEFilename )
29 {
30  DEBUG_ENTRY( "atmdat_LAMDA_readin()" );
31 
32  long int intlnct;
33 
34  FILE *ioLevData;
35  realnum fstatwt,fenergyWN,fenergy,feinsteina;
36 
37  char chLine[FILENAME_PATH_LENGTH_2];
38 
39  ASSERT( intNS >= 0 );
40 
41  static const int MAX_NUM_LEVELS = 999;
42 
43  dBaseSpecies[intNS].lgMolecular = true;
44  dBaseSpecies[intNS].lgLTE = false;
45 
46  /*Load the species name in the dBaseSpecies array structure*/
47  if(DEBUGSTATE)
48  printf("The name of the %li species is %s \n",intNS+1,dBaseSpecies[intNS].chLabel);
49 
50  /*Open the files*/
51  if( trace.lgTrace )
52  fprintf( ioQQQ," moldat_readin opening %s:",chEFilename);
53 
54  ioLevData = open_data( chEFilename, "r" );
55 
56  intlnct = 0;
57  while(intlnct < 3)
58  {
59  intlnct++;
60  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
61  {
62  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
64  }
65  }
66  /*Extracting out the molecular weight*/
67  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
68  {
69  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
71  }
72  dBaseSpecies[intNS].fmolweight = (realnum)atof(chLine);
73 
74  /*Discard this line*/
75  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
76  {
77  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
79  }
80 
81  // sections of the file are separated by line that begin with "!"
82  check_LAMDA_comment(chLine,chEFilename);
83 
84  /*Reading in the number of energy levels*/
85  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
86  {
87  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
89  }
90  long nMolLevs = (long)atoi(chLine);
91 
92  long HighestIndexInFile = nMolLevs;
93 
94  dBaseSpecies[intNS].numLevels_max = HighestIndexInFile;
95 
97 
98  /* truncate model to preset max */
99  nMolLevs = MIN3( atmdat.nLamdaMaxLevels, nMolLevs, MAX_NUM_LEVELS );
100 
101  if(nMolLevs <= 0)
102  {
103  fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEFilename );
104  fprintf( ioQQQ, "The file must be corrupted.\n" );
105  cdEXIT( EXIT_FAILURE );
106  }
107 
108  if( atmdat.lgLamdaPrint == true)
109  {
110  fprintf( ioQQQ,"Using LAMDA model %s with %li levels of %li available.\n",
111  dBaseSpecies[intNS].chLabel , nMolLevs , HighestIndexInFile );
112  }
113 
114  if (dBaseSpecies[intNS].setLevels != -1)
115  {
116  if (dBaseSpecies[intNS].setLevels > HighestIndexInFile)
117  {
118  fprintf( ioQQQ,"Using LAMDA model %s with %li requested, only %li energy levels available.\n",
119  dBaseSpecies[intNS].chLabel ,dBaseSpecies[intNS].setLevels, HighestIndexInFile );
120  nMolLevs = HighestIndexInFile;
121  }
122  else
123  {
124  nMolLevs = dBaseSpecies[intNS].setLevels;
125  }
126  }
127 
128  dBaseSpecies[intNS].numLevels_max = nMolLevs;
129  dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
130 
131  /*malloc the States array*/
132  dBaseStates[intNS].init(dBaseSpecies[intNS].chLabel,nMolLevs);
133  /* allocate the Transition array*/
134  ipdBaseTrans[intNS].reserve(nMolLevs);
135  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
136  ipdBaseTrans[intNS].reserve(ipHi,ipHi);
137  ipdBaseTrans[intNS].alloc();
138  dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
139  dBaseTrans[intNS].states() = &dBaseStates[intNS];
140  dBaseTrans[intNS].chLabel() = dBaseSpecies[intNS].chLabel;
141  dBaseSpecies[intNS].database = "LAMDA";
142 
143  int ndBase = 0;
144  for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
145  {
146  for( long ipLo = 0; ipLo < ipHi; ipLo++)
147  {
148  ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
149  dBaseTrans[intNS][ndBase].Junk();
150  dBaseTrans[intNS][ndBase].setLo(ipLo);
151  dBaseTrans[intNS][ndBase].setHi(ipHi);
152  ++ndBase;
153  }
154  }
155 
156  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
157  {
158  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
160  }
161 
162  // sections of the file are separated by line that begin with "!"
163  check_LAMDA_comment(chLine,chEFilename);
164 
165  for( long ipLev=0; ipLev<HighestIndexInFile; ipLev++)
166  {
167  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
168  {
169  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
171  }
172 
173  // skip these levels
174  if( ipLev >= nMolLevs )
175  continue;
176 
177  long i = 1;
178  bool lgEOL;
179  long index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
180  fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
181  fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
182 
183  if( index != ipLev + 1 )
184  {
185  fprintf( ioQQQ, "internal error in file %s\n", chEFilename );
187  }
188  dBaseStates[intNS][ipLev].energy().set(fenergy,"cm^-1");
189  dBaseStates[intNS][ipLev].g() = fstatwt;
190 
191  if( ipLev > 0 )
192  {
193  // Use volatile variables to ensure normalization to
194  // standard precision before comparison
195  volatile realnum enlev = dBaseStates[intNS][ipLev].energy().Ryd();
196  volatile realnum enprev = dBaseStates[intNS][ipLev-1].energy().Ryd();
197  if (enlev < enprev)
198  {
199  fprintf( ioQQQ, " The energy levels are not in order in species %s at index %li.\n",
200  dBaseSpecies[intNS].chLabel, ipLev );
202  }
203  }
204  if(DEBUGSTATE)
205  {
206  printf("The converted energy is %f \n",dBaseStates[intNS][ipLev].energy().WN());
207  printf("The value of g is %f \n",dBaseStates[intNS][ipLev].g());
208  }
209  }
210 
211  /* fill in all transition energies, can later overwrite for specific radiative transitions */
212  for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
213  tr!= dBaseTrans[intNS].end(); ++tr)
214  {
215  int ipHi = (*tr).ipHi();
216  int ipLo = (*tr).ipLo();
217  //fenergyWN = (realnum)MAX2( 1.01*rfield.emm()*RYD_INF, dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN() );
218  fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN());
219 
220  (*tr).EnergyWN() = fenergyWN;
221 
222  /* there are OH hyperfine levels where i+1 and i have exactly
223  * the same energy. The refractive index routine will FPE with
224  * an energy of zero - so we do this test */
225  if( fenergyWN>SMALLFLOAT )
226  (*tr).WLAng() = (realnum)(1e+8f/fenergyWN/RefIndex(fenergyWN));
227  else
228  (*tr).WLAng() = 1e30f;
229  }
230 
231  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
232  {
233  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
235  }
236  if(chLine[0]!='!')
237  {
238  fprintf( ioQQQ, " The number of energy levels in file %s is not correct, expected to find line starting with!.\n",chEFilename);
239  fprintf( ioQQQ , "%s\n", chLine );
241  }
242  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
243  {
244  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
246  }
247  long intgrtct = atoi(chLine);
248  /*The above gives the number of radiative transitions*/
249  if(intgrtct <= 0)
250  {
251  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
253  }
254  if(DEBUGSTATE)
255  {
256  printf("The number of radiative transitions is %li \n",intgrtct);
257  }
258  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
259  {
260  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
262  }
263 
264  /*intrtct refers to radiative transitions count*/
265  for( long intrtct=0; intrtct<intgrtct; intrtct++)
266  {
267  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
268  {
269  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
271  }
272 
273  long i = 1;
274  bool lgEOL;
275 
276  long index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
277  long ipHiInFile = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
278  long ipLoInFile = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
279  long ipHi = ipHiInFile - 1;
280  long ipLo = ipLoInFile - 1;
281 
282  if( ! ( ipHi >= 0 && ipLo >= 0 && ipHi > ipLo && index == intrtct + 1 ) )
283  {
284  fprintf( ioQQQ, "internal error in file %s\n", chEFilename );
286  }
287 
288  // skip these lines
289  if( ipLo >= nMolLevs || ipHi >= nMolLevs )
290  continue;
291 
292  feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
293  /* don't need the energy in GHz, so throw it away. */
294  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
295  fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy().WN() -dBaseStates[intNS][ipLo].energy().WN());
296  fenergyWN = MAX2( fenergyWN, 1.01 * RYD_INF * rfield.emm() );
297 
298  TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
299  ASSERT(!(*tr).hasEmis()); // Check for duplicates
300  (*tr).AddLine2Stack();
301  (*tr).Emis().Aul() = feinsteina;
302  ASSERT( !isnan( (*tr).EnergyK() ) );
303  (*tr).EnergyWN() = fenergyWN;
304  (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
305 
306  (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(),(*tr).EnergyWN(), (*(*tr).Hi()).g());
307 
308  (*tr).setComment( db_comment_tran_levels( ipLoInFile, ipHiInFile ) );
309 
310  if(DEBUGSTATE)
311  {
312  printf("The upper level is %ld \n",ipHi+1);
313  printf("The lower level is %ld \n",ipLo+1);
314  printf("The Einstein A is %E \n",(*tr).Emis().Aul());
315  printf("The Energy of the transition is %E \n",(*tr).EnergyK());
316  }
317  }
318 
319  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
320  {
321  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
323  }
324 
325  if(chLine[0]!='!')
326  {
327  fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename);
329  }
330 
331  long nCollPartners = -1;
332  /*Getting the number of collisional partners*/
333  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
334  {
335  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
337  }
338  else
339  {
340  nCollPartners = atoi(chLine);
341  }
342  /* if the number of colliders has a negative value, do the molecule in LTE */
343  if (nCollPartners < 0)
344  dBaseSpecies[intNS].lgLTE = true;
345  /*Checking the number of collisional partners does not exceed 9*/
346  if( nCollPartners > ipNCOLLIDER )
347  {
348  fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
350  }
351 
352  FunctPtr func = new FunctLAMDA();
353  // loop over partners
354  for( long ipPartner = 0; ipPartner < nCollPartners; ++ ipPartner )
355  {
356  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
357  {
358  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
360  }
361 
362  check_LAMDA_comment(chLine,chEFilename);
363 
364  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
365  {
366  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
368  }
369  /*Extract out the name of the collider*/
370  /*The following are the rules expected in the datafiles to extract the names of the collider*/
371  /*The line which displays the species and the collider starts with a number*/
372  /*This refers to the collider in the Leiden database*/
373  /*In the Leiden database 1 refers to H2,2 to para-H2,3 to ortho-H2
374  4 to electrons,5 to H and 6 to He*/
375  char *chCollName = strtok(chLine," ");
376  /*Leiden Collider Index*/
377  int intLColliderIndex = atoi(chCollName);
378  int intCollIndex = -1;
379  /*In Cloudy,We assign the following indices for the colliders*/
380  /*electron=0,proton=1,He+=2,He++=3,atomic hydrogen=4,He=5,oH2=6,pH2=7,H2=8*/
381 
382  if(intLColliderIndex == 1)
383  {
384  intCollIndex = ipH2;
385  }
386  else if(intLColliderIndex == 2)
387  {
388  intCollIndex = ipH2_PARA;
389  }
390  else if(intLColliderIndex == 3)
391  {
392  intCollIndex = ipH2_ORTHO;
393  }
394  else if(intLColliderIndex == 4)
395  {
396  intCollIndex = ipELECTRON;
397  }
398  else if(intLColliderIndex == 5)
399  {
400  intCollIndex = ipATOM_H;
401  }
402  else if(intLColliderIndex == 6)
403  {
404  intCollIndex = ipATOM_HE;
405  }
406  else
407  {
408  // this happens for some LAMDA files (as of Jan 20, 2009) because there is no integer
409  // in the datafile indicating which collider the subsequent table is for
410  TotalInsanity();
411  }
412 
413  ASSERT( intCollIndex < ipNCOLLIDER );
414 
415  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
416  {
417  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
419  }
420 
421  check_LAMDA_comment(chLine,chEFilename);
422 
423  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
424  {
425  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
427  }
428  // Get the number of transitions
429  long nCollTrans = atoi(chLine);
430 
431  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
432  {
433  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
435  }
436 
437  check_LAMDA_comment(chLine,chEFilename);
438 
439  if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
440  {
441  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
443  }
444  // Get the number of temperatures
445  long intCollTemp = atoi(chLine);
446 
447  // Read and store the table of rate coefficients
448  ReadCollisionRateTable( AtmolCollRateCoeff[intNS][intCollIndex],
449  ioLevData, func, nMolLevs, intCollTemp, nCollTrans );
450  }
451  delete func;
452 
453  fclose( ioLevData );
454 
455  return;
456 }
457 
double emm() const
Definition: mesh.h:93
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
t_atmdat atmdat
Definition: atmdat.cpp:6
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:296
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
const realnum SMALLFLOAT
Definition: cpu.h:246
long nLamdaMaxLevels
Definition: atmdat.h:397
vector< multi_arr< int, 2 > > ipdBaseTrans
Definition: taulines.cpp:17
double RefIndex(double EnergyWN)
FILE * ioQQQ
Definition: cddefines.cpp:7
multi_arr< CollRateCoeffArray, 2 > AtmolCollRateCoeff
Definition: taulines.cpp:19
static void check_LAMDA_comment(const char *chLine, char *chEFilename)
t_trace trace
Definition: trace.cpp:5
Definition: atmdat.h:461
static const bool DEBUGSTATE
bool lgLamdaPrint
Definition: atmdat.h:393
double energy(const genericState &gs)
bool lgTrace
Definition: trace.h:12
t_rfield rfield
Definition: rfield.cpp:9
void setProperties(species &sp)
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
#define ASSERT(exp)
Definition: cddefines.h:613
double GetGF(double trans_prob, double enercm, double gup)
string db_comment_tran_levels(long ipLoFile, long ipHiFile)
Definition: atmdat.h:515
void atmdat_LAMDA_readin(long intNS, char *chFileName)
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define isnan
Definition: cddefines.h:659
vector< qList > dBaseStates
Definition: taulines.cpp:16
#define MIN3(a, b, c)
Definition: cddefines.h:808
vector< species > dBaseSpecies
Definition: taulines.cpp:15
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
void ReadCollisionRateTable(CollRateCoeffArray &coll_rate_table, FILE *io, FunctPtr GetIndices, long nMolLevs, long nTemps, long nTrans)
Definition: atmdat.cpp:66
Definition: collision.h:19
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381