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 "physconst.h" 00006 #include "taulines.h" 00007 #include "trace.h" 00008 #include "string.h" 00009 #include "thirdparty.h" 00010 00011 #define DEBUGSTATE false 00012 00013 /*Separate Routine to read in the molecules*/ 00014 void atmdat_lamda_readin( long intNS ) 00015 { 00016 DEBUG_ENTRY( "atmdat_lamda_readin()" ); 00017 00018 int ngMolLevs = -10000,intCollIndex = -10000,intLColliderIndex = -10000; 00019 /* type is set to 0 for non chianti and 1 for chianti*/ 00020 long int i,j,nMolLevs,intlnct,intrtct,intgrtct,intCollPart, 00021 intDCollPart,intCollTran,nCollTrans,intCollTemp,intcollindex, 00022 ipLo,ipHi; 00023 /*intrtct refers to radiative transitions count*/ 00024 FILE *atmolLevDATA; 00025 realnum fstatwt,fenergyK,fenergyWN,fenergy,feinsteina; 00026 00027 char chLine[FILENAME_PATH_LENGTH_2] , 00028 chEFilename[FILENAME_PATH_LENGTH_2], 00029 *chColltemp,*chCollName; 00030 00031 ASSERT( lgSpeciesMolecule[intNS] ); 00032 00033 /*Load the species name in the Species array structure*/ 00034 if(DEBUGSTATE) 00035 printf("The name of the %li species is %s \n",intNS+1,Species[intNS].chptrSpName); 00036 strcpy( chEFilename , Species[intNS].chptrSpName ); 00037 strcat( chEFilename , ".dat"); 00038 uncaps( chEFilename ); 00039 00040 /*Open the files*/ 00041 if( trace.lgTrace ) 00042 fprintf( ioQQQ," moldat_readin opening %s:",chEFilename); 00043 00044 atmolLevDATA = open_data( chEFilename, "r" ); 00045 00046 nMolLevs = 0; 00047 ngMolLevs = 0; 00048 intrtct = 0; 00049 intgrtct = 0; 00050 intlnct = 0; 00051 ipHi = 0; 00052 ipLo = 0; 00053 j = 0; 00054 while(intlnct < 3) 00055 { 00056 intlnct++; 00057 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00058 { 00059 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00060 cdEXIT(EXIT_FAILURE); 00061 } 00062 } 00063 /*Extracting out the molecular weight*/ 00064 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00065 { 00066 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00067 cdEXIT(EXIT_FAILURE); 00068 } 00069 Species[intNS].fmolweight = (realnum)atof(chLine); 00070 00071 /*Discard this line*/ 00072 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00073 { 00074 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00075 cdEXIT(EXIT_FAILURE); 00076 } 00077 /*Reading in the number of energy levels*/ 00078 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00079 { 00080 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00081 cdEXIT(EXIT_FAILURE); 00082 } 00083 ngMolLevs = atoi(chLine); 00084 00085 if(ngMolLevs <= 0) 00086 { 00087 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00088 cdEXIT(EXIT_FAILURE); 00089 } 00090 Species[intNS].numLevels_max = ngMolLevs; 00091 Species[intNS].numLevels_local = ngMolLevs; 00092 if(DEBUGSTATE) 00093 { 00094 printf("The number of energy levels is %li \n",Species[intNS].numLevels_max); 00095 } 00096 /*the above refers to the number of energy levels specified in the data file*/ 00097 /*malloc the States array*/ 00098 atmolStates[intNS] = (quantumState *)MALLOC((size_t)(ngMolLevs)*sizeof(quantumState)); 00099 /*malloc the Transition array*/ 00100 atmolTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)ngMolLevs); 00101 atmolTrans[intNS][0] = NULL; 00102 for( ipHi = 1; ipHi < ngMolLevs; ipHi++) 00103 { 00104 atmolTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi); 00105 for(ipLo = 0; ipLo < ipHi; ipLo++) 00106 { 00107 atmolTrans[intNS][ipHi][ipLo].Lo = &atmolStates[intNS][ipLo]; 00108 atmolTrans[intNS][ipHi][ipLo].Hi = &atmolStates[intNS][ipHi]; 00109 atmolTrans[intNS][ipHi][ipLo].Emis = &DummyEmis; 00110 } 00111 } 00112 00113 for( intcollindex = 0; intcollindex <NUM_COLLIDERS; intcollindex++ ) 00114 { 00115 CollRatesArray[intNS][intcollindex] = NULL; 00116 } 00117 00118 00119 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00120 { 00121 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00122 cdEXIT(EXIT_FAILURE); 00123 } 00124 00125 for( nMolLevs=0; nMolLevs<ngMolLevs; nMolLevs++) 00126 { 00127 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00128 { 00129 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00130 cdEXIT(EXIT_FAILURE); 00131 } 00132 /*information needed for label*/ 00133 strcpy(atmolStates[intNS][nMolLevs].chLabel, " "); 00134 strncpy(atmolStates[intNS][nMolLevs].chLabel,Species[intNS].chptrSpName, 4); 00135 atmolStates[intNS][nMolLevs].chLabel[4] = '\0'; 00136 00137 long i = 1; 00138 bool lgEOL; 00139 long index; 00140 00141 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00142 fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00143 fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00144 00145 ASSERT( index == nMolLevs + 1 ); 00146 00147 atmolStates[intNS][nMolLevs].energy = fenergy; 00148 atmolStates[intNS][nMolLevs].g = fstatwt; 00149 00150 if (nMolLevs > 0) 00151 { 00152 if (atmolStates[intNS][nMolLevs].energy < atmolStates[intNS][nMolLevs-1].energy) 00153 { 00154 fprintf( ioQQQ, " The energy levels are not in order"); 00155 cdEXIT(EXIT_FAILURE); 00156 } 00157 } 00158 if(DEBUGSTATE) 00159 { 00160 printf("The converted energy is %f \n",atmolStates[intNS][nMolLevs].energy); 00161 printf("The value of g is %f \n",atmolStates[intNS][nMolLevs].g); 00162 } 00163 } 00164 00165 /* fill in all transition energies, can later overwrite for specific radiative transitions */ 00166 for( i=1; i<ngMolLevs; i++ ) 00167 { 00168 for( j=0; j<i; j++ ) 00169 { 00170 fenergyWN = atmolStates[intNS][i].energy - atmolStates[intNS][j].energy; 00171 fenergyK = (realnum)(fenergyWN*T1CM); 00172 00173 atmolTrans[intNS][i][j].EnergyK = fenergyK; 00174 atmolTrans[intNS][i][j].EnergyWN = fenergyWN; 00175 atmolTrans[intNS][i][j].EnergyErg = (realnum)ERG1CM *fenergyWN; 00176 00177 /* protect against division by zero. */ 00178 atmolTrans[intNS][i][j].WLAng = (realnum)(1e+8f/MAX2(1E-5f,fenergyWN)/ 00179 RefIndex(fenergyWN)); 00180 } 00181 } 00182 00183 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00184 { 00185 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00186 cdEXIT(EXIT_FAILURE); 00187 } 00188 if(chLine[0]!='!') 00189 { 00190 fprintf( ioQQQ, " The number of energy levels in file %s is not correct.\n",chEFilename); 00191 cdEXIT(EXIT_FAILURE); 00192 } 00193 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00194 { 00195 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00196 cdEXIT(EXIT_FAILURE); 00197 } 00198 intgrtct = atoi(chLine); 00199 /*The above gives the number of radiative transitions*/ 00200 if(intgrtct <= 0) 00201 { 00202 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00203 cdEXIT(EXIT_FAILURE); 00204 } 00205 if(DEBUGSTATE) 00206 { 00207 printf("The number of radiative transitions is %li \n",intgrtct); 00208 } 00209 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00210 { 00211 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00212 cdEXIT(EXIT_FAILURE); 00213 } 00214 00215 for( intrtct=0; intrtct<intgrtct; intrtct++) 00216 { 00217 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00218 { 00219 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00220 cdEXIT(EXIT_FAILURE); 00221 } 00222 00223 long i = 1; 00224 bool lgEOL; 00225 long index; 00226 00227 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00228 ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1; 00229 ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1; 00230 feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00231 /* don't need the energy in GHz, so throw it away. */ 00232 FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00233 fenergyK = (realnum)(((atmolStates[intNS][ipHi].energy) -(atmolStates[intNS][ipLo].energy))*T1CM); 00234 ASSERT( index == intrtct + 1 ); 00235 00236 atmolTrans[intNS][ipHi][ipLo].Emis = AddLine2Stack(feinsteina, &atmolTrans[intNS][ipHi][ipLo]); 00237 atmolTrans[intNS][ipHi][ipLo].EnergyK = fenergyK; 00238 fenergyWN = (realnum)((fenergyK)/T1CM); 00239 atmolTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN; 00240 atmolTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN; 00241 atmolTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN)); 00242 00243 ASSERT( !isnan( atmolTrans[intNS][ipHi][ipLo].EnergyK ) ); 00244 00245 if(DEBUGSTATE) 00246 { 00247 printf("The upper level is %ld \n",ipHi+1); 00248 printf("The lower level is %ld \n",ipLo+1); 00249 printf("The Einstein A is %E \n",atmolTrans[intNS][ipHi][ipLo].Emis->Aul); 00250 printf("The Energy of the transition is %E \n",atmolTrans[intNS][ipHi][ipLo].EnergyK); 00251 } 00252 } 00253 00254 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00255 { 00256 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00257 cdEXIT(EXIT_FAILURE); 00258 } 00259 00260 if(chLine[0]!='!') 00261 { 00262 fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename); 00263 cdEXIT(EXIT_FAILURE); 00264 } 00265 00266 /*Getting the number of collisional partners*/ 00267 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00268 { 00269 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00270 cdEXIT(EXIT_FAILURE); 00271 } 00272 else 00273 { 00274 intCollPart = atoi(chLine); 00275 } 00276 /*Checking the number of collisional partners does not exceed 9*/ 00277 if(intCollPart > NUM_COLLIDERS-1) 00278 { 00279 fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename ); 00280 cdEXIT(EXIT_FAILURE); 00281 } 00282 /*Creating the duplicate of the number of collisional partners which is reduced*/ 00283 intDCollPart = intCollPart; 00284 while(intDCollPart > 0) 00285 { 00286 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00287 { 00288 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00289 cdEXIT(EXIT_FAILURE); 00290 } 00291 00292 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00293 { 00294 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00295 cdEXIT(EXIT_FAILURE); 00296 } 00297 /*Extract out the name of the collider*/ 00298 /*The following are the rules expected in the datafiles to extract the names of the collider*/ 00299 /*The line which displays the species and the collider starts with a number*/ 00300 /*This refers to the collider in the Leiden database*/ 00301 /*In the Leiden database 1 referes to H2,2 to para-H2,3 to ortho-H2 00302 4 to electrons,5 to H and 6 to He*/ 00303 chCollName = strtok(chLine," "); 00304 /*Leiden Collider Index*/ 00305 intLColliderIndex = atoi(chCollName); 00306 /*In Cloudy,We assign the following indices for the colliders*/ 00307 /*electron=0,proton=1,atomic hydrogen=2,He=3,He+=4,He++=5,oH2=6,pH2=7,H2=8*/ 00308 00309 if(intLColliderIndex == 1) 00310 { 00311 intCollIndex = 8; 00312 } 00313 else if(intLColliderIndex == 2) 00314 { 00315 intCollIndex = 7; 00316 } 00317 else if(intLColliderIndex == 3) 00318 { 00319 intCollIndex = 6; 00320 } 00321 else if(intLColliderIndex == 4) 00322 { 00323 intCollIndex = 0; 00324 } 00325 else if(intLColliderIndex == 5) 00326 { 00327 intCollIndex = 2; 00328 } 00329 else if(intLColliderIndex == 6) 00330 { 00331 intCollIndex = 3; 00332 } 00333 else 00334 { 00335 TotalInsanity(); 00336 } 00337 /*This is where we allocate memory if the collider exists*/ 00338 /*Needed to take care of the he collisions*/ 00339 CollRatesArray[intNS][intCollIndex] = (double**)MALLOC((unsigned long)(ngMolLevs) * sizeof(double*)); 00340 for( ipHi = 0; ipHi<ngMolLevs; ipHi++ ) 00341 { 00342 CollRatesArray[intNS][intCollIndex][ipHi] = (double*)MALLOC((unsigned long)(ngMolLevs) * sizeof(double)); 00343 for( ipLo = 0; ipLo<ngMolLevs; ipLo++ ) 00344 { 00345 CollRatesArray[intNS][intCollIndex][ipHi][ipLo] = 0.0; 00346 } 00347 } 00348 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00349 { 00350 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00351 cdEXIT(EXIT_FAILURE); 00352 } 00353 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00354 { 00355 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00356 cdEXIT(EXIT_FAILURE); 00357 } 00358 /*Number of coll trans*/ 00359 intCollTran = atoi(chLine); 00360 00361 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00362 { 00363 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00364 cdEXIT(EXIT_FAILURE); 00365 } 00366 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00367 { 00368 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00369 cdEXIT(EXIT_FAILURE); 00370 } 00371 /*Number of coll temps*/ 00372 intCollTemp = atoi(chLine); 00373 /*Storing the number of collisional temperatures*/ 00374 AtmolCollRateCoeff[intNS][intCollIndex].ntemps = intCollTemp; 00375 /*Mallocing*/ 00376 AtmolCollRateCoeff[intNS][intCollIndex].temps = 00377 (double *)MALLOC((unsigned long)intCollTemp*sizeof(double)); 00378 AtmolCollRateCoeff[intNS][intCollIndex].collrates = 00379 (double***)MALLOC((unsigned long)(ngMolLevs)*sizeof(double**)); 00380 for( i=0; i<ngMolLevs; i++ ) 00381 { 00382 AtmolCollRateCoeff[intNS][intCollIndex].collrates[i] = 00383 (double **)MALLOC((unsigned long)(ngMolLevs)*sizeof(double*)); 00384 for( j=0; j<ngMolLevs; j++ ) 00385 { 00386 AtmolCollRateCoeff[intNS][intCollIndex].collrates[i][j] = 00387 (double *)MALLOC((unsigned long)(intCollTemp)*sizeof(double)); 00388 } 00389 } 00390 00391 /*Discard the header line*/ 00392 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00393 { 00394 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00395 cdEXIT(EXIT_FAILURE); 00396 } 00397 /*Getting the collisional Temperatures*/ 00398 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00399 { 00400 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00401 cdEXIT(EXIT_FAILURE); 00402 } 00403 /*Filling the collisional temps array*/ 00404 chColltemp =strtok(chLine," "); 00405 AtmolCollRateCoeff[intNS][intCollIndex].temps[0] =(realnum) atof(chColltemp); 00406 for( i=1; i<intCollTemp; i++ ) 00407 { 00408 chColltemp =strtok(NULL," "); 00409 AtmolCollRateCoeff[intNS][intCollIndex].temps[i] = atof(chColltemp); 00410 } 00411 /*Discard the header line*/ 00412 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00413 { 00414 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00415 cdEXIT(EXIT_FAILURE); 00416 } 00417 /*Getting all the collisional transitions data*/ 00418 for( nCollTrans=0; nCollTrans<intCollTran; nCollTrans++ ) 00419 { 00420 if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL ) 00421 { 00422 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename); 00423 cdEXIT(EXIT_FAILURE); 00424 } 00425 00426 long i = 1; 00427 bool lgEOL; 00428 long index; 00429 00430 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00431 ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1; 00432 ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1; 00433 00434 /* Indices between the very highest levels seem to be reversed */ 00435 if( ipHi < ipLo ) 00436 { 00437 ASSERT( ipLo == ngMolLevs - 1); 00438 long temp = ipHi; 00439 ipHi = ipLo; 00440 ipLo = temp; 00441 } 00442 00443 for( long j=0; j<intCollTemp; j++ ) 00444 { 00445 AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][j] = 00446 (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ); 00447 } 00448 00449 ASSERT( index == nCollTrans + 1 ); 00450 00451 if(DEBUGSTATE) 00452 { 00453 printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo); 00454 printf("The collision rates are"); 00455 for(i=0;i<intCollTemp;i++) 00456 { 00457 printf("\n %e",AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][i]); 00458 } 00459 printf("\n"); 00460 } 00461 } 00462 00463 intDCollPart = intDCollPart -1; 00464 } 00465 00466 return; 00467 }