cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "lines_service.h" 00005 #include "taulines.h" 00006 #include "trace.h" 00007 #include "string.h" 00008 #include "thirdparty.h" 00009 #include "dense.h" 00010 #include "atmdat.h" 00011 /*File nemala.cpp was developed by Humeshkar B Nemala as a part of his thesis work during the Summer of 2007*/ 00012 /* Initially the code has been developed to read in energy levels,radiative and 00013 * collisional data from the CHIANTI and LEIDEN databases. The idea is to extend it to more databases. 00014 * In the case of the Leiden database there is a single .dat file which has the energy levels information, 00015 * radiative and collisional data, with the data corresponding to each collider coming one after the other. 00016 * In the case of CHIANTI, the energy levels data, radiative data and collision data are present in seperate files. 00017 * While LEIDEN gives collisional rate coefficients, CHIANTI gives collisional strengths. 00018 * In the case of CHIANTI only two colliders are used:electrons and protons. They appear as separate files. 00019 * The electron collision strengths files are always expected to be there. A flag is set and data processed 00020 * if the file on proton collision strengths is available.*/ 00021 00022 /* There is an initialization file called species.ini which tells Cloudy what kind of data is to be used */ 00023 /* Structures are created separately to hold the transition data,radiative and collisional data */ 00024 /* The collisional structures are different for different databases depending upon whether */ 00025 /* collisional strengths or collisional rate coefficients are used.Finally a superstructure is constructed to hold */ 00026 /* the total collisional rate obtained by considering all the colliders */ 00027 /* The colliders considered are electron,proton,Atomic Hydrogen,He,He+,He++,Ortho Molecular Hydrogen,Para Molecular Hydrogen and Molecular Hydrogen */ 00028 void database_readin(void); 00029 void states_popfill( void); 00030 void states_nelemfill(void); 00031 void database_prep(int); 00032 void emislines_fillredis(void); 00033 STATIC void states_propprint(void); 00034 STATIC int getAtNo(char []); 00035 00036 #define DEBUGSTATE false 00037 00038 emission *AddLine2Stack( realnum Aul, transition *trans ) 00039 { 00040 00041 DEBUG_ENTRY( "AddLine2Stack()" ); 00042 00043 ASSERT(linesAdded2 < MAX_NUM_LINES); 00044 atmolEmis[linesAdded2].Aul = Aul; 00045 atmolEmis[linesAdded2].tran = trans; 00046 linesAdded2++; 00047 ASSERT(linesAdded2 <= MAX_NUM_LINES); 00048 return &atmolEmis[linesAdded2-1]; 00049 } 00050 00051 void Nemala_Start( void ) 00052 { 00053 int i,j,los,intNoSp; 00054 00055 FILE *atmolDATA; 00056 char *chToken,*chAtNo,*chIonStg; 00057 00058 char chLine[FILENAME_PATH_LENGTH_2], 00059 chDLine[FILENAME_PATH_LENGTH_2]; 00060 00061 static int nCalled = 0; 00062 00063 DEBUG_ENTRY( "Nemala_Start()" ); 00064 00065 linesAdded2 = 0; 00066 00067 /* only do this once. */ 00068 if( nCalled > 0 ) 00069 { 00070 return; 00071 } 00072 00073 /* this is first call, increment the nCalled counterso never do this again */ 00074 ++nCalled; 00075 00076 // Reading in the species.ini file: Humeshkar B Nemala 00077 if( trace.lgTrace ) 00078 fprintf( ioQQQ," Nemala_Start opening species.ini:"); 00079 00080 atmolDATA = open_data( "species.ini", "r" ); 00081 00082 if( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) == NULL ) 00083 { 00084 fprintf( ioQQQ, " Nemala_Start could not read first line of species.ini.\n"); 00085 cdEXIT(EXIT_FAILURE); 00086 } 00087 00088 /* count how many lines are in the file, ignoring all lines 00089 * starting with '#':This would give the number of molecules */ 00090 nSpecies = 0; 00091 while( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) != NULL ) 00092 { 00093 /* we want to count the lines that do not start with % 00094 * since these contain data */ 00095 if( (chLine[0] != '#')&&(chLine[0] != '\n')&&(chLine[0] != ' ')) 00096 ++nSpecies; 00097 } 00098 if(DEBUGSTATE) 00099 printf("The number of species is %li \n",nSpecies); 00100 00101 /* now rewind the file so we can read it a second time*/ 00102 if( fseek( atmolDATA , 0 , SEEK_SET ) != 0 ) 00103 { 00104 fprintf( ioQQQ, " Nemala_Start could not rewind species.ini.\n"); 00105 cdEXIT(EXIT_FAILURE); 00106 } 00107 00108 /*Once we know the number of species we initialize a pointer to a one 00109 * dimensional array which contains information about the DB Type */ 00110 lgSpeciesMolecule = (bool *)MALLOC( (unsigned long)(nSpecies)*sizeof(bool)); 00111 00112 /*Initialization of the Species Structure*/ 00113 Species = (species *)MALLOC( (unsigned long)nSpecies*sizeof(species)); 00114 00115 /*Initialization of the collisional rates array structure*/ 00116 /*Assume that there are 9 colliders*/ 00117 CollRatesArray = (double****)MALLOC((unsigned long)nSpecies * sizeof(double***)); 00118 AtmolCollRateCoeff = (CollRateCoeffArray **)MALLOC((unsigned long)nSpecies * sizeof(CollRateCoeffArray*)); 00119 AtmolCollSplines = (CollSplinesArray****)MALLOC((unsigned long)nSpecies *sizeof(CollSplinesArray***)); 00120 00121 /*Mallocing here takes care of the number of colliders*/ 00122 for( i=0; i<nSpecies; i++ ) 00123 { 00124 AtmolCollRateCoeff[i] = (CollRateCoeffArray *)MALLOC(NUM_COLLIDERS * sizeof(CollRateCoeffArray)); 00125 //AtmolCollSplines[i] = (CollSplinesArray***)MALLOC(NUM_COLLIDERS *sizeof(CollSplinesArray**)); 00126 CollRatesArray[i] = (double***)MALLOC(NUM_COLLIDERS * sizeof(double**)); 00127 } 00128 /*Initializing all the elements of the matrix to 0*/ 00129 /*For each species and collider*/ 00130 for( i=0; i<nSpecies; i++ ) 00131 { 00132 for( j=0; j<NUM_COLLIDERS; j++ ) 00133 { 00134 AtmolCollRateCoeff[i][j].ntemps = 0; 00135 AtmolCollRateCoeff[i][j].temps = NULL; 00136 AtmolCollRateCoeff[i][j].collrates = NULL; 00137 } 00138 } 00139 00140 /*****************************************/ 00141 i=0; 00142 while( read_whole_line( chLine , (int)sizeof(chLine) , atmolDATA ) != NULL ) 00143 { 00144 caps( chLine ); 00145 if ((chLine[0]!='#') && (chLine[0]!='\n')&&(chLine[0]!='\t')&&(chLine[0]!='\r')) 00146 { 00147 ASSERT(i<nSpecies); 00148 if(DEBUGSTATE) 00149 printf("%s",chLine); 00150 strcpy(chDLine, chLine); 00151 chToken = strtok(chDLine," "); 00152 los = (int)strlen(chToken); 00153 if(DEBUGSTATE) 00154 printf("The length of the string is %d\n",los); 00155 Species[i].chptrSpName = (char *)MALLOC((unsigned long)(los+1) * sizeof(char)); 00156 strcpy(Species[i].chptrSpName,chToken); 00157 Species[i].intAtNo = 0; 00158 Species[i].intIonStage = 0; 00159 if(DEBUGSTATE) 00160 printf("The name of the species is%s \n",Species[i].chptrSpName); 00161 if( nMatch("LEID",chLine) ) 00162 { 00163 lgSpeciesMolecule[i] = true; 00164 } 00165 else if( nMatch("CHIA",chLine)) 00166 { 00167 lgSpeciesMolecule[i] = false; 00168 chAtNo = strtok(chToken,"_"); 00169 Species[i].intAtNo = getAtNo(chAtNo); 00170 chIonStg = strtok(NULL,"_"); 00171 Species[i].intIonStage = atoi(chIonStg)-1; 00172 } 00173 else 00174 TotalInsanity(); 00175 00176 if(DEBUGSTATE) 00177 printf("The species is molecular? %c \n",TorF(lgSpeciesMolecule[i]) ); 00178 ++i; 00179 } 00180 } 00181 00182 database_readin(); 00183 00184 /*Setting the population of states to 0*/ 00185 states_popfill(); 00186 /*Setting the values of nelem arbitrarily to 1*/ 00187 states_nelemfill(); 00188 00189 /*Setting nelem of the states to an arbitrary value*/ 00190 /*Loop over species*/ 00191 for( intNoSp=0; intNoSp<nSpecies; intNoSp++ ) 00192 { 00193 database_prep(intNoSp); 00194 } 00195 00196 /*Looping over the emission lines and filling in the redistribution function value*/ 00197 emislines_fillredis(); 00198 00199 /*To print the states*/ 00200 if(DEBUGSTATE) 00201 states_propprint(); 00202 return; 00203 } 00204 00205 /*Filling the redistribution function value*/ 00206 void emislines_fillredis(void) 00207 { 00208 DEBUG_ENTRY("emislines_fillredis()"); 00209 for( long i=0; i<linesAdded2; i++ ) 00210 { 00211 atmolEmis[i].iRedisFun = ipPRD; 00212 } 00213 return; 00214 } 00215 00216 /*This function fills the nelem value arbitrarily*/ 00217 void states_nelemfill(void) 00218 { 00219 DEBUG_ENTRY( "states_nelemfill()" ); 00220 00221 for( long i=0; i<nSpecies; i++ ) 00222 { 00223 for( long j=0; j<Species[i].numLevels_max; j++ ) 00224 { 00225 if( lgSpeciesMolecule[i] ) 00226 { 00227 fixit(); /* this is needed because right now, all lines test on 00228 * dense.xIonDense[nelem][IonStg] */ 00229 atmolStates[i][j].nelem = 1; 00230 atmolStates[i][j].IonStg = 1; 00231 } 00232 else 00233 { 00234 atmolStates[i][j].nelem = Species[i].intAtNo + 1; 00235 atmolStates[i][j].IonStg = Species[i].intIonStage; 00236 } 00237 } 00238 } 00239 return; 00240 } 00241 00242 /*This function prints the various properties of states*/ 00243 STATIC void states_propprint(void) 00244 { 00245 DEBUG_ENTRY( "states_propprint()" ); 00246 00247 for( long i=0; i<nSpecies; i++ ) 00248 { 00249 printf("The species is %s \n",Species[i].chptrSpName); 00250 printf("The data output is in the following format \n"); 00251 printf("Label Energy St.wt Pop Lifetime\n"); 00252 00253 for( long j=0; j<Species[i].numLevels_max; j++ ) 00254 { 00255 printf("This is the %ld state \n",j); 00256 printf("%s %f %f %f %e \n",atmolStates[i][j].chLabel, 00257 atmolStates[i][j].energy, 00258 atmolStates[i][j].g, 00259 atmolStates[i][j].Pop, 00260 atmolStates[i][j].lifetime); 00261 } 00262 } 00263 return; 00264 } 00265 00266 00267 void database_prep(int intSpIndex) 00268 { 00269 realnum fsumAs; 00270 int ipHi,ipLo; 00271 00272 DEBUG_ENTRY( "database_prep()" ); 00273 00274 /*Get the lifetimes*/ 00275 atmolStates[intSpIndex][0].lifetime= BIGFLOAT; 00276 for( ipHi=1; ipHi < Species[intSpIndex].numLevels_max; ipHi++ ) 00277 { 00278 fsumAs = SMALLFLOAT; 00279 for( ipLo=0; ipLo<ipHi; ipLo++ ) 00280 { 00281 if(atmolTrans[intSpIndex][ipHi][ipLo].Emis!=NULL) 00282 fsumAs = fsumAs + atmolTrans[intSpIndex][ipHi][ipLo].Emis->Aul; 00283 } 00284 atmolStates[intSpIndex][ipHi].lifetime = 1./fsumAs; 00285 } 00286 return; 00287 } 00288 00289 /*Separate Routine to read in the molecules*/ 00290 void database_readin(void) 00291 { 00292 DEBUG_ENTRY( "database_readin()" ); 00293 00294 long intNS; 00295 00296 atmolStates = (quantumState **)MALLOC((unsigned long)nSpecies * sizeof(quantumState*)); 00297 atmolTrans = (transition ***)MALLOC((unsigned long)nSpecies * sizeof(transition**)); 00298 00299 /* nSpecies should be global*/ 00300 /*The data array containing the species names should also be global*/ 00301 for( intNS = 0; intNS < nSpecies; intNS++ ) 00302 { 00303 if ( lgSpeciesMolecule[intNS] ) 00304 { 00305 atmdat_lamda_readin( intNS ); 00306 } 00307 else 00308 { 00309 atmdat_Chianti_readin( intNS ); 00310 } 00311 } 00312 return; 00313 } 00314 00315 STATIC int getAtNo(char *p) 00316 { 00317 DEBUG_ENTRY("getAtNo()"); 00318 if(strcmp(p,"H")==0) 00319 { 00320 return(ipHYDROGEN); 00321 } 00322 else if(strcmp(p,"HE")== 0) 00323 { 00324 return(ipHELIUM); 00325 } 00326 else if(strcmp(p,"C")== 0) 00327 { 00328 return(ipCARBON); 00329 } 00330 else if(strcmp(p,"N")== 0) 00331 { 00332 return(ipNITROGEN); 00333 } 00334 else if(strcmp(p,"O")== 0) 00335 { 00336 return(ipOXYGEN); 00337 } 00338 else if(strcmp(p,"NE")== 0) 00339 { 00340 return(ipNEON); 00341 } 00342 else if(strcmp(p,"NA")== 0) 00343 { 00344 return(ipSODIUM); 00345 } 00346 else if(strcmp(p,"MG")== 0) 00347 { 00348 return(ipMAGNESIUM); 00349 } 00350 else if(strcmp(p,"AL")== 0) 00351 { 00352 return(ipALUMINIUM); 00353 } 00354 else if(strcmp(p,"SI")== 0) 00355 { 00356 return(ipSILICON); 00357 } 00358 else if(strcmp(p,"P")== 0) 00359 { 00360 return(ipPHOSPHORUS); 00361 } 00362 else if(strcmp(p,"S")== 0) 00363 { 00364 return(ipSULPHUR); 00365 } 00366 else if(strcmp(p,"CL")== 0) 00367 { 00368 return(ipCHLORINE); 00369 } 00370 else if(strcmp(p,"AR")== 0) 00371 { 00372 return(ipARGON ); 00373 } 00374 else if(strcmp(p,"K")== 0) 00375 { 00376 return(ipPOTASSIUM); 00377 } 00378 else if(strcmp(p,"CA")== 0) 00379 { 00380 return(ipCALCIUM); 00381 } 00382 else if(strcmp(p,"SC")== 0) 00383 { 00384 return(ipSCANDIUM); 00385 } 00386 else if(strcmp(p,"TI")== 0) 00387 { 00388 return(ipTITANIUM); 00389 } 00390 else if(strcmp(p,"V")== 0) 00391 { 00392 return(ipVANADIUM); 00393 } 00394 else if(strcmp(p,"CR")== 0) 00395 { 00396 return(ipCHROMIUM ); 00397 } 00398 else if(strcmp(p,"MN")== 0) 00399 { 00400 return(ipMANGANESE); 00401 } 00402 else if(strcmp(p,"FE")== 0) 00403 { 00404 return(ipIRON); 00405 } 00406 else if(strcmp(p,"CO")== 0) 00407 { 00408 return(ipCOBALT); 00409 } 00410 else if(strcmp(p,"NI")== 0) 00411 { 00412 return(ipNICKEL); 00413 } 00414 else if(strcmp(p,"CU")== 0) 00415 { 00416 return(ipCOPPER); 00417 } 00418 else if(strcmp(p,"ZN")== 0) 00419 { 00420 return(ipZINC); 00421 } 00422 else 00423 { 00424 TotalInsanity(); 00425 } 00426 }