cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
nemala.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 "taulines.h"
6 #include "trace.h"
7 #include "string.h"
8 #include "thirdparty.h"
9 #include "dense.h"
10 #include "atmdat.h"
11 /*File nemala.cpp was developed by Humeshkar B Nemala as a part of his thesis work during the Summer of 2007*/
12 /* Initially the code has been developed to read in energy levels,radiative and
13  * collisional data from the CHIANTI and LEIDEN databases. The idea is to extend it to more databases.
14  * In the case of the Leiden database there is a single .dat file which has the energy levels information,
15  * radiative and collisional data, with the data corresponding to each collider coming one after the other.
16  * In the case of CHIANTI, the energy levels data, radiative data and collision data are present in seperate files.
17  * While LEIDEN gives collisional rate coefficients, CHIANTI gives collisional strengths.
18  * In the case of CHIANTI only two colliders are used:electrons and protons. They appear as separate files.
19  * The electron collision strengths files are always expected to be there. A flag is set and data processed
20  * if the file on proton collision strengths is available.*/
21 
22 /* There is an initialization file called species.ini which tells Cloudy what kind of data is to be used */
23 /* Structures are created separately to hold the transition data,radiative and collisional data */
24 /* The collisional structures are different for different databases depending upon whether */
25 /* collisional strengths or collisional rate coefficients are used.Finally a superstructure is constructed to hold */
26 /* the total collisional rate obtained by considering all the colliders */
27 /* The colliders considered are electron,proton,Atomic Hydrogen,He,He+,He++,Ortho Molecular Hydrogen,Para Molecular Hydrogen and Molecular Hydrogen */
28 void database_readin(void);
29 void states_popfill( void);
30 void states_nelemfill(void);
31 void database_prep(int);
32 void emislines_fillredis(void);
33 STATIC void states_propprint(void);
34 STATIC int getAtNo(char []);
35 
36 #define DEBUGSTATE false
37 
39 {
40 
41  DEBUG_ENTRY( "AddLine2Stack()" );
42 
44  atmolEmis[linesAdded2].Aul = Aul;
45  atmolEmis[linesAdded2].tran = trans;
46  linesAdded2++;
48  return &atmolEmis[linesAdded2-1];
49 }
50 
51 void Nemala_Start( void )
52 {
53  int i,j,los,intNoSp;
54 
55  FILE *atmolDATA;
56  char *chToken,*chAtNo,*chIonStg;
57 
58  char chLine[FILENAME_PATH_LENGTH_2],
59  chDLine[FILENAME_PATH_LENGTH_2];
60 
61  static int nCalled = 0;
62 
63  DEBUG_ENTRY( "Nemala_Start()" );
64 
65  linesAdded2 = 0;
66 
67  /* only do this once. */
68  if( nCalled > 0 )
69  {
70  return;
71  }
72 
73  /* this is first call, increment the nCalled counterso never do this again */
74  ++nCalled;
75 
76  // Reading in the species.ini file: Humeshkar B Nemala
77  if( trace.lgTrace )
78  fprintf( ioQQQ," Nemala_Start opening species.ini:");
79 
80  atmolDATA = open_data( "species.ini", "r" );
81 
82  if( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) == NULL )
83  {
84  fprintf( ioQQQ, " Nemala_Start could not read first line of species.ini.\n");
85  cdEXIT(EXIT_FAILURE);
86  }
87 
88  /* count how many lines are in the file, ignoring all lines
89  * starting with '#':This would give the number of molecules */
90  nSpecies = 0;
91  while( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) != NULL )
92  {
93  /* we want to count the lines that do not start with %
94  * since these contain data */
95  if( (chLine[0] != '#')&&(chLine[0] != '\n')&&(chLine[0] != ' '))
96  ++nSpecies;
97  }
98  if(DEBUGSTATE)
99  printf("The number of species is %li \n",nSpecies);
100 
101  /* now rewind the file so we can read it a second time*/
102  if( fseek( atmolDATA , 0 , SEEK_SET ) != 0 )
103  {
104  fprintf( ioQQQ, " Nemala_Start could not rewind species.ini.\n");
105  cdEXIT(EXIT_FAILURE);
106  }
107 
108  /*Once we know the number of species we initialize a pointer to a one
109  * dimensional array which contains information about the DB Type */
110  lgSpeciesMolecule = (bool *)MALLOC( (unsigned long)(nSpecies)*sizeof(bool));
111 
112  /*Initialization of the Species Structure*/
113  Species = (species *)MALLOC( (unsigned long)nSpecies*sizeof(species));
114 
115  /*Initialization of the collisional rates array structure*/
116  /*Assume that there are 9 colliders*/
117  CollRatesArray = (double****)MALLOC((unsigned long)nSpecies * sizeof(double***));
119  AtmolCollSplines = (CollSplinesArray****)MALLOC((unsigned long)nSpecies *sizeof(CollSplinesArray***));
120 
121  /*Mallocing here takes care of the number of colliders*/
122  for( i=0; i<nSpecies; i++ )
123  {
125  //AtmolCollSplines[i] = (CollSplinesArray***)MALLOC(NUM_COLLIDERS *sizeof(CollSplinesArray**));
126  CollRatesArray[i] = (double***)MALLOC(NUM_COLLIDERS * sizeof(double**));
127  }
128  /*Initializing all the elements of the matrix to 0*/
129  /*For each species and collider*/
130  for( i=0; i<nSpecies; i++ )
131  {
132  for( j=0; j<NUM_COLLIDERS; j++ )
133  {
134  AtmolCollRateCoeff[i][j].ntemps = 0;
135  AtmolCollRateCoeff[i][j].temps = NULL;
136  AtmolCollRateCoeff[i][j].collrates = NULL;
137  }
138  }
139 
140  /*****************************************/
141  i=0;
142  while( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) != NULL )
143  {
144  caps( chLine );
145  if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
146  {
147  ASSERT(i<nSpecies);
148  if(DEBUGSTATE)
149  printf("%s",chLine);
150  strcpy(chDLine, chLine);
151  chToken = strtok(chDLine," ");
152  los = (int)strlen(chToken);
153  if(DEBUGSTATE)
154  printf("The length of the string is %d\n",los);
155  Species[i].chptrSpName = (char *)MALLOC((unsigned long)(los+1) * sizeof(char));
156  strcpy(Species[i].chptrSpName,chToken);
157  Species[i].intAtNo = 0;
158  Species[i].intIonStage = 0;
159  if(DEBUGSTATE)
160  printf("The name of the species is%s \n",Species[i].chptrSpName);
161  if( nMatch("LEID",chLine) )
162  {
163  lgSpeciesMolecule[i] = true;
164  }
165  else if( nMatch("CHIA",chLine))
166  {
167  lgSpeciesMolecule[i] = false;
168  chAtNo = strtok(chToken,"_");
169  Species[i].intAtNo = getAtNo(chAtNo);
170  chIonStg = strtok(NULL,"_");
171  Species[i].intIonStage = atoi(chIonStg)-1;
172  }
173  else
174  TotalInsanity();
175 
176  if(DEBUGSTATE)
177  printf("The species is molecular? %c \n",TorF(lgSpeciesMolecule[i]) );
178  ++i;
179  }
180  }
181 
182  database_readin();
183 
184  /*Setting the population of states to 0*/
185  states_popfill();
186  /*Setting the values of nelem arbitrarily to 1*/
188 
189  /*Setting nelem of the states to an arbitrary value*/
190  /*Loop over species*/
191  for( intNoSp=0; intNoSp<nSpecies; intNoSp++ )
192  {
193  database_prep(intNoSp);
194  }
195 
196  /*Looping over the emission lines and filling in the redistribution function value*/
198 
199  /*To print the states*/
200  if(DEBUGSTATE)
202  return;
203 }
204 
205 /*Filling the redistribution function value*/
207 {
208  DEBUG_ENTRY("emislines_fillredis()");
209  for( long i=0; i<linesAdded2; i++ )
210  {
212  }
213  return;
214 }
215 
216 /*This function fills the nelem value arbitrarily*/
218 {
219  DEBUG_ENTRY( "states_nelemfill()" );
220 
221  for( long i=0; i<nSpecies; i++ )
222  {
223  for( long j=0; j<Species[i].numLevels_max; j++ )
224  {
225  if( lgSpeciesMolecule[i] )
226  {
227  fixit(); /* this is needed because right now, all lines test on
228  * dense.xIonDense[nelem][IonStg] */
229  atmolStates[i][j].nelem = 1;
230  atmolStates[i][j].IonStg = 1;
231  }
232  else
233  {
234  atmolStates[i][j].nelem = Species[i].intAtNo + 1;
236  }
237  }
238  }
239  return;
240 }
241 
242 /*This function prints the various properties of states*/
244 {
245  DEBUG_ENTRY( "states_propprint()" );
246 
247  for( long i=0; i<nSpecies; i++ )
248  {
249  printf("The species is %s \n",Species[i].chptrSpName);
250  printf("The data output is in the following format \n");
251  printf("Label Energy St.wt Pop Lifetime\n");
252 
253  for( long j=0; j<Species[i].numLevels_max; j++ )
254  {
255  printf("This is the %ld state \n",j);
256  printf("%s %f %f %f %e \n",atmolStates[i][j].chLabel,
257  atmolStates[i][j].energy,
258  atmolStates[i][j].g,
259  atmolStates[i][j].Pop,
260  atmolStates[i][j].lifetime);
261  }
262  }
263  return;
264 }
265 
266 
267 void database_prep(int intSpIndex)
268 {
269  realnum fsumAs;
270  int ipHi,ipLo;
271 
272  DEBUG_ENTRY( "database_prep()" );
273 
274  /*Get the lifetimes*/
275  atmolStates[intSpIndex][0].lifetime= BIGFLOAT;
276  for( ipHi=1; ipHi < Species[intSpIndex].numLevels_max; ipHi++ )
277  {
278  fsumAs = SMALLFLOAT;
279  for( ipLo=0; ipLo<ipHi; ipLo++ )
280  {
281  if(atmolTrans[intSpIndex][ipHi][ipLo].Emis!=NULL)
282  fsumAs = fsumAs + atmolTrans[intSpIndex][ipHi][ipLo].Emis->Aul;
283  }
284  atmolStates[intSpIndex][ipHi].lifetime = 1./fsumAs;
285  }
286  return;
287 }
288 
289 /*Separate Routine to read in the molecules*/
290 void database_readin(void)
291 {
292  DEBUG_ENTRY( "database_readin()" );
293 
294  long intNS;
295 
296  atmolStates = (quantumState **)MALLOC((unsigned long)nSpecies * sizeof(quantumState*));
297  atmolTrans = (transition ***)MALLOC((unsigned long)nSpecies * sizeof(transition**));
298 
299  /* nSpecies should be global*/
300  /*The data array containing the species names should also be global*/
301  for( intNS = 0; intNS < nSpecies; intNS++ )
302  {
303  if ( lgSpeciesMolecule[intNS] )
304  {
305  atmdat_lamda_readin( intNS );
306  }
307  else
308  {
309  atmdat_Chianti_readin( intNS );
310  }
311  }
312  return;
313 }
314 
315 STATIC int getAtNo(char *p)
316 {
317  DEBUG_ENTRY("getAtNo()");
318  if(strcmp(p,"H")==0)
319  {
320  return(ipHYDROGEN);
321  }
322  else if(strcmp(p,"HE")== 0)
323  {
324  return(ipHELIUM);
325  }
326  else if(strcmp(p,"C")== 0)
327  {
328  return(ipCARBON);
329  }
330  else if(strcmp(p,"N")== 0)
331  {
332  return(ipNITROGEN);
333  }
334  else if(strcmp(p,"O")== 0)
335  {
336  return(ipOXYGEN);
337  }
338  else if(strcmp(p,"NE")== 0)
339  {
340  return(ipNEON);
341  }
342  else if(strcmp(p,"NA")== 0)
343  {
344  return(ipSODIUM);
345  }
346  else if(strcmp(p,"MG")== 0)
347  {
348  return(ipMAGNESIUM);
349  }
350  else if(strcmp(p,"AL")== 0)
351  {
352  return(ipALUMINIUM);
353  }
354  else if(strcmp(p,"SI")== 0)
355  {
356  return(ipSILICON);
357  }
358  else if(strcmp(p,"P")== 0)
359  {
360  return(ipPHOSPHORUS);
361  }
362  else if(strcmp(p,"S")== 0)
363  {
364  return(ipSULPHUR);
365  }
366  else if(strcmp(p,"CL")== 0)
367  {
368  return(ipCHLORINE);
369  }
370  else if(strcmp(p,"AR")== 0)
371  {
372  return(ipARGON );
373  }
374  else if(strcmp(p,"K")== 0)
375  {
376  return(ipPOTASSIUM);
377  }
378  else if(strcmp(p,"CA")== 0)
379  {
380  return(ipCALCIUM);
381  }
382  else if(strcmp(p,"SC")== 0)
383  {
384  return(ipSCANDIUM);
385  }
386  else if(strcmp(p,"TI")== 0)
387  {
388  return(ipTITANIUM);
389  }
390  else if(strcmp(p,"V")== 0)
391  {
392  return(ipVANADIUM);
393  }
394  else if(strcmp(p,"CR")== 0)
395  {
396  return(ipCHROMIUM );
397  }
398  else if(strcmp(p,"MN")== 0)
399  {
400  return(ipMANGANESE);
401  }
402  else if(strcmp(p,"FE")== 0)
403  {
404  return(ipIRON);
405  }
406  else if(strcmp(p,"CO")== 0)
407  {
408  return(ipCOBALT);
409  }
410  else if(strcmp(p,"NI")== 0)
411  {
412  return(ipNICKEL);
413  }
414  else if(strcmp(p,"CU")== 0)
415  {
416  return(ipCOPPER);
417  }
418  else if(strcmp(p,"ZN")== 0)
419  {
420  return(ipZINC);
421  }
422  else
423  {
424  TotalInsanity();
425  }
426 }

Generated for cloudy by doxygen 1.8.1.2