41 #define DEBUGSTATE false
46 FILE *ioMASTERLIST, *ioVERSION;
54 const int MAX_NUM_SPECIES = 1000;
60 static int nCalled = 0;
61 long nSpeciesLAMDA, nSpeciesSTOUT, nSpeciesCHIANTI;
89 long numModelsNotUsed = 0;
90 strcpy( chPath,
"lamda" );
92 strcat( chPath,
"masterlist" );
99 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) == NULL )
101 fprintf(
ioQQQ,
" database_readin could not read first line of LAMDA masterlist.\n");
107 if ((chLine[0]!=
'#') && (chLine[0]!=
'\n')&&(chLine[0]!=
'\t')&&(chLine[0]!=
'\r'))
109 strcpy(chDLine, chLine);
110 chToken = strtok(chDLine,
" \t\n");
115 ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
117 strcpy( chLabels[
nSpecies], chToken );
121 strcpy( chPaths[nSpecies],
"lamda" );
123 chToken = strtok( NULL,
" \t\n" );
124 strcat( chPaths[nSpecies], chToken );
132 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) != NULL );
138 fclose(ioMASTERLIST);
146 for(
int i=0; i<nSpeciesLAMDA; i++)
167 strcpy( chPath,
"cdms+jpl" );
169 strcat( chPath,
"masterlist" );
174 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) == NULL )
176 fprintf(
ioQQQ,
" database_readin could not read first line of CDMS/JPL masterlist.\n");
182 if ((chLine[0]!=
'#') && (chLine[0]!=
'\n')&&(chLine[0]!=
'\t')&&(chLine[0]!=
'\r'))
184 strcpy(chDLine, chLine);
185 chToken = strtok(chDLine,
" \t\n");
187 if( strcmp( chToken,
"SH" ) == 0 )
188 strcpy( chToken,
"HS" );
189 if( strcmp( chToken,
"SH+" ) == 0 )
190 strcpy( chToken,
"HS+" );
191 if( strcmp( chToken,
"SD" ) == 0 )
192 strcpy( chToken,
"DS" );
193 if( strcmp( chToken,
"CCH" ) == 0 )
194 strcpy( chToken,
"C2H" );
195 if( strcmp( chToken,
"CCD" ) == 0 )
196 strcpy( chToken,
"C2D" );
197 if( strcmp( chToken,
"^17OO" ) == 0 )
198 strcpy( chToken,
"O^17O" );
199 if( strcmp( chToken,
"H^18O" ) == 0 )
200 strcpy( chToken,
"^18OH" );
201 if( strcmp( chToken,
"HCCD" ) == 0 )
202 strcpy( chToken,
"C2HD" );
203 if( strcmp( chToken,
"^13CCCH" ) == 0 )
204 strcpy( chToken,
"^13CC2H" );
205 if( strcmp( chToken,
"CC^13CH" ) == 0 )
206 strcpy( chToken,
"C2^13CH" );
207 if( strcmp( chToken,
"H^13CCCN" ) == 0 )
208 strcpy( chToken,
"H^13CC2N" );
209 if( strcmp( chToken,
"HCC^13CN" ) == 0 )
210 strcpy( chToken,
"HC2^13CN" );
211 if( strcmp( chToken,
"HCCC^15N" ) == 0 )
212 strcpy( chToken,
"HC3^15N" );
214 if( strcmp( chToken,
"Si^13CC" ) == 0 )
215 strcpy( chToken,
"SiC^13C" );
220 ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
221 strcpy( chLabels[
nSpecies], chToken );
224 strcpy( chPaths[nSpecies],
"cdms+jpl" );
226 chToken = strtok( NULL,
" \t\n" );
227 strcat( chPaths[nSpecies], chToken );
234 fprintf(
ioQQQ,
"Warning: CDMS/JPL species %s not found\n", chToken );
238 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) != NULL );
240 fclose(ioMASTERLIST);
251 vector<long> numLevels(MAX_NUM_SPECIES,0L);
256 strcpy( chPath,
"stout" );
258 strcat( chPath,
"masterlist" );
271 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) == NULL )
273 fprintf(
ioQQQ,
" database_readin could not read first line of stout.ini.\n");
277 bool lgEOL1=
true, lgEOL2=
true, lgEOL3=
true;
278 long int nMonRdST=-1, nDayRdST=-1;
280 long int nYrRdST = (long)
FFmtRead(chLine,&ipST,
sizeof(chLine),&lgEOL1);
283 nMonRdST = (long)
FFmtRead(chLine,&ipST,
sizeof(chLine),&lgEOL2);
285 nDayRdST = (long)
FFmtRead(chLine,&ipST,
sizeof(chLine),&lgEOL3);
289 fprintf(
ioQQQ,
"PROBLEM, there must be three magic numbers on the first line of the stout masterlist file.\n");
293 static long int nYrST =11 , nMonST = 10, nDayST = 25;
294 if( ( nYrRdST != nYrST ) || ( nMonRdST != nMonST ) || ( nDayRdST != nDayST ) )
297 " I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
298 nYrST , nMonST , nDayST , nYrRdST , nMonRdST , nDayRdST );
299 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
302 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) == NULL )
304 fprintf(
ioQQQ,
" database_readin could not read first line of CHIANTI masterlist.\n");
310 strcpy(chDLine, chLine);
315 chToken = strtok(chDLine,
" \t\n");
317 if ((chLine[0]!=
'#') && (chLine[0]!=
'\n')&&(chLine[0]!=
'\t')&&(chLine[0]!=
'\r'))
320 ASSERT( nSpeciesSTOUT + 1 <= MAX_NUM_SPECIES );
323 strcpy( chLabels[
nSpecies], chToken );
324 strcpy( chLabelsOrig[nSpecies], chLabels[nSpecies] );
328 char *chNumLevs = strtok(NULL,
"\n");
329 if( chNumLevs != NULL )
333 long numLevs = (long)
FFmtRead(chNumLevs,&i,
sizeof(chLine),&lgEOL);
343 fprintf(
ioQQQ,
"PROBLEM the limit to the number of levels must be positive, it was %li\n", numLevs);
350 bool skipSpecies =
false;
353 for(
int j = nSpeciesLAMDA; j <
nSpecies; j++)
355 if( strcmp( chLabelsOrig[j], chLabelsOrig[nSpecies] ) == 0)
366 char *chElement, chTokenTemp[7];
367 strcpy( chTokenTemp, chToken );
368 (void) strtok(chTokenTemp,
" \n");
369 chElement = strtok(chTokenTemp,
"_");
374 strcpy( chPaths[nSpecies],
"stout" );
376 strcat( chPaths[nSpecies], chElement );
378 strcat( chPaths[nSpecies], chLabels[nSpecies] );
380 strcat( chPaths[nSpecies], chLabels[nSpecies] );
382 ASSERT( isalpha(chToken[0]) );
385 if( isalpha(chToken[1]) )
396 ASSERT( chToken[cursor]==
'_' );
398 ASSERT( isdigit(chToken[cursor]) );
400 if( isdigit(chToken[cursor+1]) )
402 chLabels[
nSpecies][2] = chToken[cursor++];
403 chLabels[
nSpecies][3] = chToken[cursor++];
408 chLabels[
nSpecies][3] = chToken[cursor++];
411 ASSERT( chToken[cursor]==
'\0' || chToken[cursor]==
'd' );
419 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) != NULL );
420 fclose(ioMASTERLIST);
435 strcpy( chPath,
"chianti" );
439 strcpy( chPathSave , chPath );
440 strcat(chPath,
"VERSION");
442 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioVERSION ) == NULL )
444 fprintf(
ioQQQ,
" database_readin could not read first line of the Chianti VERSION.\n");
457 strcpy(chPath,chPathSave);
459 strcat( chPath,
"masterlist" );
462 strcpy( chPathSave , chPath );
475 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) == NULL )
477 fprintf(
ioQQQ,
" database_readin could not read first line of CloudyChianti.ini.\n");
481 bool lgEOL1=
true, lgEOL2=
true, lgEOL3=
true;
482 long int nMonRd=-1, nDayRd=-1;
484 long int nYrRd = (long)
FFmtRead(chLine,&ipST,
sizeof(chLine),&lgEOL1);
487 nMonRd = (long)
FFmtRead(chLine,&ipST,
sizeof(chLine),&lgEOL2);
489 nDayRd = (long)
FFmtRead(chLine,&ipST,
sizeof(chLine),&lgEOL3);
493 static long int nYr=11 , nMon = 10, nDay = 3;
496 fprintf(
ioQQQ,
"PROBLEM, there must be three magic numbers on the first line of the chianti masterlist file.\n");
498 " I expected to find the numbers %2.2li %2.2li %2.2li.\n" ,
503 if( ( nYrRd != nYr ) || ( nMonRd != nMon ) || ( nDayRd != nDay ) )
506 " database_readin: the Chianti masterlist file is not the current version.\n" );
508 " database_readin obtain the current version from the Cloudy web site.\n" );
510 " I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
511 nYr , nMon , nDay , nYrRd , nMonRd , nDayRd );
512 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
516 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) == NULL )
518 fprintf(
ioQQQ,
" database_readin could not read first line of CHIANTI masterlist.\n");
525 if ((chLine[0]!=
'#') && (chLine[0]!=
'\n')&&(chLine[0]!=
'\t')&&(chLine[0]!=
'\r'))
529 strcpy(chDLine, chLine);
530 chToken = strtok(chDLine,
" \t\n");
532 fixit(
"insert logic here to exclude some ions");
535 if( chToken[strlen(chToken)-1] !=
'd' )
538 ASSERT( nSpeciesCHIANTI + 1 <= MAX_NUM_SPECIES );
539 strcpy( chLabels[
nSpecies], chToken );
540 strcpy( chLabelsOrig[nSpecies], chLabels[nSpecies]);
544 char *chNumLevs = strtok(NULL,
"\n");
545 if( chNumLevs != NULL )
549 long numLevs = (long)
FFmtRead(chNumLevs,&i,
sizeof(chLine),&lgEOL);
559 fprintf(
ioQQQ,
"PROBLEM the limit to the number of levels must be positive, it was %li\n", numLevs);
566 bool skipSpecies =
false;
569 for(
int j = nSpeciesLAMDA; j < (nSpecies - nSpeciesCHIANTI); j++)
571 if( strcmp( chLabelsOrig[j], chLabelsOrig[nSpecies] ) == 0)
573 fprintf(
ioQQQ,
"Skipping the Chianti version of %s, using Stout version\n",chLabels[nSpecies]);
579 for(
int j = nSpecies - nSpeciesCHIANTI; j <
nSpecies; j++)
581 if( strcmp( chLabelsOrig[j], chLabelsOrig[nSpecies] ) == 0)
591 char *chElement, chTokenTemp[7];
592 strcpy( chTokenTemp, chToken );
593 (void) strtok(chTokenTemp,
" \n");
594 chElement = strtok(chTokenTemp,
"_");
599 strcpy( chPaths[nSpecies],
"chianti" );
601 strcat( chPaths[nSpecies], chElement );
603 strcat( chPaths[nSpecies], chLabels[nSpecies] );
605 strcat( chPaths[nSpecies], chLabels[nSpecies] );
607 ASSERT( isalpha(chToken[0]) );
610 if( isalpha(chToken[1]) )
621 ASSERT( chToken[cursor]==
'_' );
623 ASSERT( isdigit(chToken[cursor]) );
625 if( isdigit(chToken[cursor+1]) )
627 chLabels[
nSpecies][2] = chToken[cursor++];
628 chLabels[
nSpecies][3] = chToken[cursor++];
633 chLabels[
nSpecies][3] = chToken[cursor++];
636 ASSERT( chToken[cursor]==
'\0' || chToken[cursor]==
'd' );
645 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioMASTERLIST ) != NULL );
647 fclose(ioMASTERLIST);
677 size_t los = strlen(chLabels[i]);
686 if( numLevels[i] > 0 )
698 if( i<nSpeciesLAMDA )
703 else if( i < nSpeciesLAMDA + nSpeciesSTOUT )
708 else if( i < nSpeciesLAMDA + nSpeciesSTOUT + nSpeciesCHIANTI )
722 for(
long i=nSpeciesLAMDA; i<nSpeciesLAMDA+nSpeciesSTOUT; i++ )
727 for(
long i=nSpeciesLAMDA+nSpeciesSTOUT; i<nSpeciesLAMDA+nSpeciesSTOUT+nSpeciesCHIANTI; i++ )
737 fprintf(
save.
ipSDSFile,
"Chianti (C), Stout(S), Iso-sequences (I), old internal treatment( ).\n\n");
740 for(
int i=0; i<
LIMELM; i++)
746 for(
int i=0; i<
LIMELM; i++)
749 for(
int j=0; j<i+1; j++)
759 fprintf(
ioQQQ,
"\n\nDEBUG: Below are the contents of chdBaseSources[][]. It should contain a database name for each species.\n");
761 for(
int i=0; i<
LIMELM; i++)
766 for(
int i = 0; i <
LIMELM; i++ )
769 for(
int j = 0; j < LIMELM+1; j++ )
781 for( intNoSp=0; intNoSp<
nSpecies; intNoSp++ )
801 bool lgLevelsToTrim =
true;
804 string speciesLabel =
dBaseStates[ipSpecies].chLabel();
806 long totalNumLevels =
dBaseSpecies[ipSpecies].numLevels_max;
808 while( lgLevelsToTrim )
811 lgLevelsToTrim =
true;
815 fprintf(
ioQQQ,
"PROBLEM: Spectrum %s (species: %s) has no transition probabilities out of the first %li levels.\n",
816 spectralLabel, speciesLabel.c_str(), totalNumLevels);
817 fprintf(
ioQQQ,
"Consider allowing Cloudy to use more levels (see Hazy 1 SPECIES STOUT/CHIANTI LEVELS MAX), add more low-level"
818 " transition probabilities, or disable %s in the masterlist.\n\n", spectralLabel);
822 for(
int ipLo = 0; ipLo < ipHi; ipLo++)
825 aul = tr->Emis().Aul();
828 fprintf(
ioQQQ,
"trim_levels():\t%s\t%i\t%li\t%e\n", spectralLabel, ipLo+1, ipHi+1, aul);
833 lgLevelsToTrim =
false;
844 fprintf(
ioQQQ,
"Spectrum %s (species: %s) trimmed to %li levels (original %li) to have positive Aul.\n",
846 speciesLabel.c_str(),
862 strncpy( chToken, sp->
chLabel, 2 );
864 if( strcmp(
"p-", chToken )==0 )
866 else if( strcmp(
"o-", chToken )==0 )
868 else if( strcmp(
"e-", chToken )==0 )
870 else if( strcmp(
"a-", chToken )==0 )
875 fixit(
"what fraction should e-type and a-type Methanol have? Assume 50/50 for now.");
904 if ( strlen(chLabel) != 4 ||
905 ! (isalpha(chLabel[0])) ||
906 ! (chLabel[1] ==
' ' || isalpha(chLabel[1])) ||
907 ! (chLabel[2] ==
' ' || isdigit(chLabel[2])) ||
908 ! (chLabel[3] ==
' ' || isdigit(chLabel[3])))
914 strncpy( chToken, chLabel, 2 );
916 for(
long ipElement=0; ipElement<
LIMELM; ipElement++ )
924 strncpy( chToken, chLabel + 2, 2 );
925 IonStg = atoi(chToken);
936 char chStage[5] = {
'\0'};
940 sprintf( chStage,
"+%li", ion );
942 return chLabelChemical + chStage;
949 strncpy( chLabelChemical, chemLab.c_str(), size_t(
CHARS_SPECIES) );
958 ASSERT( IonStg >= 1 && IonStg <= nelem+2 );
961 if( nelem - (IonStg-1) <
NISO )
964 fprintf(
ioQQQ,
" Iso-sequences are handled by our own model.\n");
985 size_t plus_sign_pos = chLabelChem.find_first_of(
'+' );
987 if( plus_sign_pos == string::npos )
990 chLabelSpec = chLabelChem;
992 if( chLabelSpec.length() == 1 )
995 if( chLabelSpec.length() == 2 &&
1004 string elm = chLabelChem.substr( 0, plus_sign_pos );
1006 if( elm.length() == 1 )
1013 chLabelSpec = chLabelChem;
1018 int ionstg = atoi( chLabelChem.substr( plus_sign_pos+1 ).c_str() );
1026 chLabelSpec += tmp.str();
1038 long nelem = 0, IonStg;
1043 fixit(
"should never be used if lgMolecular");
1098 fprintf(
ioQQQ,
" PROBLEM: could not find species %li - %s\n",i,
1125 printf(
"The species is %s \n",
dBaseSpecies[i].chLabel);
1126 printf(
"The data output is in the following format \n");
1127 printf(
"Label Energy St.wt Pop Lifetime\n");
1131 printf(
"This is the %ld state \n",j);
1132 printf(
"%s %f %f %f %e \n",
dBaseStates[i][j].chLabel().c_str(),
1150 em !=
dBaseTrans[intSpIndex].Emis().end(); ++em)
1152 fsumAs[(*em).Tran().ipHi()] += (*em).Aul();
1158 (*em).iRedisFun() =
ipPRD;
1163 if( em->Tran().ipLo() == 0 )
1166 em->iRedisFun() =
ipCRD;
1172 em->iRedisFun() =
ipCRDW;
1178 for(
int ipHi=1; ipHi <
dBaseSpecies[intSpIndex].numLevels_max; ipHi++ )
1180 dBaseStates[intSpIndex][ipHi].lifetime() = 1./fsumAs[ipHi];
char chLamdaFile[FILENAME_PATH_LENGTH]
STATIC void database_prep(int)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
vector< StoutCollArray > StoutCollData
const int FILENAME_PATH_LENGTH_2
char chStoutFile[FILENAME_PATH_LENGTH]
NORETURN void TotalInsanity(void)
STATIC void set_fractionation(species *sp)
vector< multi_arr< int, 2 > > ipdBaseTrans
void atmdat_STOUT_readin(long intNS, char *chFileName)
multi_arr< CollRateCoeffArray, 2 > AtmolCollRateCoeff
void trimTrailingWhiteSpace(string &str)
STATIC void states_propprint(void)
char chVersion[iVersionLength]
static t_version & Inst()
STATIC void spectral_to_chemical(char *chLabelChemical, char *chLabel, long &nelem, long &IonStg)
t_elementnames elementnames
bool lgIonStoutOn[LIMELM][LIMELM+1]
void uncaps(char *chCard)
bool lgIonChiantiOn[LIMELM][LIMELM+1]
bool lgdBaseSourceExists[LIMELM][LIMELM+1]
static const double aulThreshold
double energy(const genericState &gs)
char chdBaseSources[LIMELM][LIMELM+1][10]
void chemical_to_spectral(const string chLabelChem, string &chLabelSpec)
char chCloudyChiantiFile[FILENAME_PATH_LENGTH]
molecule * findspecies(const char buf[])
valarray< class molezone > species
char species[CHARS_SPECIES]
void parsespect(char *chLabel, long &nelem, long &IonStg)
STATIC void states_nelemfill(void)
vector< multi_arr< CollSplinesArray, 3 > > AtmolCollSplines
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
realnum AtomicWeight[LIMELM]
string makeChemical(long nelem, long ion)
STATIC void trim_levels(long)
vector< vector< long > > ipSpecIon
void atmdat_LAMDA_readin(long intNS, char *chFileName)
#define DEBUG_ENTRY(funcname)
STATIC void states_popfill(void)
vector< qList > dBaseStates
vector< species > dBaseSpecies
void database_readin(void)
int fprintf(const Output &stream, const char *format,...)
void atmdat_CHIANTI_readin(long intNS, char *chFileName)
double maxWN[LIMELM][LIMELM+1]
vector< TransitionList > AllTransitions
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
bool isElementSym(const char *chSym)
vector< TransitionList > dBaseTrans
static const int iVersionLength
void check_data(const char *fname, access_scheme scheme)
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)