00001
00002
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 #include "mole.h"
00011
00012 #define DEBUGSTATE false
00013
00014
00015 void atmdat_LAMDA_readin( long intNS, char *chEFilename )
00016 {
00017 DEBUG_ENTRY( "atmdat_LAMDA_readin()" );
00018
00019 int nMolLevs = -10000,intCollIndex = -10000,intLColliderIndex = -10000;
00020
00021 long int intlnct,intrtct,intgrtct,intCollPart,
00022 intDCollPart,intCollTran,nCollTrans,intCollTemp,intcollindex;
00023
00024 FILE *ioLevData;
00025 realnum fstatwt,fenergyK,fenergyWN,fenergy,feinsteina;
00026
00027 char chLine[FILENAME_PATH_LENGTH_2] ,
00028
00029 *chColltemp,*chCollName;
00030
00031 ASSERT( intNS >= 0 );
00032
00033 const int MAX_NUM_LEVELS = 70;
00034
00035 Species[intNS].lgMolecular = true;
00036
00037
00038 if(DEBUGSTATE)
00039 printf("The name of the %li species is %s \n",intNS+1,Species[intNS].chLabel);
00040
00041
00042 if( trace.lgTrace )
00043 fprintf( ioQQQ," moldat_readin opening %s:",chEFilename);
00044
00045 ioLevData = open_data( chEFilename, "r" );
00046
00047 nMolLevs = 0;
00048 intrtct = 0;
00049 intgrtct = 0;
00050 intlnct = 0;
00051 while(intlnct < 3)
00052 {
00053 intlnct++;
00054 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00055 {
00056 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00057 cdEXIT(EXIT_FAILURE);
00058 }
00059 }
00060
00061 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00062 {
00063 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00064 cdEXIT(EXIT_FAILURE);
00065 }
00066 Species[intNS].fmolweight = (realnum)atof(chLine);
00067
00068
00069 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00070 {
00071 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00072 cdEXIT(EXIT_FAILURE);
00073 }
00074
00075
00076 ASSERT( chLine[0] == '!' );
00077
00078
00079 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00080 {
00081 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00082 cdEXIT(EXIT_FAILURE);
00083 }
00084 nMolLevs = atoi(chLine);
00085
00086 long HighestIndexInFile = nMolLevs;
00087
00088 nMolLevs = MIN2( nMolLevs, MAX_NUM_LEVELS );
00089
00090 if(nMolLevs <= 0)
00091 {
00092 fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEFilename );
00093 fprintf( ioQQQ, "The file must be corrupted.\n" );
00094 cdEXIT( EXIT_FAILURE );
00095 }
00096
00097 Species[intNS].numLevels_max = nMolLevs;
00098 Species[intNS].numLevels_local = Species[intNS].numLevels_max;
00099
00100
00101 dBaseStates[intNS] = (quantumState *)MALLOC((size_t)(nMolLevs)*sizeof(quantumState));
00102
00103 dBaseTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)nMolLevs);
00104 dBaseTrans[intNS][0] = NULL;
00105 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
00106 {
00107 dBaseTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi);
00108 for( long ipLo = 0; ipLo < ipHi; ipLo++)
00109 {
00110 dBaseTrans[intNS][ipHi][ipLo].Junk();
00111 dBaseTrans[intNS][ipHi][ipLo].Lo = &dBaseStates[intNS][ipLo];
00112 dBaseTrans[intNS][ipHi][ipLo].Hi = &dBaseStates[intNS][ipHi];
00113 }
00114 }
00115
00116 for( intcollindex = 0; intcollindex <NUM_COLLIDERS; intcollindex++ )
00117 {
00118 CollRatesArray[intNS][intcollindex] = NULL;
00119 }
00120
00121
00122 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00123 {
00124 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00125 cdEXIT(EXIT_FAILURE);
00126 }
00127
00128
00129 ASSERT( chLine[0] == '!' );
00130
00131 for( long ipLev=0; ipLev<HighestIndexInFile; ipLev++)
00132 {
00133 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00134 {
00135 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00136 cdEXIT(EXIT_FAILURE);
00137 }
00138
00139
00140 if( ipLev >= nMolLevs )
00141 continue;
00142
00143
00144 strcpy(dBaseStates[intNS][ipLev].chLabel, " ");
00145 strncpy(dBaseStates[intNS][ipLev].chLabel,Species[intNS].chLabel, 4);
00146 dBaseStates[intNS][ipLev].chLabel[4] = '\0';
00147
00148 if( dBaseStates[intNS][ipLev].chLabel[2]=='\0' )
00149 {
00150 dBaseStates[intNS][ipLev].chLabel[2]=' ';
00151 dBaseStates[intNS][ipLev].chLabel[3]=' ';
00152 }
00153 else if( dBaseStates[intNS][ipLev].chLabel[3]=='\0' )
00154 {
00155 dBaseStates[intNS][ipLev].chLabel[3]=' ';
00156 }
00157
00158
00159 dBaseStates[intNS][ipLev].sp = &Species[intNS];
00160
00161 long i = 1;
00162 bool lgEOL;
00163 long index;
00164
00165 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00166 fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00167 fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00168
00169 ASSERT( index == ipLev + 1 );
00170 dBaseStates[intNS][ipLev].energy.set(fenergy,"cm^-1");
00171 dBaseStates[intNS][ipLev].g = fstatwt;
00172
00173 if (ipLev > 0)
00174 {
00175 if (dBaseStates[intNS][ipLev].energy.WN() < dBaseStates[intNS][ipLev-1].energy.WN())
00176 {
00177 fprintf( ioQQQ, " The energy levels are not in order in species %s at index %li.\n",
00178 Species[intNS].chLabel, ipLev );
00179 cdEXIT(EXIT_FAILURE);
00180 }
00181 }
00182 if(DEBUGSTATE)
00183 {
00184 printf("The converted energy is %f \n",dBaseStates[intNS][ipLev].energy.WN());
00185 printf("The value of g is %f \n",dBaseStates[intNS][ipLev].g);
00186 }
00187 }
00188
00189
00190 for( long ipHi=1; ipHi<nMolLevs; ipHi++ )
00191 {
00192 for( long ipLo=0; ipLo<ipHi; ipLo++ )
00193 {
00194 fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy.WN() - dBaseStates[intNS][ipLo].energy.WN());
00195 fenergyK = (realnum)(fenergyWN*T1CM);
00196
00197 dBaseTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
00198 dBaseTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
00199 dBaseTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
00200
00201
00202
00203
00204 if( fenergyWN>SMALLFLOAT )
00205 dBaseTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8f/fenergyWN/
00206 RefIndex(fenergyWN));
00207 else
00208 dBaseTrans[intNS][ipHi][ipLo].WLAng = 1e30f;
00209 }
00210 }
00211
00212 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00213 {
00214 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00215 cdEXIT(EXIT_FAILURE);
00216 }
00217 if(chLine[0]!='!')
00218 {
00219 fprintf( ioQQQ, " The number of energy levels in file %s is not correct.\n",chEFilename);
00220 cdEXIT(EXIT_FAILURE);
00221 }
00222 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00223 {
00224 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00225 cdEXIT(EXIT_FAILURE);
00226 }
00227 intgrtct = atoi(chLine);
00228
00229 if(intgrtct <= 0)
00230 {
00231 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00232 cdEXIT(EXIT_FAILURE);
00233 }
00234 if(DEBUGSTATE)
00235 {
00236 printf("The number of radiative transitions is %li \n",intgrtct);
00237 }
00238 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00239 {
00240 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00241 cdEXIT(EXIT_FAILURE);
00242 }
00243
00244 for( intrtct=0; intrtct<intgrtct; intrtct++)
00245 {
00246 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00247 {
00248 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00249 cdEXIT(EXIT_FAILURE);
00250 }
00251
00252 long i = 1;
00253 bool lgEOL;
00254 long index, ipHi, ipLo;
00255
00256 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00257 ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
00258 ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
00259
00260 ASSERT( ipHi >= 0 );
00261 ASSERT( ipLo >= 0 );
00262
00263
00264 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
00265 continue;
00266
00267 feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00268
00269 FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00270 fenergyK = (realnum)((dBaseStates[intNS][ipHi].energy.WN() -dBaseStates[intNS][ipLo].energy.WN())*T1CM);
00271 ASSERT( index == intrtct + 1 );
00272
00273 dBaseTrans[intNS][ipHi][ipLo].Emis = AddLine2Stack(feinsteina, &dBaseTrans[intNS][ipHi][ipLo]);
00274 dBaseTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
00275 ASSERT( !isnan( dBaseTrans[intNS][ipHi][ipLo].EnergyK ) );
00276 fenergyWN = (realnum)((fenergyK)/T1CM);
00277 dBaseTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
00278 dBaseTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
00279 dBaseTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
00280
00281 dBaseTrans[intNS][ipHi][ipLo].Emis->gf = (realnum)GetGF(dBaseTrans[intNS][ipHi][ipLo].Emis->Aul,
00282 dBaseTrans[intNS][ipHi][ipLo].EnergyWN, dBaseTrans[intNS][ipHi][ipLo].Hi->g);
00283
00284 if(DEBUGSTATE)
00285 {
00286 printf("The upper level is %ld \n",ipHi+1);
00287 printf("The lower level is %ld \n",ipLo+1);
00288 printf("The Einstein A is %E \n",dBaseTrans[intNS][ipHi][ipLo].Emis->Aul);
00289 printf("The Energy of the transition is %E \n",dBaseTrans[intNS][ipHi][ipLo].EnergyK);
00290 }
00291 }
00292
00293 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00294 {
00295 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00296 cdEXIT(EXIT_FAILURE);
00297 }
00298
00299 if(chLine[0]!='!')
00300 {
00301 fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename);
00302 cdEXIT(EXIT_FAILURE);
00303 }
00304
00305
00306 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00307 {
00308 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00309 cdEXIT(EXIT_FAILURE);
00310 }
00311 else
00312 {
00313 intCollPart = atoi(chLine);
00314 }
00315
00316 if(intCollPart > NUM_COLLIDERS-1)
00317 {
00318 fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
00319 cdEXIT(EXIT_FAILURE);
00320 }
00321
00322 intDCollPart = intCollPart;
00323 while(intDCollPart > 0)
00324 {
00325 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00326 {
00327 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00328 cdEXIT(EXIT_FAILURE);
00329 }
00330
00331 ASSERT( chLine[0] == '!' );
00332
00333 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00334 {
00335 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00336 cdEXIT(EXIT_FAILURE);
00337 }
00338
00339
00340
00341
00342
00343
00344 chCollName = strtok(chLine," ");
00345
00346 intLColliderIndex = atoi(chCollName);
00347
00348
00349
00350 if(intLColliderIndex == 1)
00351 {
00352 intCollIndex = 8;
00353 }
00354 else if(intLColliderIndex == 2)
00355 {
00356 intCollIndex = 7;
00357 }
00358 else if(intLColliderIndex == 3)
00359 {
00360 intCollIndex = 6;
00361 }
00362 else if(intLColliderIndex == 4)
00363 {
00364 intCollIndex = 0;
00365 }
00366 else if(intLColliderIndex == 5)
00367 {
00368 intCollIndex = 2;
00369 }
00370 else if(intLColliderIndex == 6)
00371 {
00372 intCollIndex = 3;
00373 }
00374 else
00375 {
00376 TotalInsanity();
00377 }
00378
00379
00380 CollRatesArray[intNS][intCollIndex] = (double**)MALLOC((unsigned long)(nMolLevs) * sizeof(double*));
00381 for( long ipHi = 0; ipHi<nMolLevs; ipHi++ )
00382 {
00383 CollRatesArray[intNS][intCollIndex][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs) * sizeof(double));
00384 for( long ipLo = 0; ipLo<nMolLevs; ipLo++ )
00385 {
00386 CollRatesArray[intNS][intCollIndex][ipHi][ipLo] = 0.0;
00387 }
00388 }
00389 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00390 {
00391 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00392 cdEXIT(EXIT_FAILURE);
00393 }
00394
00395 ASSERT( chLine[0] == '!' );
00396
00397 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00398 {
00399 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00400 cdEXIT(EXIT_FAILURE);
00401 }
00402
00403 intCollTran = atoi(chLine);
00404
00405 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00406 {
00407 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00408 cdEXIT(EXIT_FAILURE);
00409 }
00410
00411 ASSERT( chLine[0] == '!' );
00412
00413 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00414 {
00415 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00416 cdEXIT(EXIT_FAILURE);
00417 }
00418
00419 intCollTemp = atoi(chLine);
00420
00421 AtmolCollRateCoeff[intNS][intCollIndex].ntemps = intCollTemp;
00422
00423 AtmolCollRateCoeff[intNS][intCollIndex].temps =
00424 (double *)MALLOC((unsigned long)intCollTemp*sizeof(double));
00425 AtmolCollRateCoeff[intNS][intCollIndex].collrates =
00426 (double***)MALLOC((unsigned long)(nMolLevs)*sizeof(double**));
00427 for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
00428 {
00429 AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi] =
00430 (double **)MALLOC((unsigned long)(nMolLevs)*sizeof(double*));
00431 for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
00432 {
00433 AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo] =
00434 (double *)MALLOC((unsigned long)(intCollTemp)*sizeof(double));
00435
00436 memset( AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo], 0, (unsigned long)(intCollTemp)*sizeof(double) );
00437 }
00438 }
00439
00440
00441 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00442 {
00443 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00444 cdEXIT(EXIT_FAILURE);
00445 }
00446
00447 ASSERT( chLine[0] == '!' );
00448
00449
00450 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00451 {
00452 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00453 cdEXIT(EXIT_FAILURE);
00454 }
00455
00456 chColltemp =strtok(chLine," ");
00457 AtmolCollRateCoeff[intNS][intCollIndex].temps[0] =(realnum) atof(chColltemp);
00458 for( long ipTe=1; ipTe<intCollTemp; ipTe++ )
00459 {
00460 chColltemp =strtok(NULL," ");
00461 AtmolCollRateCoeff[intNS][intCollIndex].temps[ipTe] = atof(chColltemp);
00462 }
00463
00464 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00465 {
00466 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00467 cdEXIT(EXIT_FAILURE);
00468 }
00469
00470 for( nCollTrans=0; nCollTrans<intCollTran; nCollTrans++ )
00471 {
00472 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00473 {
00474 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00475 cdEXIT(EXIT_FAILURE);
00476 }
00477
00478 long i = 1;
00479 bool lgEOL;
00480 long index, ipHi, ipLo;
00481
00482 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00483 ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
00484 ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
00485
00486
00487 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
00488 continue;
00489
00490
00491 if( ipHi < ipLo )
00492 {
00493 ASSERT( ipLo == nMolLevs - 1);
00494 long temp = ipHi;
00495 ipHi = ipLo;
00496 ipLo = temp;
00497 }
00498
00499 for( long j=0; j<intCollTemp; j++ )
00500 {
00501 AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][j] =
00502 (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00503 }
00504
00505 ASSERT( index == nCollTrans + 1 );
00506
00507 if(DEBUGSTATE)
00508 {
00509 printf("The values of up and lo are %ld & %ld \n",ipHi,ipLo);
00510 printf("The collision rates are");
00511 for(i=0;i<intCollTemp;i++)
00512 {
00513 printf("\n %e",AtmolCollRateCoeff[intNS][intCollIndex].collrates[ipHi][ipLo][i]);
00514 }
00515 printf("\n");
00516 }
00517 }
00518
00519 intDCollPart = intDCollPart -1;
00520 }
00521
00522 fclose( ioLevData );
00523
00524 return;
00525 }