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