cloudy trunk
Loading...
Searching...
No Matches
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 */
33STATIC void states_popfill(void);
34STATIC void states_nelemfill(void);
35STATIC void database_prep(int);
37STATIC void states_propprint(void);
38/*SpeciesJunk set all elements of species struc to dangerous values */
39STATIC void SpeciesJunk( species *sp );
40
41#define DEBUGSTATE false
42void 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],
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");
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");
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");
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 );
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");
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");
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");
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 );
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");
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]));
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 */
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
565 }
566
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)
581 return;
582}
583
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 */
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 */
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");
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 {
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*/
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
758STATIC 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 */
782{
783 sp->chLabel = NULL;
784 set_NaN(sp->fmolweight);
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 atmdat
Definition atmdat.cpp:6
void atmdat_STOUT_readin(long intNS, char *chFileName)
void atmdat_CHIANTI_readin(long intNS, char *chFileName)
void atmdat_LAMDA_readin(long intNS, char *chFileName)
#define DEBUGSTATE
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
struct t_CollSplinesArray CollSplinesArray
@ CHARS_SPECIES
Definition cddefines.h:274
#define MALLOC(exp)
Definition cddefines.h:501
struct t_species species
Definition cddefines.h:1224
const int LIMELM
Definition cddefines.h:258
struct t_StoutColls StoutColls
#define EXIT_FAILURE
Definition cddefines.h:140
char toupper(char c)
Definition cddefines.h:700
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
long max(int a, long b)
Definition cddefines.h:775
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void uncaps(char *chCard)
Definition service.cpp:263
void fixit(void)
Definition service.cpp:991
const int ipPRD
Definition cddefines.h:290
EmissionProxy::iterator iterator
Definition emission.h:317
int index
Definition mole.h:169
@ ipNCOLLIDER
Definition collision.h:18
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
void set_NaN(sys_float &x)
Definition cpu.cpp:682
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_input input
Definition input.cpp:12
t_mole_local mole
Definition mole.cpp:7
molecule * findspecies(const char buf[])
molecule * null_mole
static char chElement[NPUNLM][5]
static double * g
Definition species2.cpp:28
STATIC void states_nelemfill(void)
Definition species.cpp:630
STATIC void database_prep(int)
Definition species.cpp:758
STATIC void set_fractionation(species *sp)
Definition species.cpp:584
STATIC void states_propprint(void)
Definition species.cpp:734
void database_readin(void)
Definition species.cpp:42
STATIC void states_popfill(void)
Definition species.cpp:615
STATIC void SpeciesJunk(species *sp)
Definition species.cpp:781
long numLevels_max
Definition cddefines.h:1239
bool lgMolecular
Definition cddefines.h:1245
long numLevels_local
Definition cddefines.h:1241
realnum fracIsotopologue
Definition cddefines.h:1251
char * chLabel
Definition cddefines.h:1235
realnum fracType
Definition cddefines.h:1249
realnum fmolweight
Definition cddefines.h:1243
long int nSpecies
Definition taulines.cpp:21
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
StoutColls **** StoutCollData
Definition taulines.cpp:20
species * dBaseSpecies
Definition taulines.cpp:14
CollSplinesArray **** AtmolCollSplines
Definition taulines.cpp:19
vector< TransitionList > AllTransitions
Definition taulines.cpp:8