cloudy  trunk
 All Data Structures 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-2008 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 
11 #define DEBUGSTATE false
12 
13 /*Separate Routine to read in the molecules*/
14 void atmdat_lamda_readin( long intNS )
15 {
16  DEBUG_ENTRY( "atmdat_lamda_readin()" );
17 
18  int ngMolLevs = -10000,intCollIndex = -10000,intLColliderIndex = -10000;
19  /* type is set to 0 for non chianti and 1 for chianti*/
20  long int i,j,nMolLevs,intlnct,intrtct,intgrtct,intCollPart,
21  intDCollPart,intCollTran,nCollTrans,intCollTemp,intcollindex,
22  ipLo,ipHi;
23  /*intrtct refers to radiative transitions count*/
24  FILE *atmolLevDATA;
25  realnum fstatwt,fenergyK,fenergyWN,fenergy,feinsteina;
26 
27  char chLine[FILENAME_PATH_LENGTH_2] ,
28  chEFilename[FILENAME_PATH_LENGTH_2],
29  *chColltemp,*chCollName;
30 
31  ASSERT( lgSpeciesMolecule[intNS] );
32 
33  /*Load the species name in the Species array structure*/
34  if(DEBUGSTATE)
35  printf("The name of the %li species is %s \n",intNS+1,Species[intNS].chptrSpName);
36  strcpy( chEFilename , Species[intNS].chptrSpName );
37  strcat( chEFilename , ".dat");
38  uncaps( chEFilename );
39 
40  /*Open the files*/
41  if( trace.lgTrace )
42  fprintf( ioQQQ," moldat_readin opening %s:",chEFilename);
43 
44  atmolLevDATA = open_data( chEFilename, "r" );
45 
46  nMolLevs = 0;
47  ngMolLevs = 0;
48  intrtct = 0;
49  intgrtct = 0;
50  intlnct = 0;
51  ipHi = 0;
52  ipLo = 0;
53  j = 0;
54  while(intlnct < 3)
55  {
56  intlnct++;
57  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
58  {
59  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
60  cdEXIT(EXIT_FAILURE);
61  }
62  }
63  /*Extracting out the molecular weight*/
64  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
65  {
66  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
67  cdEXIT(EXIT_FAILURE);
68  }
69  Species[intNS].fmolweight = (realnum)atof(chLine);
70 
71  /*Discard this line*/
72  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
73  {
74  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
75  cdEXIT(EXIT_FAILURE);
76  }
77  /*Reading in the number of energy levels*/
78  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
79  {
80  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
81  cdEXIT(EXIT_FAILURE);
82  }
83  ngMolLevs = atoi(chLine);
84 
85  if(ngMolLevs <= 0)
86  {
87  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
88  cdEXIT(EXIT_FAILURE);
89  }
90  Species[intNS].numLevels_max = ngMolLevs;
91  Species[intNS].numLevels_local = ngMolLevs;
92  if(DEBUGSTATE)
93  {
94  printf("The number of energy levels is %li \n",Species[intNS].numLevels_max);
95  }
96  /*the above refers to the number of energy levels specified in the data file*/
97  /*malloc the States array*/
98  atmolStates[intNS] = (quantumState *)MALLOC((size_t)(ngMolLevs)*sizeof(quantumState));
99  /*malloc the Transition array*/
100  atmolTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)ngMolLevs);
101  atmolTrans[intNS][0] = NULL;
102  for( ipHi = 1; ipHi < ngMolLevs; ipHi++)
103  {
104  atmolTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi);
105  for(ipLo = 0; ipLo < ipHi; ipLo++)
106  {
107  atmolTrans[intNS][ipHi][ipLo].Lo = &atmolStates[intNS][ipLo];
108  atmolTrans[intNS][ipHi][ipLo].Hi = &atmolStates[intNS][ipHi];
109  atmolTrans[intNS][ipHi][ipLo].Emis = &DummyEmis;
110  }
111  }
112 
113  for( intcollindex = 0; intcollindex <NUM_COLLIDERS; intcollindex++ )
114  {
115  CollRatesArray[intNS][intcollindex] = NULL;
116  }
117 
118 
119  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
120  {
121  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
122  cdEXIT(EXIT_FAILURE);
123  }
124 
125  for( nMolLevs=0; nMolLevs<ngMolLevs; nMolLevs++)
126  {
127  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
128  {
129  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
130  cdEXIT(EXIT_FAILURE);
131  }
132  /*information needed for label*/
133  strcpy(atmolStates[intNS][nMolLevs].chLabel, " ");
134  strncpy(atmolStates[intNS][nMolLevs].chLabel,Species[intNS].chptrSpName, 4);
135  atmolStates[intNS][nMolLevs].chLabel[4] = '\0';
136 
137  long i = 1;
138  bool lgEOL;
139  long index;
140 
141  index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
142  fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
143  fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
144 
145  ASSERT( index == nMolLevs + 1 );
146 
147  atmolStates[intNS][nMolLevs].energy = fenergy;
148  atmolStates[intNS][nMolLevs].g = fstatwt;
149 
150  if (nMolLevs > 0)
151  {
152  if (atmolStates[intNS][nMolLevs].energy < atmolStates[intNS][nMolLevs-1].energy)
153  {
154  fprintf( ioQQQ, " The energy levels are not in order");
155  cdEXIT(EXIT_FAILURE);
156  }
157  }
158  if(DEBUGSTATE)
159  {
160  printf("The converted energy is %f \n",atmolStates[intNS][nMolLevs].energy);
161  printf("The value of g is %f \n",atmolStates[intNS][nMolLevs].g);
162  }
163  }
164 
165  /* fill in all transition energies, can later overwrite for specific radiative transitions */
166  for( i=1; i<ngMolLevs; i++ )
167  {
168  for( j=0; j<i; j++ )
169  {
170  fenergyWN = atmolStates[intNS][i].energy - atmolStates[intNS][j].energy;
171  fenergyK = (realnum)(fenergyWN*T1CM);
172 
173  atmolTrans[intNS][i][j].EnergyK = fenergyK;
174  atmolTrans[intNS][i][j].EnergyWN = fenergyWN;
175  atmolTrans[intNS][i][j].EnergyErg = (realnum)ERG1CM *fenergyWN;
176 
177  /* protect against division by zero. */
178  atmolTrans[intNS][i][j].WLAng = (realnum)(1e+8f/MAX2(1E-5f,fenergyWN)/
179  RefIndex(fenergyWN));
180  }
181  }
182 
183  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
184  {
185  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
186  cdEXIT(EXIT_FAILURE);
187  }
188  if(chLine[0]!='!')
189  {
190  fprintf( ioQQQ, " The number of energy levels in file %s is not correct.\n",chEFilename);
191  cdEXIT(EXIT_FAILURE);
192  }
193  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
194  {
195  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
196  cdEXIT(EXIT_FAILURE);
197  }
198  intgrtct = atoi(chLine);
199  /*The above gives the number of radiative transitions*/
200  if(intgrtct <= 0)
201  {
202  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
203  cdEXIT(EXIT_FAILURE);
204  }
205  if(DEBUGSTATE)
206  {
207  printf("The number of radiative transitions is %li \n",intgrtct);
208  }
209  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
210  {
211  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
212  cdEXIT(EXIT_FAILURE);
213  }
214 
215  for( intrtct=0; intrtct<intgrtct; intrtct++)
216  {
217  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
218  {
219  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
220  cdEXIT(EXIT_FAILURE);
221  }
222 
223  long i = 1;
224  bool lgEOL;
225  long index;
226 
227  index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
228  ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
229  ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
230  feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
231  /* don't need the energy in GHz, so throw it away. */
232  FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
233  fenergyK = (realnum)(((atmolStates[intNS][ipHi].energy) -(atmolStates[intNS][ipLo].energy))*T1CM);
234  ASSERT( index == intrtct + 1 );
235 
236  atmolTrans[intNS][ipHi][ipLo].Emis = AddLine2Stack(feinsteina, &atmolTrans[intNS][ipHi][ipLo]);
237  atmolTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
238  fenergyWN = (realnum)((fenergyK)/T1CM);
239  atmolTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
240  atmolTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
241  atmolTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
242 
243  ASSERT( !isnan( atmolTrans[intNS][ipHi][ipLo].EnergyK ) );
244 
245  if(DEBUGSTATE)
246  {
247  printf("The upper level is %ld \n",ipHi+1);
248  printf("The lower level is %ld \n",ipLo+1);
249  printf("The Einstein A is %E \n",atmolTrans[intNS][ipHi][ipLo].Emis->Aul);
250  printf("The Energy of the transition is %E \n",atmolTrans[intNS][ipHi][ipLo].EnergyK);
251  }
252  }
253 
254  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
255  {
256  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
257  cdEXIT(EXIT_FAILURE);
258  }
259 
260  if(chLine[0]!='!')
261  {
262  fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename);
263  cdEXIT(EXIT_FAILURE);
264  }
265 
266  /*Getting the number of collisional partners*/
267  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
268  {
269  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
270  cdEXIT(EXIT_FAILURE);
271  }
272  else
273  {
274  intCollPart = atoi(chLine);
275  }
276  /*Checking the number of collisional partners does not exceed 9*/
277  if(intCollPart > NUM_COLLIDERS-1)
278  {
279  fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
280  cdEXIT(EXIT_FAILURE);
281  }
282  /*Creating the duplicate of the number of collisional partners which is reduced*/
283  intDCollPart = intCollPart;
284  while(intDCollPart > 0)
285  {
286  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
287  {
288  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
289  cdEXIT(EXIT_FAILURE);
290  }
291 
292  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
293  {
294  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
295  cdEXIT(EXIT_FAILURE);
296  }
297  /*Extract out the name of the collider*/
298  /*The following are the rules expected in the datafiles to extract the names of the collider*/
299  /*The line which displays the species and the collider starts with a number*/
300  /*This refers to the collider in the Leiden database*/
301  /*In the Leiden database 1 referes to H2,2 to para-H2,3 to ortho-H2
302  4 to electrons,5 to H and 6 to He*/
303  chCollName = strtok(chLine," ");
304  /*Leiden Collider Index*/
305  intLColliderIndex = atoi(chCollName);
306  /*In Cloudy,We assign the following indices for the colliders*/
307  /*electron=0,proton=1,atomic hydrogen=2,He=3,He+=4,He++=5,oH2=6,pH2=7,H2=8*/
308 
309  if(intLColliderIndex == 1)
310  {
311  intCollIndex = 8;
312  }
313  else if(intLColliderIndex == 2)
314  {
315  intCollIndex = 7;
316  }
317  else if(intLColliderIndex == 3)
318  {
319  intCollIndex = 6;
320  }
321  else if(intLColliderIndex == 4)
322  {
323  intCollIndex = 0;
324  }
325  else if(intLColliderIndex == 5)
326  {
327  intCollIndex = 2;
328  }
329  else if(intLColliderIndex == 6)
330  {
331  intCollIndex = 3;
332  }
333  else
334  {
335  TotalInsanity();
336  }
337  /*This is where we allocate memory if the collider exists*/
338  /*Needed to take care of the he collisions*/
339  CollRatesArray[intNS][intCollIndex] = (double**)MALLOC((unsigned long)(ngMolLevs) * sizeof(double*));
340  for( ipHi = 0; ipHi<ngMolLevs; ipHi++ )
341  {
342  CollRatesArray[intNS][intCollIndex][ipHi] = (double*)MALLOC((unsigned long)(ngMolLevs) * sizeof(double));
343  for( ipLo = 0; ipLo<ngMolLevs; ipLo++ )
344  {
345  CollRatesArray[intNS][intCollIndex][ipHi][ipLo] = 0.0;
346  }
347  }
348  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
349  {
350  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
351  cdEXIT(EXIT_FAILURE);
352  }
353  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
354  {
355  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
356  cdEXIT(EXIT_FAILURE);
357  }
358  /*Number of coll trans*/
359  intCollTran = atoi(chLine);
360 
361  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
362  {
363  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
364  cdEXIT(EXIT_FAILURE);
365  }
366  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
367  {
368  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
369  cdEXIT(EXIT_FAILURE);
370  }
371  /*Number of coll temps*/
372  intCollTemp = atoi(chLine);
373  /*Storing the number of collisional temperatures*/
374  AtmolCollRateCoeff[intNS][intCollIndex].ntemps = intCollTemp;
375  /*Mallocing*/
376  AtmolCollRateCoeff[intNS][intCollIndex].temps =
377  (double *)MALLOC((unsigned long)intCollTemp*sizeof(double));
378  AtmolCollRateCoeff[intNS][intCollIndex].collrates =
379  (double***)MALLOC((unsigned long)(ngMolLevs)*sizeof(double**));
380  for( i=0; i<ngMolLevs; i++ )
381  {
382  AtmolCollRateCoeff[intNS][intCollIndex].collrates[i] =
383  (double **)MALLOC((unsigned long)(ngMolLevs)*sizeof(double*));
384  for( j=0; j<ngMolLevs; j++ )
385  {
386  AtmolCollRateCoeff[intNS][intCollIndex].collrates[i][j] =
387  (double *)MALLOC((unsigned long)(intCollTemp)*sizeof(double));
388  }
389  }
390 
391  /*Discard the header line*/
392  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
393  {
394  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
395  cdEXIT(EXIT_FAILURE);
396  }
397  /*Getting the collisional Temperatures*/
398  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
399  {
400  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
401  cdEXIT(EXIT_FAILURE);
402  }
403  /*Filling the collisional temps array*/
404  chColltemp =strtok(chLine," ");
405  AtmolCollRateCoeff[intNS][intCollIndex].temps[0] =(realnum) atof(chColltemp);
406  for( i=1; i<intCollTemp; i++ )
407  {
408  chColltemp =strtok(NULL," ");
409  AtmolCollRateCoeff[intNS][intCollIndex].temps[i] = atof(chColltemp);
410  }
411  /*Discard the header line*/
412  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
413  {
414  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
415  cdEXIT(EXIT_FAILURE);
416  }
417  /*Getting all the collisional transitions data*/
418  for( nCollTrans=0; nCollTrans<intCollTran; nCollTrans++ )
419  {
420  if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
421  {
422  fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
423  cdEXIT(EXIT_FAILURE);
424  }
425 
426  long i = 1;
427  bool lgEOL;
428  long index;
429 
430  index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
431  ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
432  ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
433 
434  /* Indices between the very highest levels seem to be reversed */
435  if( ipHi < ipLo )
436  {
437  ASSERT( ipLo == ngMolLevs - 1);
438  long temp = ipHi;
439  ipHi = ipLo;
440  ipLo = temp;
441  }
442 
443  for( long j=0; j<intCollTemp; j++ )
444  {
445  AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][j] =
446  (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
447  }
448 
449  ASSERT( index == nCollTrans + 1 );
450 
451  if(DEBUGSTATE)
452  {
453  printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo);
454  printf("The collision rates are");
455  for(i=0;i<intCollTemp;i++)
456  {
457  printf("\n %e",AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][i]);
458  }
459  printf("\n");
460  }
461  }
462 
463  intDCollPart = intDCollPart -1;
464  }
465 
466  return;
467 }

Generated for cloudy by doxygen 1.8.4