00001
00002
00003 #include "cddefines.h"
00004 #include "lines_service.h"
00005 #include "taulines.h"
00006 #include "trace.h"
00007 #include "string.h"
00008 #include "input.h"
00009 #include "thirdparty.h"
00010 #include "dense.h"
00011 #include "atmdat.h"
00012 #include "mole.h"
00013 #include "elementnames.h"
00014 #include "version.h"
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 void states_popfill( void);
00034 void states_nelemfill(void);
00035 void database_prep(int);
00036 void emislines_fillredis(void);
00037 void set_fractionation( species *sp );
00038
00039 STATIC void states_propprint(void);
00040
00041 STATIC void SpeciesJunk( species *sp );
00042
00043 #define DEBUGSTATE false
00044
00045 emission *AddLine2Stack( realnum Aul, transition *trans )
00046 {
00047
00048 DEBUG_ENTRY( "AddLine2Stack()" );
00049
00050 if(linesAdded2 >= MAX_NUM_LINES)
00051 {
00052 fprintf(ioQQQ,"DISASTER 1 the number of lines linesAdded2=%li is greater than"\
00053 " the largest number allowed, MAX_NUM_LINES=%li\n Consider increasing "\
00054 "MAX_NUM_LINES\n", linesAdded2 , MAX_NUM_LINES );
00055 cdEXIT(EXIT_FAILURE);
00056 }
00057 dBaseLines[linesAdded2].Aul = Aul;
00058 dBaseLines[linesAdded2].tran = trans;
00059 linesAdded2++;
00060 if(linesAdded2 >= MAX_NUM_LINES)
00061 {
00062 fprintf(ioQQQ,"DISASTER 2 the number of lines linesAdded2=%li is greater than"\
00063 " the largest number allowed, MAX_NUM_LINES=%li\n Consider increasing "\
00064 "MAX_NUM_LINES\n", linesAdded2 , MAX_NUM_LINES );
00065 cdEXIT(EXIT_FAILURE);
00066 }
00067 return &dBaseLines[linesAdded2-1];
00068 }
00069
00070 void database_readin( void )
00071 {
00072 int i,j,intNoSp;
00073
00074 FILE *ioMASTERLIST, *ioVERSION;
00075
00076 char *chToken;
00077
00078 char chLine[FILENAME_PATH_LENGTH_2],
00079 chDLine[FILENAME_PATH_LENGTH_2],
00080 chPath[FILENAME_PATH_LENGTH_2] = "";
00081
00082 const int MAX_NUM_SPECIES = 1000;
00083
00084 fixit();
00085 char chLabels[MAX_NUM_SPECIES][7];
00086 char chPaths[MAX_NUM_SPECIES][FILENAME_PATH_LENGTH_2];
00087
00088 static int nCalled = 0;
00089 long nSpeciesLAMDA, nSpeciesCHIANTI;
00090
00091 DEBUG_ENTRY( "database_readin()" );
00092
00093 linesAdded2 = 0;
00094
00095
00096 if( nCalled > 0 )
00097 {
00098 return;
00099 }
00100
00101
00102 ++nCalled;
00103
00104
00105 nSpecies = 0;
00106
00108
00109
00110
00112
00113
00114
00115 nSpeciesLAMDA = 0;
00116
00117 if( atmdat.lgLamdaOn )
00118 {
00119 long numModelsNotUsed = 0;
00120 strcpy( chPath, "lamda" );
00121 strcat( chPath, input.chDelimiter );
00122 strcat( chPath, "masterlist" );
00123 ioMASTERLIST = open_data( chPath, "r" );
00124
00125 if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
00126 {
00127 fprintf( ioQQQ, " database_readin could not read first line of LAMDA masterlist.\n");
00128 cdEXIT(EXIT_FAILURE);
00129 }
00130
00131 do
00132 {
00133 if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
00134 {
00135 strcpy(chDLine, chLine);
00136 chToken = strtok(chDLine," \t\n");
00137 if( findspecies( chToken ) != &null_mole ||
00138 ( chToken[1]=='-' && findspecies( chToken+2 ) != &null_mole ) )
00139 {
00140 ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
00141 ASSERT( nSpeciesLAMDA + 1 <= MAX_NUM_SPECIES );
00142 ASSERT( strlen(chToken) < CHARS_SPECIES );
00143 strcpy( chLabels[nSpecies], chToken );
00144
00145
00146 strcpy( chPaths[nSpecies], "lamda" );
00147 strcat( chPaths[nSpecies], input.chDelimiter );
00148 chToken = strtok( NULL," \t\n" );
00149 strcat( chPaths[nSpecies], chToken );
00150 ++nSpecies;
00151 ++nSpeciesLAMDA;
00152 }
00153 else
00154 ++numModelsNotUsed;
00155 }
00156 }
00157 while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
00158
00159
00160
00161
00162
00163 fclose(ioMASTERLIST);
00164 }
00165
00167
00168
00169
00171
00172 nSpeciesCHIANTI = 0;
00173
00174 if( atmdat.lgChiantiOn )
00175 {
00176 char chPathSave[FILENAME_PATH_LENGTH_2];
00177 strcpy( chPath, "chianti" );
00178 strcat( chPath, input.chDelimiter );
00179
00180
00181 strcpy( chPathSave , chPath );
00182 strcat(chPath,"VERSION");
00183 ioVERSION = open_data(chPath,"r");
00184 if( read_whole_line( chLine , (int)sizeof(chLine) , ioVERSION ) == NULL )
00185 {
00186 fprintf( ioQQQ, " database_readin could not read first line of the Chianti VERSION.\n");
00187 cdEXIT(EXIT_FAILURE);
00188 }
00189 fclose(ioVERSION);
00190
00191 strncpy(atmdat.chVersion,chLine,atmdat.iVersionLength);
00192
00193 atmdat.chVersion[atmdat.iVersionLength-1] = '\0';
00194
00195 strcpy(chPath,chPathSave);
00196
00197 strcat( chPath, "masterlist" );
00198 strcat( chPath, input.chDelimiter );
00199
00200 strcpy( chPathSave , chPath );
00201
00202
00203 strcat( chPath, "CloudyChianti.ini" );
00204
00205
00206 if( (ioMASTERLIST = open_data( chPath, "r" ,AS_DATA_ONLY_TRY )) ==NULL )
00207 {
00208 strcat( chPathSave, "masterlist.ions.v6" );
00209 ioMASTERLIST = open_data( chPathSave, "r" );
00210 }
00211 else
00212 {
00213
00214 if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
00215 {
00216 fprintf( ioQQQ, " database_readin could not read first line of CloudyChianti.ini.\n");
00217 cdEXIT(EXIT_FAILURE);
00218 }
00219
00220 bool lgEOL;
00221 long int ip = 1;
00222 long int nYrRd = (long)FFmtRead(chLine,&ip,INPUT_LINE_LENGTH,&lgEOL);
00223 long int nMonRd = (long)FFmtRead(chLine,&ip,INPUT_LINE_LENGTH,&lgEOL);
00224 long int nDayRd = (long)FFmtRead(chLine,&ip,INPUT_LINE_LENGTH,&lgEOL);
00225
00226 static long int nYr=10 , nMon = 7, nDay = 13;
00227 if( ( nYrRd != nYr ) || ( nMonRd != nMon ) || ( nDayRd != nDay ) )
00228 {
00229 fprintf( ioQQQ,
00230 " database_readin: the version of CloudyChianti.ini is not the current version.\n" );
00231 fprintf( ioQQQ,
00232 " database_readin obtain the current version from the Cloudy web site.\n" );
00233 fprintf( ioQQQ,
00234 " I expected to find the number %2.2li %2.2li %2.2li and got %2.2li %2.2li %2.2li instead.\n" ,
00235 nYr , nMon , nDay , nYrRd , nMonRd , nDayRd );
00236 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00237 cdEXIT(EXIT_FAILURE);
00238 }
00239 }
00240
00241 if( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) == NULL )
00242 {
00243 fprintf( ioQQQ, " database_readin could not read first line of CHIANTI masterlist.\n");
00244 cdEXIT(EXIT_FAILURE);
00245 }
00246
00247 do
00248 {
00249
00250 if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r'))
00251 {
00252 strcpy(chDLine, chLine);
00253 chToken = strtok(chDLine," \n");
00254
00255 fixit();
00256
00257
00258 if( chToken[3]!='d' && chToken[4]!='d' && chToken[5]!='d' )
00259 {
00260 ASSERT( nSpecies + 1 <= MAX_NUM_SPECIES );
00261 ASSERT( nSpeciesCHIANTI + 1 <= MAX_NUM_SPECIES );
00262 strcpy( chLabels[nSpecies], chToken );
00263
00264 char *chElement, chTokenTemp[7];
00265 strcpy( chTokenTemp, chToken );
00266 chElement = strtok(chTokenTemp," \n");
00267 chElement = strtok(chTokenTemp,"_");
00268 uncaps( chElement );
00269
00270
00271
00272 strcpy( chPaths[nSpecies], "chianti" );
00273 strcat( chPaths[nSpecies], input.chDelimiter );
00274 strcat( chPaths[nSpecies], chElement );
00275 strcat( chPaths[nSpecies], input.chDelimiter );
00276 strcat( chPaths[nSpecies], chLabels[nSpecies] );
00277 strcat( chPaths[nSpecies], input.chDelimiter );
00278 strcat( chPaths[nSpecies], chLabels[nSpecies] );
00279
00280 ASSERT( isalpha(chToken[0]) );
00281 long cursor=0;
00282 chLabels[nSpecies][0] = chToken[0];
00283 if( isalpha(chToken[1]) )
00284 {
00285 chLabels[nSpecies][1] = chToken[1];
00286 cursor = 2;
00287 }
00288 else
00289 {
00290 chLabels[nSpecies][1] = ' ';
00291 cursor = 1;
00292 }
00293
00294 ASSERT( chToken[cursor++]=='_' );
00295 ASSERT( isdigit(chToken[cursor]) );
00296
00297 if( isdigit(chToken[cursor+1]) )
00298 {
00299 chLabels[nSpecies][2] = chToken[cursor++];
00300 chLabels[nSpecies][3] = chToken[cursor++];
00301 }
00302 else
00303 {
00304 chLabels[nSpecies][2] = ' ';
00305 chLabels[nSpecies][3] = chToken[cursor++];
00306 }
00307 chLabels[nSpecies][4] = '\0';
00308 ASSERT( chToken[cursor]=='\0' || chToken[cursor]=='d' );
00309
00310
00311 chLabels[nSpecies][0] = toupper( chLabels[nSpecies][0] );
00312 fprintf( ioQQQ," Using CHIANTI model %7s.\n", chLabels[nSpecies] );
00313 ++nSpecies;
00314 ++nSpeciesCHIANTI;
00315 }
00316 else
00317 {
00318 fprintf( ioQQQ," NOT Using CHIANTI model %7s.\n", chToken );
00319 }
00320 }
00321 }
00322 while( read_whole_line( chLine , (int)sizeof(chLine) , ioMASTERLIST ) != NULL );
00323
00324 fclose(ioMASTERLIST);
00325 }
00326
00327
00328 if( nSpecies==0 )
00329 return;
00330
00331
00332 Species = (species *)MALLOC( (unsigned long)nSpecies*sizeof(species));
00333
00334
00335
00336 CollRatesArray = (double****)MALLOC((unsigned long)nSpecies * sizeof(double***));
00337 AtmolCollRateCoeff = (CollRateCoeffArray **)MALLOC((unsigned long)nSpecies * sizeof(CollRateCoeffArray*));
00338 AtmolCollSplines = (CollSplinesArray****)MALLOC((unsigned long)nSpecies *sizeof(CollSplinesArray***));
00339
00340
00341 for( i=0; i<nSpecies; i++ )
00342 {
00343 AtmolCollRateCoeff[i] = (CollRateCoeffArray *)MALLOC(NUM_COLLIDERS * sizeof(CollRateCoeffArray));
00344 CollRatesArray[i] = (double***)MALLOC(NUM_COLLIDERS * sizeof(double**));
00345 }
00346
00347
00348 for( i=0; i<nSpecies; i++ )
00349 {
00350 for( j=0; j<NUM_COLLIDERS; j++ )
00351 {
00352 AtmolCollRateCoeff[i][j].ntemps = 0;
00353 AtmolCollRateCoeff[i][j].temps = NULL;
00354 AtmolCollRateCoeff[i][j].collrates = NULL;
00355 }
00356 }
00357
00358
00359 dBaseStates = (quantumState **)MALLOC((unsigned long)nSpecies * sizeof(quantumState*));
00360 dBaseTrans = (transition ***)MALLOC((unsigned long)nSpecies * sizeof(transition**));
00361
00362 for( i = 0; i < nSpecies; i++ )
00363 {
00364 SpeciesJunk( &Species[i] );
00365
00366 size_t los = max(4,strlen(chLabels[i]));
00367 ASSERT( los >= 4 && los <= 7 );
00368 Species[i].chLabel = new char[los+1];
00369 strcpy(Species[i].chLabel,chLabels[i]);
00370
00371
00372 set_fractionation( &Species[i] );
00373
00374
00375 if( Species[i].chLabel[2]=='\0' )
00376 {
00377 Species[i].chLabel[2]=' ';
00378 Species[i].chLabel[3]=' ';
00379 Species[i].chLabel[4]='\0';
00380 }
00381 else if( Species[i].chLabel[3]=='\0' )
00382 {
00383 Species[i].chLabel[3]=' ';
00384 Species[i].chLabel[4]='\0';
00385 }
00386
00387 if( i<nSpeciesLAMDA )
00388 {
00389
00390 atmdat_LAMDA_readin( i, chPaths[i] );
00391 }
00392 else if( i < nSpeciesLAMDA + nSpeciesCHIANTI )
00393 {
00394
00395 atmdat_CHIANTI_readin( i, chPaths[i] );
00396 }
00397 else
00398 TotalInsanity();
00399 }
00400
00401
00402
00403 states_popfill();
00404
00405 states_nelemfill();
00406
00407
00408
00409 for( intNoSp=0; intNoSp<nSpecies; intNoSp++ )
00410 {
00411 database_prep(intNoSp);
00412 }
00413
00414
00415 emislines_fillredis();
00416
00417
00418 if(DEBUGSTATE)
00419 states_propprint();
00420 return;
00421 }
00422
00423 void set_fractionation( species *sp )
00424 {
00425 DEBUG_ENTRY("set_fractionation()");
00426
00427 char chToken[3];
00428
00429 sp->fracIsotopologue = 1.f;
00430
00431 strncpy( chToken, sp->chLabel, 2 );
00432 chToken[2] = '\0';
00433 if( strcmp( "p-", chToken )==0 )
00434 sp->fracType = 0.25f;
00435 else if( strcmp( "o-", chToken )==0 )
00436 sp->fracType = 0.75f;
00437 else if( strcmp( "e-", chToken )==0 )
00438 sp->fracType = 0.5f;
00439 else if( strcmp( "a-", chToken )==0 )
00440 sp->fracType = 0.5f;
00441 else
00442 sp->fracType = 1.0f;
00443
00444 fixit();
00445
00446
00447 if( sp->chLabel[1]=='-')
00448 memmove(sp->chLabel,sp->chLabel+2,strlen(sp->chLabel+2)+1);
00449
00450 return;
00451 }
00452
00453
00454 void emislines_fillredis(void)
00455 {
00456 DEBUG_ENTRY("emislines_fillredis()");
00457 for( long i=0; i<linesAdded2; i++ )
00458 {
00459 dBaseLines[i].iRedisFun = ipPRD;
00460 }
00461 return;
00462 }
00463
00464
00465 void states_nelemfill(void)
00466 {
00467 DEBUG_ENTRY( "states_nelemfill()" );
00468
00469 for( long i=0; i<nSpecies; i++ )
00470 {
00471 long nelem = 0, IonStg;
00472
00473 if( Species[i].lgMolecular )
00474 {
00475 fixit();
00476
00477
00478 nelem = -1;
00479 IonStg = -1;
00480 }
00481 else
00482 {
00483 char chToken[3];
00484 strncpy( chToken, Species[i].chLabel, 2 );
00485 chToken[2] = '\0';
00486 for( long ipElement=0; ipElement<LIMELM; ipElement++ )
00487 {
00488 if( strcmp( elementnames.chElementSym[ipElement], chToken )==0 )
00489 {
00490 nelem = ipElement + 1;
00491 break;
00492 }
00493 }
00494 ASSERT( nelem > 0 && nelem <= LIMELM );
00495 strncpy( chToken, Species[i].chLabel + 2, 2 );
00496 IonStg = atoi(chToken);
00497 ASSERT( IonStg >= 1 && IonStg <= nelem+1 );
00498 Species[i].fmolweight = dense.AtomicWeight[nelem-1];
00499
00500 dense.lgIonChiantiOn[nelem-1][IonStg-1] = true;
00501 }
00502
00503 for( long j=0; j<Species[i].numLevels_max; j++ )
00504 {
00505 dBaseStates[i][j].nelem = nelem;
00506 dBaseStates[i][j].IonStg = IonStg;
00507 }
00508 }
00509 return;
00510 }
00511
00512
00513 STATIC void states_propprint(void)
00514 {
00515 DEBUG_ENTRY( "states_propprint()" );
00516
00517 for( long i=0; i<nSpecies; i++ )
00518 {
00519 printf("The species is %s \n",Species[i].chLabel);
00520 printf("The data output is in the following format \n");
00521 printf("Label Energy St.wt Pop Lifetime\n");
00522
00523 for( long j=0; j<Species[i].numLevels_max; j++ )
00524 {
00525 printf("This is the %ld state \n",j);
00526 printf("%s %f %f %f %e \n",dBaseStates[i][j].chLabel,
00527 dBaseStates[i][j].energy.WN(),
00528 dBaseStates[i][j].g,
00529 dBaseStates[i][j].Pop,
00530 dBaseStates[i][j].lifetime);
00531 }
00532 }
00533 return;
00534 }
00535
00536
00537 void database_prep(int intSpIndex)
00538 {
00539 realnum fsumAs;
00540 int ipHi,ipLo;
00541
00542 DEBUG_ENTRY( "database_prep()" );
00543
00544
00545 dBaseStates[intSpIndex][0].lifetime= BIGFLOAT;
00546 for( ipHi=1; ipHi < Species[intSpIndex].numLevels_max; ipHi++ )
00547 {
00548 fsumAs = SMALLFLOAT;
00549 for( ipLo=0; ipLo<ipHi; ipLo++ )
00550 {
00551 if(dBaseTrans[intSpIndex][ipHi][ipLo].Emis!=NULL)
00552 fsumAs = fsumAs + dBaseTrans[intSpIndex][ipHi][ipLo].Emis->Aul;
00553 }
00554 dBaseStates[intSpIndex][ipHi].lifetime = 1./fsumAs;
00555 }
00556 return;
00557 }
00558
00559
00560 STATIC void SpeciesJunk( species *sp )
00561 {
00562 sp->chLabel = NULL;
00563 set_NaN(sp->fmolweight);
00564 set_NaN(sp->fracIsotopologue );
00565 set_NaN(sp->fracType);
00566 sp->lgMolecular = false;
00567 sp->numLevels_local = -INT_MAX;
00568 sp->numLevels_max = -INT_MAX;
00569 sp->lgActive = false;
00570
00571 return;
00572 }