cloudy  trunk
species.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 "taulines.h"
6 #include "trace.h"
7 #include "string.h"
8 #include "input.h"
9 #include "thirdparty.h"
10 #include "dense.h"
11 #include "atmdat.h"
12 #include "mole.h"
13 #include "elementnames.h"
14 #include "version.h"
15 
16 /*File nemala.cpp was developed by Humeshkar B Nemala as a part of his thesis work during the Summer of 2007*/
17 /* Initially the code has been developed to read in energy levels,radiative and
18  * collisional data from the CHIANTI and LEIDEN databases. The idea is to extend it to more databases.
19  * In the case of the Leiden database there is a single .dat file which has the energy levels information,
20  * radiative and collisional data, with the data corresponding to each collider coming one after the other.
21  * In the case of CHIANTI, the energy levels data, radiative data and collision data are present in seperate files.
22  * While LEIDEN gives collisional rate coefficients, CHIANTI gives collisional strengths.
23  * In the case of CHIANTI only two colliders are used:electrons and protons. They appear as separate files.
24  * The electron collision strengths files are always expected to be there. A flag is set and data processed
25  * if the file on proton collision strengths is available.*/
26 
27 /* There is an initialization file called species.ini which tells Cloudy what kind of data is to be used */
28 /* Structures are created separately to hold the transition data,radiative and collisional data */
29 /* The collisional structures are different for different databases depending upon whether */
30 /* collisional strengths or collisional rate coefficients are used.Finally a superstructure is constructed to hold */
31 /* the total collisional rate obtained by considering all the colliders */
32 /* The colliders considered are electron,proton,Atomic Hydrogen,He,He+,He++,Ortho Molecular Hydrogen,Para Molecular Hydrogen and Molecular Hydrogen */
33 STATIC void states_popfill(void);
34 STATIC void states_nelemfill(void);
35 STATIC void database_prep(int);
36 STATIC void set_fractionation( species *sp );
37 STATIC void states_propprint(void);
38 /*SpeciesJunk set all elements of species struc to dangerous values */
39 STATIC void SpeciesJunk( species *sp );
40 
41 #define DEBUGSTATE false
42 void database_readin( void )
43 {
44  int i,intNoSp;
45 
46  FILE *ioMASTERLIST, *ioVERSION;
47 
48  char *chToken;
49 
50  char chLine[FILENAME_PATH_LENGTH_2],
51  chDLine[FILENAME_PATH_LENGTH_2],
52  chPath[FILENAME_PATH_LENGTH_2] = "";
53 
54  const int MAX_NUM_SPECIES = 1000;
55 
56  char chLabels[MAX_NUM_SPECIES][CHARS_SPECIES];
57  char chLabelsOrig[MAX_NUM_SPECIES][CHARS_SPECIES];
58  char chPaths[MAX_NUM_SPECIES][FILENAME_PATH_LENGTH_2];
59 
60  static int nCalled = 0;
61  long nSpeciesLAMDA, nSpeciesSTOUT, nSpeciesCHIANTI;
62 
63  DEBUG_ENTRY( "database_readin()" );
64 
65  /* only do this once. */
66  if( nCalled > 0 )
67  {
68  return;
69  }
70 
71  /* this is first call, increment the nCalled counterso never do this again */
72  ++nCalled;
73 
74  // read masterlists, count number of species
75  nSpecies = 0;
76 
78  //
79  // Read LAMDA masterlist
80  //
82 
83  /* count how many lines are in the file, ignoring all lines
84  * starting with '#':This would give the number of molecules */
85  nSpeciesLAMDA = 0;
86 
87  if( atmdat.lgLamdaOn )
88  {
89  long numModelsNotUsed = 0;
90  strcpy( chPath, "lamda" );
91  strcat( chPath, input.chDelimiter );
92  strcat( chPath, "masterlist" );
93  ioMASTERLIST = open_data( chPath, "r" );
94 
95  if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
96  {
97  fprintf( ioQQQ, " database_readin could not read first line of LAMDA masterlist.\n");
98  cdEXIT(EXIT_FAILURE);
99  }
100 
101  do
102  {
103  if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
104  {
105  strcpy(chDLine, chLine);
106  chToken = strtok(chDLine," \t\n");
107  if( findspecies( chToken ) != null_mole ||
108  ( chToken[1]=='-' && findspecies( chToken+2 ) != null_mole ) )
109  {
110  ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
111  ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
112  ASSERT( strlen(chToken) < CHARS_SPECIES );
113  strcpy( chLabels[nSpecies], chToken );
114  chLabels[nSpecies][CHARS_SPECIES-1] = '\0';
115 
116  // path is, for example, LAMDA/no.dat
117  strcpy( chPaths[nSpecies], "lamda" );
118  strcat( chPaths[nSpecies], input.chDelimiter );
119  chToken = strtok( NULL," \t\n" );
120  strcat( chPaths[nSpecies], chToken );
121  ++nSpecies;
122  ++nSpeciesLAMDA;
123  }
124  else
125  ++numModelsNotUsed;
126  }
127  }
128  while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
129 
130  /* \todo 1 - save this and stuff as note since not really a "PROBLEM" but worth reporting */
131  //if( !t_version::Inst().lgRelease && numModelsNotUsed > 0 )
132  // fprintf( ioQQQ, "\n PROBLEM - %li LAMDA models could not be found in chemistry network.\n\n\n", numModelsNotUsed );
133 
134  fclose(ioMASTERLIST);
135  }
136 
138  //
139  // Read CDMS/JPL masterlist
140  //
141  // These data files are in LAMDA format
142  //
144 
145  if( atmdat.lgCalpgmOn )
146  {
147  long numModelsNotUsed = 0;
148  strcpy( chPath, "cdms+jpl" );
149  strcat( chPath, input.chDelimiter );
150  strcat( chPath, "masterlist" );
151  ioMASTERLIST = open_data( chPath, "r" );
152 
153  if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
154  {
155  fprintf( ioQQQ, " database_readin could not read first line of CDMS/JPL masterlist.\n");
156  cdEXIT(EXIT_FAILURE);
157  }
158 
159  do
160  {
161  if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
162  {
163  strcpy(chDLine, chLine);
164  chToken = strtok(chDLine," \t\n");
165  // hacks for alternative dialects...
166  if( strcmp( chToken, "SH" ) == 0 )
167  strcpy( chToken, "HS" );
168  if( strcmp( chToken, "SH+" ) == 0 )
169  strcpy( chToken, "HS+" );
170  if( strcmp( chToken, "CCH" ) == 0 )
171  strcpy( chToken, "C2H" );
172  if( findspecies( chToken ) != null_mole ||
173  ( chToken[1]=='-' && findspecies( chToken+2 ) != null_mole ) )
174  {
175  ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
176  ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
177  strcpy( chLabels[nSpecies], chToken );
178  chLabels[nSpecies][CHARS_SPECIES-1] = '\0';
179 
180  strcpy( chPaths[nSpecies], "cdms+jpl" );
181  strcat( chPaths[nSpecies], input.chDelimiter );
182  chToken = strtok( NULL," \t\n" );
183  strcat( chPaths[nSpecies], chToken );
184  ++nSpecies;
185  ++nSpeciesLAMDA;
186  }
187  else
188  ++numModelsNotUsed;
189  }
190  }
191  while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
192 
193  // if( !t_version::Inst().lgRelease && numModelsNotUsed > 0 )
194  // fprintf( ioQQQ, "\nPROBLEM - %li CDMS/JPL models could not be found in chemistry network.\n\n",
195  // numModelsNotUsed );
196 
197  fclose(ioMASTERLIST);
198  }
199 
201  //
202  // Read STOUT masterlist and VERSION
203  //
205  nSpeciesSTOUT = 0;
206  if( atmdat.lgStoutOn )
207  {
208  // default location of Stout masterlist file
209  strcpy( chPath, "stout" );
210  strcat( chPath, input.chDelimiter );
211  strcat( chPath, "masterlist" );
212  strcat( chPath, input.chDelimiter );
213 
214  strcat( chPath, atmdat.chStoutFile );
215 
216  // first try local directory, then data/SED
217  if( (ioMASTERLIST = fopen( atmdat.chStoutFile , "r" ) ) == NULL )
218  {
219  ioMASTERLIST = open_data( chPath, "r" );
220  }
221 
222  // magic number
223  if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
224  {
225  fprintf( ioQQQ, " database_readin could not read first line of stout.ini.\n");
226  cdEXIT(EXIT_FAILURE);
227  }
228 
229  bool lgEOLST;
230  long int ipST = 1;
231  long int nYrRdST = (long)FFmtRead(chLine,&ipST,sizeof(chLine),&lgEOLST);
232  long int nMonRdST = (long)FFmtRead(chLine,&ipST,sizeof(chLine),&lgEOLST);
233  long int nDayRdST = (long)FFmtRead(chLine,&ipST,sizeof(chLine),&lgEOLST);
234 
235  static long int nYrST =11 , nMonST = 10, nDayST = 25;
236  if( ( nYrRdST != nYrST ) || ( nMonRdST != nMonST ) || ( nDayRdST != nDayST ) )
237  {
238  fprintf( ioQQQ,
239  " I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
240  nYrST , nMonST , nDayST , nYrRdST , nMonRdST , nDayRdST );
241  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
242  cdEXIT(EXIT_FAILURE);
243  }
244  if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
245  {
246  fprintf( ioQQQ, " database_readin could not read first line of CHIANTI masterlist.\n");
247  cdEXIT(EXIT_FAILURE);
248  }
249 
250  do
251  {
252  strcpy(chDLine, chLine);
253  chToken = strtok(chDLine," \n");
254  if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
255  {
256  ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
257  ASSERT( nSpeciesSTOUT + 1 <= MAX_NUM_SPECIES );
258  strcpy( chLabels[nSpecies], chToken );
259  strcpy( chLabelsOrig[nSpecies], chLabels[nSpecies] );
260 
261  char *chElement, chTokenTemp[7];
262  strcpy( chTokenTemp, chToken );
263  chElement = strtok(chTokenTemp," \n");
264  chElement = strtok(chTokenTemp,"_");
265  uncaps( chElement );
266 
267 // printf("STOUT:%s\n",chLabels[nSpecies]);
268 
269  // path is, for example, CHIANTI/ar/ar_10/ar_10
270  // we will append extensions later
271  strcpy( chPaths[nSpecies], "stout" );
272  strcat( chPaths[nSpecies], input.chDelimiter );
273  strcat( chPaths[nSpecies], chElement );
274  strcat( chPaths[nSpecies], input.chDelimiter );
275  strcat( chPaths[nSpecies], chLabels[nSpecies] );
276  strcat( chPaths[nSpecies], input.chDelimiter );
277  strcat( chPaths[nSpecies], chLabels[nSpecies] );
278 
279  ASSERT( isalpha(chToken[0]) );
280  long cursor=0;
281  chLabels[nSpecies][0] = chToken[0];
282  if( isalpha(chToken[1]) )
283  {
284  chLabels[nSpecies][1] = chToken[1];
285  cursor = 2;
286  }
287  else
288  {
289  chLabels[nSpecies][1] = ' ';
290  cursor = 1;
291  }
292 
293  ASSERT( chToken[cursor]=='_' );
294  ++cursor;
295  ASSERT( isdigit(chToken[cursor]) );
296 
297  if( isdigit(chToken[cursor+1]) )
298  {
299  chLabels[nSpecies][2] = chToken[cursor++];
300  chLabels[nSpecies][3] = chToken[cursor++];
301  }
302  else
303  {
304  chLabels[nSpecies][2] = ' ';
305  chLabels[nSpecies][3] = chToken[cursor++];
306  }
307  chLabels[nSpecies][4] = '\0';
308  ASSERT( chToken[cursor]=='\0' || chToken[cursor]=='d' );
309 
310  // now capitalize the first letter
311  chLabels[nSpecies][0] = toupper( chLabels[nSpecies][0] );
312  ++nSpecies;
313  ++nSpeciesSTOUT;
314  }
315  }
316  while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
317  fclose(ioMASTERLIST);
318  }
319 
320 
322  //
323  // Read CHIANTI masterlist and VERSION
324  //
326 
327  nSpeciesCHIANTI = 0;
328 
329  if( atmdat.lgChiantiOn )
330  {
331  char chPathSave[FILENAME_PATH_LENGTH_2];
332  strcpy( chPath, "chianti" );
333  strcat( chPath, input.chDelimiter );
334  //Preserve the path /chianti/ with chPathSave
335  //Start reading in the chianti version number
336  strcpy( chPathSave , chPath );
337  strcat(chPath,"VERSION");
338  ioVERSION = open_data(chPath,"r");
339  if( read_whole_line( chLine , (int)sizeof(chLine) , ioVERSION ) == NULL )
340  {
341  fprintf( ioQQQ, " database_readin could not read first line of the Chianti VERSION.\n");
342  cdEXIT(EXIT_FAILURE);
343  }
344  fclose(ioVERSION);
345  // chianti version - string since can contain letters
346  strncpy(atmdat.chVersion,chLine,atmdat.iVersionLength);
347  // remove newline we may captured
348  long len = strlen(atmdat.chVersion);
349  if( atmdat.chVersion[len-1] == '\n' )
350  atmdat.chVersion[len-1] = '\0';
351  // may have been null earlier in string, but make sure terminated at limit
352  atmdat.chVersion[atmdat.iVersionLength-1] = '\0';
353  //Restore the previous chPath
354  strcpy(chPath,chPathSave);
355  // Read in the masterlist
356  strcat( chPath, "masterlist" );
357  strcat( chPath, input.chDelimiter );
358  // save copy
359  strcpy( chPathSave , chPath );
360 
361  // our subset of Chianti
362  strcat( chPath, atmdat.chCloudyChiantiFile );
363 
364  // first try local directory, then data/chianti
365  if( (ioMASTERLIST = fopen( atmdat.chCloudyChiantiFile , "r" ) ) == NULL )
366  {
367  ioMASTERLIST = open_data( chPath, "r" );
368  }
369 
370  // magic number
371  if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
372  {
373  fprintf( ioQQQ, " database_readin could not read first line of CloudyChianti.ini.\n");
374  cdEXIT(EXIT_FAILURE);
375  }
376 
377  bool lgEOL;
378  long int ip = 1;
379  long int nYrRd = (long)FFmtRead(chLine,&ip,sizeof(chLine),&lgEOL);
380  long int nMonRd = (long)FFmtRead(chLine,&ip,sizeof(chLine),&lgEOL);
381  long int nDayRd = (long)FFmtRead(chLine,&ip,sizeof(chLine),&lgEOL);
382 
383  static long int nYr=11 , nMon = 10, nDay = 3;
384  if( ( nYrRd != nYr ) || ( nMonRd != nMon ) || ( nDayRd != nDay ) )
385  {
386  fprintf( ioQQQ,
387  " database_readin: the version of CloudyChianti.ini is not the current version.\n" );
388  fprintf( ioQQQ,
389  " database_readin obtain the current version from the Cloudy web site.\n" );
390  fprintf( ioQQQ,
391  " I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
392  nYr , nMon , nDay , nYrRd , nMonRd , nDayRd );
393  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
394  cdEXIT(EXIT_FAILURE);
395  }
396 
397  if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
398  {
399  fprintf( ioQQQ, " database_readin could not read first line of CHIANTI masterlist.\n");
400  cdEXIT(EXIT_FAILURE);
401  }
402 
403  do
404  {
405 
406  if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
407  {
408  strcpy(chDLine, chLine);
409  chToken = strtok(chDLine," \n");
410 
411  fixit(); //insert logic here to exclude some ions (for example, iso sequences)
412  // exclude for now the satellite lines (denoted by a "d" after the label
413  //if( chToken[3]=='d' || chToken[4]=='d' || chToken[5]=='d' )
414  if( chToken[3]!='d' && chToken[4]!='d' && chToken[5]!='d' )
415  {
416  ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
417  ASSERT( nSpeciesCHIANTI + 1 <= MAX_NUM_SPECIES );
418  strcpy( chLabels[nSpecies], chToken );
419  strcpy( chLabelsOrig[nSpecies], chLabels[nSpecies]);
420 
421  bool skipSpecies = false;
422 
423  //Check for duplicate species with Stout
424  for( int j = nSpeciesLAMDA; j < nSpecies; j++)
425  {
426  if( strcmp( chLabelsOrig[j], chLabelsOrig[nSpecies] ) == 0)
427  {
428  printf("Skipping the Chianti version of %s, using Stout version\n",chLabels[nSpecies]);
429  skipSpecies = true;
430  break;
431  }
432  }
433  if( skipSpecies)
434  continue;
435 
436  char *chElement, chTokenTemp[7];
437  strcpy( chTokenTemp, chToken );
438  chElement = strtok(chTokenTemp," \n");
439  chElement = strtok(chTokenTemp,"_");
440  uncaps( chElement );
441 
442  // path is, for example, CHIANTI/ar/ar_10/ar_10
443  // we will append extensions later
444  strcpy( chPaths[nSpecies], "chianti" );
445  strcat( chPaths[nSpecies], input.chDelimiter );
446  strcat( chPaths[nSpecies], chElement );
447  strcat( chPaths[nSpecies], input.chDelimiter );
448  strcat( chPaths[nSpecies], chLabels[nSpecies] );
449  strcat( chPaths[nSpecies], input.chDelimiter );
450  strcat( chPaths[nSpecies], chLabels[nSpecies] );
451 
452  ASSERT( isalpha(chToken[0]) );
453  long cursor=0;
454  chLabels[nSpecies][0] = chToken[0];
455  if( isalpha(chToken[1]) )
456  {
457  chLabels[nSpecies][1] = chToken[1];
458  cursor = 2;
459  }
460  else
461  {
462  chLabels[nSpecies][1] = ' ';
463  cursor = 1;
464  }
465 
466  ASSERT( chToken[cursor]=='_' );
467  ++cursor;
468  ASSERT( isdigit(chToken[cursor]) );
469 
470  if( isdigit(chToken[cursor+1]) )
471  {
472  chLabels[nSpecies][2] = chToken[cursor++];
473  chLabels[nSpecies][3] = chToken[cursor++];
474  }
475  else
476  {
477  chLabels[nSpecies][2] = ' ';
478  chLabels[nSpecies][3] = chToken[cursor++];
479  }
480  chLabels[nSpecies][4] = '\0';
481  ASSERT( chToken[cursor]=='\0' || chToken[cursor]=='d' );
482 
483  // now capitalize the first letter
484  chLabels[nSpecies][0] = toupper( chLabels[nSpecies][0] );
485  ++nSpecies;
486  ++nSpeciesCHIANTI;
487  }
488  }
489  }
490  while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
491 
492  fclose(ioMASTERLIST);
493  }
494 
495  /* no species found, nothing to do */
496  if( nSpecies==0 )
497  return;
498 
499  /*Initialization of the dBaseSpecies Structure*/
500  dBaseSpecies = (species *)MALLOC( (unsigned long)nSpecies*sizeof(species));
501 
502  /*Initialization of the collisional rates array structure*/
503  AtmolCollRateCoeff.reserve( nSpecies );
504  AtmolCollSplines = (CollSplinesArray****)MALLOC((unsigned long)nSpecies *sizeof(CollSplinesArray***));
505  StoutCollData = (StoutColls****)MALLOC((unsigned long)nSpecies *sizeof(StoutColls***));
506 
507  /*Mallocing here takes care of the number of colliders*/
508  for( i=0; i<nSpecies; i++ )
509  {
510  AtmolCollRateCoeff.reserve( i, ipNCOLLIDER );
511  }
512  AtmolCollRateCoeff.alloc();
513 
514  // malloc state and transition arrays
515  dBaseStates.resize(nSpecies);
516  ipdBaseTrans.resize(nSpecies);
517 
518  for( i = 0; i < nSpecies; i++ )
519  {
520  dBaseTrans.push_back(TransitionList("dBaseTrans",&dBaseStates[i]));
521  SpeciesJunk( &dBaseSpecies[i] );
522  // label should be a minimum of 4 characters long
523  size_t los = max(4,strlen(chLabels[i]));
524  ASSERT( los >= 4 && los <= 7 );
525  dBaseSpecies[i].chLabel = new char[los+1];
526  strcpy(dBaseSpecies[i].chLabel,chLabels[i]);
527  dBaseSpecies[i].lgActive = true;
528 
529  /* set type and isotopologue fractions */
530  set_fractionation( &dBaseSpecies[i] );
531 
532  // set_fractionation trims off "p-","o-", etc. Now have set label. Check size.
533  los = (int)strlen( dBaseSpecies[i].chLabel );
534  ASSERT( los < CHARS_SPECIES );
535 
536  // pad label to at least four characters.
537  if( dBaseSpecies[i].chLabel[2]=='\0' )
538  {
539  dBaseSpecies[i].chLabel[2]=' ';
540  dBaseSpecies[i].chLabel[3]=' ';
541  dBaseSpecies[i].chLabel[4]='\0';
542  }
543  else if( dBaseSpecies[i].chLabel[3]=='\0' )
544  {
545  dBaseSpecies[i].chLabel[3]=' ';
546  dBaseSpecies[i].chLabel[4]='\0';
547  }
548 
549  if( i<nSpeciesLAMDA )
550  {
551  // Read in LAMDA data files
552  atmdat_LAMDA_readin( i, chPaths[i] );
553  }
554  else if( i < nSpeciesLAMDA + nSpeciesSTOUT )
555  {
556  atmdat_STOUT_readin( i, chPaths[i] );
557  }
558  else if( i < nSpeciesLAMDA + nSpeciesSTOUT + nSpeciesCHIANTI )
559  {
560  // Read in CHIANTI data files
561  atmdat_CHIANTI_readin( i, chPaths[i] );
562  }
563  else
564  TotalInsanity();
565  }
566 
567  states_popfill();
568  states_nelemfill();
569 
570  /*Setting nelem of the states to an arbitrary value*/
571  /*Loop over species*/
572  for( intNoSp=0; intNoSp<nSpecies; intNoSp++ )
573  {
574  database_prep(intNoSp);
575  AllTransitions.push_back(dBaseTrans[intNoSp]);
576  }
577 
578  /*To print the states*/
579  if(DEBUGSTATE)
580  states_propprint();
581  return;
582 }
583 
584 STATIC void set_fractionation( species *sp )
585 {
586  DEBUG_ENTRY("set_fractionation()");
587 
588  char chToken[3];
589 
590  sp->fracIsotopologue = 1.f;
591  //types include "p-", "o-", "e-", and "a-"
592  strncpy( chToken, sp->chLabel, 2 );
593  chToken[2] = '\0';
594  if( strcmp( "p-", chToken )==0 )
595  sp->fracType = 0.25f;
596  else if( strcmp( "o-", chToken )==0 )
597  sp->fracType = 0.75f;
598  else if( strcmp( "e-", chToken )==0 )
599  sp->fracType = 0.5f;
600  else if( strcmp( "a-", chToken )==0 )
601  sp->fracType = 0.5f;
602  else
603  sp->fracType = 1.0f;
604 
605  fixit(); // what fraction should e-type and a-type Methanol have? Assume 50/50 for now.
606 
607  // Now scrape the type specifier off the label.
608  if( sp->chLabel[1]=='-')
609  memmove(sp->chLabel,sp->chLabel+2,strlen(sp->chLabel+2)+1);
610 
611  return;
612 }
613 
614 /*This function zeros the population of all states */
615 STATIC void states_popfill( void)
616 {
617  DEBUG_ENTRY( "states_popfill()" );
618 
619  for( long i=0; i<nSpecies; i++)
620  {
621  for( long j=0; j<dBaseSpecies[i].numLevels_max; j++)
622  {
623  dBaseStates[i][j].Pop() = 0.;
624  }
625  }
626  return;
627 }
628 
629 /*This function fills the nelem and IonStg fields */
630 STATIC void states_nelemfill(void)
631 {
632  DEBUG_ENTRY( "states_nelemfill()" );
633 
634  for( long i=0; i<nSpecies; i++ )
635  {
636  long nelem = 0, IonStg;
637  char chLabelChemical[CHARS_SPECIES];
638 
639  if( dBaseSpecies[i].lgMolecular )
640  {
641  fixit();
642  /* these should never be used if lgMolecular
643  *set to dangerous values instead of unity. */
644  nelem = -1;
645  IonStg = -1;
646  strcpy( chLabelChemical, dBaseSpecies[i].chLabel );
647  }
648  else
649  {
650  char chToken[3];
651  strncpy( chToken, dBaseSpecies[i].chLabel, 2 );
652  chToken[2] = '\0';
653  strcpy( chLabelChemical, chToken );
654  if( chLabelChemical[1]==' ' )
655  chLabelChemical[1] = '\0';
656  for( long ipElement=0; ipElement<LIMELM; ipElement++ )
657  {
658  if( strcmp( elementnames.chElementSym[ipElement], chToken )==0 )
659  {
660  nelem = ipElement + 1;
661  break;
662  }
663  }
664  ASSERT( nelem > 0 && nelem <= LIMELM );
665  strncpy( chToken, dBaseSpecies[i].chLabel + 2, 2 );
666  IonStg = atoi(chToken);
667  char chStage[5] = {'\0'};
668  if( IonStg==2 )
669  sprintf( chStage, "+" );
670  else if( IonStg>1 )
671  sprintf( chStage, "+%li", IonStg-1 );
672  strcat( chLabelChemical, chStage );
673  ASSERT( IonStg >= 1 && IonStg <= nelem+1 );
674  //Prevent importing of iso-sequences from Chianti
675  if( nelem - IonStg < NISO )
676  {
677  fprintf(ioQQQ, " PROBLEM: Cannot use Chianti model for %s%li\n",elementnames.chElementSym[nelem-1],IonStg);
678  fprintf(ioQQQ, " Iso-sequences are handled by our own model.\n");
679  cdEXIT(EXIT_FAILURE);
680  }
681  dBaseSpecies[i].fmolweight = dense.AtomicWeight[nelem-1];
682  // do not evaluate our cooling if we are using Chianti for this species
683 
684  if( dBaseTrans[i].chLabel() == "Chianti" )
685  {
686  dense.lgIonChiantiOn[nelem-1][IonStg-1] = true;
687  }
688  else if( dBaseTrans[i].chLabel() == "Stout" )
689  {
690  dense.lgIonStoutOn[nelem-1][IonStg-1] = true;
691  }
692  else
693  {
694  TotalInsanity();
695  }
696 
697  if( atmdat.lgChiantiHybrid || atmdat.lgStoutHybrid )
698  {
699  // used in cool_dima to indicate whether to include line
700  // with shorter wl than these databases
701  dense.maxWN[nelem-1][IonStg-1] = dBaseSpecies[i].maxWN;
702  }
703  else
704  {
705  dense.maxWN[nelem-1][IonStg-1] = 0.;
706  }
707  }
708 
709  molecule *sp = findspecies(chLabelChemical);
710  if( sp == null_mole )
711  {
712  dBaseSpecies[i].index = INT_MAX;
713  if( nelem-1 >= ipHYDROGEN && dense.lgElmtOn[nelem-1] )
714  fprintf(ioQQQ," PROBLEM: could not find species %li - %s\n",i,
715  chLabelChemical );
716  }
717  else
718  {
719  dBaseSpecies[i].index = sp->index;
720  mole.species[ sp->index ].levels = &dBaseStates[i];
721  mole.species[ sp->index ].lines = &dBaseTrans[i];
722  }
723 
724  for( long j=0; j<dBaseSpecies[i].numLevels_max; j++ )
725  {
726  dBaseStates[i][j].nelem() = nelem;
727  dBaseStates[i][j].IonStg() = IonStg;
728  }
729  }
730  return;
731 }
732 
733 /*This function prints the various properties of states*/
734 STATIC void states_propprint(void)
735 {
736  DEBUG_ENTRY( "states_propprint()" );
737 
738  for( long i=0; i<nSpecies; i++ )
739  {
740  printf("The species is %s \n",dBaseSpecies[i].chLabel);
741  printf("The data output is in the following format \n");
742  printf("Label Energy St.wt Pop Lifetime\n");
743 
744  for( long j=0; j<dBaseSpecies[i].numLevels_max; j++ )
745  {
746  printf("This is the %ld state \n",j);
747  printf("%s %f %f %f %e \n",dBaseStates[i][j].chLabel(),
748  dBaseStates[i][j].energy().WN(),
749  dBaseStates[i][j].g(),
750  dBaseStates[i][j].Pop(),
751  dBaseStates[i][j].lifetime());
752  }
753  }
754  return;
755 }
756 
757 
758 STATIC void database_prep(int intSpIndex)
759 {
760  vector<realnum> fsumAs(dBaseSpecies[intSpIndex].numLevels_max,SMALLFLOAT);
761 
762  DEBUG_ENTRY( "database_prep()" );
763 
764  /*Get the lifetimes*/
765  for( EmissionList::iterator em = dBaseTrans[intSpIndex].Emis().begin();
766  em != dBaseTrans[intSpIndex].Emis().end(); ++em)
767  {
768  fsumAs[(*em).Tran().ipHi()] += (*em).Aul();
769  (*em).iRedisFun() = ipPRD;
770  }
771 
772  dBaseStates[intSpIndex][0].lifetime()= BIGFLOAT;
773  for( int ipHi=1; ipHi < dBaseSpecies[intSpIndex].numLevels_max; ipHi++ )
774  {
775  dBaseStates[intSpIndex][ipHi].lifetime() = 1./fsumAs[ipHi];
776  }
777  return;
778 }
779 
780 /*SpeciesJunk set all elements of species struc to dangerous values */
781 STATIC void SpeciesJunk( species *sp )
782 {
783  sp->chLabel = NULL;
784  set_NaN(sp->fmolweight);
785  set_NaN(sp->fracIsotopologue );
786  set_NaN(sp->fracType);
787  sp->lgMolecular = false;
788  sp->numLevels_local = -INT_MAX;
789  sp->numLevels_max = -INT_MAX;
790 
791  return;
792 }
t_atmdat::lgLamdaOn
bool lgLamdaOn
Definition: atmdat.h:239
AtmolCollSplines
CollSplinesArray **** AtmolCollSplines
Definition: taulines.cpp:19
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
t_dense::lgIonChiantiOn
bool lgIonChiantiOn[LIMELM][LIMELM+1]
Definition: dense.h:128
dense
t_dense dense
Definition: dense.cpp:24
chElement
static char chElement[NPUNLM][5]
Definition: save_colden.cpp:12
elementnames.h
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
database_readin
void database_readin(void)
Definition: species.cpp:42
AllTransitions
vector< TransitionList > AllTransitions
Definition: taulines.cpp:8
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
t_atmdat::chStoutFile
char chStoutFile[FILENAME_PATH_LENGTH]
Definition: atmdat.h:268
STATIC
#define STATIC
Definition: cddefines.h:97
mole.h
molecule::index
int index
Definition: mole.h:169
DEBUGSTATE
#define DEBUGSTATE
Definition: species.cpp:41
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
thirdparty.h
trace.h
ProxyIterator
Definition: proxy_iterator.h:58
lines_service.h
input
t_input input
Definition: input.cpp:12
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
version.h
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
atmdat.h
toupper
char toupper(char c)
Definition: cddefines.h:700
t_species::chLabel
char * chLabel
Definition: cddefines.h:1235
t_atmdat::chCloudyChiantiFile
char chCloudyChiantiFile[FILENAME_PATH_LENGTH]
Definition: atmdat.h:271
t_species::maxWN
double maxWN
Definition: cddefines.h:1257
CHARS_SPECIES
@ CHARS_SPECIES
Definition: cddefines.h:274
t_atmdat::lgStoutOn
bool lgStoutOn
Definition: atmdat.h:241
SpeciesJunk
STATIC void SpeciesJunk(species *sp)
Definition: species.cpp:781
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_atmdat::lgChiantiHybrid
bool lgChiantiHybrid
Definition: atmdat.h:233
t_CollSplinesArray
Definition: cddefines.h:1275
t_species::numLevels_local
long numLevels_local
Definition: cddefines.h:1241
t_species::index
long index
Definition: cddefines.h:1237
dense.h
t_input::chDelimiter
char chDelimiter[3]
Definition: input.h:42
database_prep
STATIC void database_prep(int)
Definition: species.cpp:758
mole
t_mole_local mole
Definition: mole.cpp:7
t_atmdat::lgChiantiOn
bool lgChiantiOn
Definition: atmdat.h:231
species
struct t_species species
Definition: cddefines.h:1224
cddefines.h
null_mole
molecule * null_mole
Definition: mole_species.cpp:64
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
uncaps
void uncaps(char *chCard)
Definition: service.cpp:263
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
set_NaN
void set_NaN(sys_float &x)
Definition: cpu.cpp:682
t_species::lgActive
bool lgActive
Definition: cddefines.h:1255
ipNCOLLIDER
@ ipNCOLLIDER
Definition: collision.h:18
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_atmdat::iVersionLength
static const int iVersionLength
Definition: atmdat.h:263
t_StoutColls
Definition: cddefines.h:1294
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
BIGFLOAT
const UNUSED realnum BIGFLOAT
Definition: cpu.h:189
StoutCollData
StoutColls **** StoutCollData
Definition: taulines.cpp:20
t_species::fmolweight
realnum fmolweight
Definition: cddefines.h:1243
ipPRD
const int ipPRD
Definition: cddefines.h:290
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
t_species::fracType
realnum fracType
Definition: cddefines.h:1249
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
TransitionList
Definition: transition.h:274
states_popfill
STATIC void states_popfill(void)
Definition: species.cpp:615
t_atmdat::lgCalpgmOn
bool lgCalpgmOn
Definition: atmdat.h:247
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
ipdBaseTrans
vector< multi_arr< int, 2 > > ipdBaseTrans
Definition: taulines.cpp:16
atmdat_CHIANTI_readin
void atmdat_CHIANTI_readin(long intNS, char *chFileName)
Definition: atmdat_chianti.cpp:595
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
states_propprint
STATIC void states_propprint(void)
Definition: species.cpp:734
fixit
void fixit(void)
Definition: service.cpp:991
t_species
Definition: cddefines.h:1232
t_dense::maxWN
double maxWN[LIMELM][LIMELM+1]
Definition: dense.h:134
t_atmdat::lgStoutHybrid
bool lgStoutHybrid
Definition: atmdat.h:243
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
taulines.h
molecule
Definition: mole.h:132
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
t_dense::lgIonStoutOn
bool lgIonStoutOn[LIMELM][LIMELM+1]
Definition: dense.h:131
set_fractionation
STATIC void set_fractionation(species *sp)
Definition: species.cpp:584
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:15
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
t_species::fracIsotopologue
realnum fracIsotopologue
Definition: cddefines.h:1251
AtmolCollRateCoeff
multi_arr< CollRateCoeffArray, 2 > AtmolCollRateCoeff
Definition: taulines.cpp:18
t_species::numLevels_max
long numLevels_max
Definition: cddefines.h:1239
t_species::lgMolecular
bool lgMolecular
Definition: cddefines.h:1245
atmdat_STOUT_readin
void atmdat_STOUT_readin(long intNS, char *chFileName)
Definition: atmdat_chianti.cpp:17
nSpecies
long int nSpecies
Definition: taulines.cpp:21
NISO
const int NISO
Definition: cddefines.h:261
input.h
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
atmdat_LAMDA_readin
void atmdat_LAMDA_readin(long intNS, char *chFileName)
Definition: atmdat_lamda.cpp:18
t_atmdat::chVersion
char chVersion[iVersionLength]
Definition: atmdat.h:265
states_nelemfill
STATIC void states_nelemfill(void)
Definition: species.cpp:630
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17