00001
00002
00003 #include "cddefines.h"
00004 #include "atmdat.h"
00005 #include "lines_service.h"
00006 #include "physconst.h"
00007 #include "rfield.h"
00008 #include "taulines.h"
00009 #include "trace.h"
00010 #include "string.h"
00011 #include "thirdparty.h"
00012 #include "mole.h"
00013 #include "rfield.h"
00014
00015 #define DEBUGSTATE false
00016
00017
00018 void atmdat_LAMDA_readin( long intNS, char *chEFilename )
00019 {
00020 DEBUG_ENTRY( "atmdat_LAMDA_readin()" );
00021
00022 int nMolLevs = -10000;
00023 long int intlnct,intrtct,intgrtct;
00024
00025
00026 FILE *ioLevData;
00027 realnum fstatwt,fenergyK,fenergyWN,fenergy,feinsteina;
00028
00029 char chLine[FILENAME_PATH_LENGTH_2];
00030
00031 ASSERT( intNS >= 0 );
00032
00033 const int MAX_NUM_LEVELS = 70;
00034
00035 dBaseSpecies[intNS].lgMolecular = true;
00036 dBaseSpecies[intNS].lgLAMDA = true;
00037
00038
00039 if(DEBUGSTATE)
00040 printf("The name of the %li species is %s \n",intNS+1,dBaseSpecies[intNS].chLabel);
00041
00042
00043 if( trace.lgTrace )
00044 fprintf( ioQQQ," moldat_readin opening %s:",chEFilename);
00045
00046 ioLevData = open_data( chEFilename, "r" );
00047
00048 nMolLevs = 0;
00049 intrtct = 0;
00050 intgrtct = 0;
00051 intlnct = 0;
00052 while(intlnct < 3)
00053 {
00054 intlnct++;
00055 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00056 {
00057 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00058 cdEXIT(EXIT_FAILURE);
00059 }
00060 }
00061
00062 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00063 {
00064 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00065 cdEXIT(EXIT_FAILURE);
00066 }
00067 dBaseSpecies[intNS].fmolweight = (realnum)atof(chLine);
00068
00069
00070 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00071 {
00072 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00073 cdEXIT(EXIT_FAILURE);
00074 }
00075
00076
00077 ASSERT( chLine[0] == '!' );
00078
00079
00080 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00081 {
00082 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00083 cdEXIT(EXIT_FAILURE);
00084 }
00085 nMolLevs = atoi(chLine);
00086
00087 long HighestIndexInFile = nMolLevs;
00088
00089 nMolLevs = MIN2( nMolLevs, MAX_NUM_LEVELS );
00090
00091 if(nMolLevs <= 0)
00092 {
00093 fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEFilename );
00094 fprintf( ioQQQ, "The file must be corrupted.\n" );
00095 cdEXIT( EXIT_FAILURE );
00096 }
00097
00098 dBaseSpecies[intNS].numLevels_max = nMolLevs;
00099 dBaseSpecies[intNS].numLevels_local = dBaseSpecies[intNS].numLevels_max;
00100
00101
00102 dBaseStates[intNS].resize(nMolLevs);
00103
00104 ipdBaseTrans[intNS].reserve(nMolLevs);
00105 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
00106 ipdBaseTrans[intNS].reserve(ipHi,ipHi);
00107 ipdBaseTrans[intNS].alloc();
00108 dBaseTrans[intNS].resize(ipdBaseTrans[intNS].size());
00109 dBaseTrans[intNS].states() = &dBaseStates[intNS];
00110 dBaseTrans[intNS].chLabel() = "LAMDA";
00111
00112 int ndBase = 0;
00113 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
00114 {
00115 for( long ipLo = 0; ipLo < ipHi; ipLo++)
00116 {
00117 ipdBaseTrans[intNS][ipHi][ipLo] = ndBase;
00118 dBaseTrans[intNS][ndBase].Junk();
00119 dBaseTrans[intNS][ndBase].setLo(ipLo);
00120 dBaseTrans[intNS][ndBase].setHi(ipHi);
00121 ++ndBase;
00122 }
00123 }
00124
00125 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00126 {
00127 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00128 cdEXIT(EXIT_FAILURE);
00129 }
00130
00131
00132 ASSERT( chLine[0] == '!' );
00133
00134 for( long ipLev=0; ipLev<HighestIndexInFile; ipLev++)
00135 {
00136 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00137 {
00138 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00139 cdEXIT(EXIT_FAILURE);
00140 }
00141
00142
00143 if( ipLev >= nMolLevs )
00144 continue;
00145
00146
00147 strcpy(dBaseStates[intNS][ipLev].chLabel(), " ");
00148 strncpy(dBaseStates[intNS][ipLev].chLabel(),dBaseSpecies[intNS].chLabel, 4);
00149
00150
00151 if( nMatch( "^13C", dBaseStates[intNS][ipLev].chLabel() ) )
00152 strcpy(dBaseStates[intNS][ipLev].chLabel(), "13CO");
00153
00154 dBaseStates[intNS][ipLev].chLabel()[4] = '\0';
00155
00156 if( dBaseStates[intNS][ipLev].chLabel()[2]=='\0' )
00157 {
00158 dBaseStates[intNS][ipLev].chLabel()[2]=' ';
00159 dBaseStates[intNS][ipLev].chLabel()[3]=' ';
00160 }
00161 else if( dBaseStates[intNS][ipLev].chLabel()[3]=='\0' )
00162 {
00163 dBaseStates[intNS][ipLev].chLabel()[3]=' ';
00164 }
00165
00166 long i = 1;
00167 bool lgEOL;
00168 long index;
00169
00170 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00171 fenergy = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00172 fstatwt = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00173
00174 ASSERT( index == ipLev + 1 );
00175 dBaseStates[intNS][ipLev].energy().set(fenergy,"cm^-1");
00176 dBaseStates[intNS][ipLev].g() = fstatwt;
00177
00178 if( ipLev > 0 )
00179 {
00180
00181
00182 volatile realnum enlev = dBaseStates[intNS][ipLev].energy().Ryd();
00183 volatile realnum enprev = dBaseStates[intNS][ipLev-1].energy().Ryd();
00184 if (enlev < enprev)
00185 {
00186 fprintf( ioQQQ, " The energy levels are not in order in species %s at index %li.\n",
00187 dBaseSpecies[intNS].chLabel, ipLev );
00188 cdEXIT(EXIT_FAILURE);
00189 }
00190 }
00191 if(DEBUGSTATE)
00192 {
00193 printf("The converted energy is %f \n",dBaseStates[intNS][ipLev].energy().WN());
00194 printf("The value of g is %f \n",dBaseStates[intNS][ipLev].g());
00195 }
00196 }
00197
00198
00199 for(TransitionList::iterator tr=dBaseTrans[intNS].begin();
00200 tr!= dBaseTrans[intNS].end(); ++tr)
00201 {
00202 int ipHi = (*tr).ipHi();
00203 int ipLo = (*tr).ipLo();
00204
00205 fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy().WN() - dBaseStates[intNS][ipLo].energy().WN());
00206 fenergyK = (realnum)(fenergyWN*T1CM);
00207
00208 (*tr).EnergyWN() = fenergyWN;
00209
00210
00211
00212
00213 if( fenergyWN>SMALLFLOAT )
00214 (*tr).WLAng() = (realnum)(1e+8f/fenergyWN/RefIndex(fenergyWN));
00215 else
00216 (*tr).WLAng() = 1e30f;
00217 }
00218
00219 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00220 {
00221 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00222 cdEXIT(EXIT_FAILURE);
00223 }
00224 if(chLine[0]!='!')
00225 {
00226 fprintf( ioQQQ, " The number of energy levels in file %s is not correct, expected to find line starting with!.\n",chEFilename);
00227 fprintf( ioQQQ , "%s\n", chLine );
00228 cdEXIT(EXIT_FAILURE);
00229 }
00230 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00231 {
00232 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00233 cdEXIT(EXIT_FAILURE);
00234 }
00235 intgrtct = atoi(chLine);
00236
00237 if(intgrtct <= 0)
00238 {
00239 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00240 cdEXIT(EXIT_FAILURE);
00241 }
00242 if(DEBUGSTATE)
00243 {
00244 printf("The number of radiative transitions is %li \n",intgrtct);
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 for( intrtct=0; intrtct<intgrtct; intrtct++)
00253 {
00254 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00255 {
00256 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00257 cdEXIT(EXIT_FAILURE);
00258 }
00259
00260 long i = 1;
00261 bool lgEOL;
00262 long index, ipHi, ipLo;
00263
00264 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00265 ipHi = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
00266 ipLo = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL ) - 1;
00267
00268 ASSERT( ipHi >= 0 );
00269 ASSERT( ipLo >= 0 );
00270 ASSERT( index == intrtct + 1 );
00271
00272
00273 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
00274 continue;
00275
00276 feinsteina = (realnum)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00277
00278 FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
00279 fenergyWN = (realnum)(dBaseStates[intNS][ipHi].energy().WN() -dBaseStates[intNS][ipLo].energy().WN());
00280 fenergyWN = MAX2( fenergyWN, 1.01 * RYD_INF * rfield.emm );
00281 fenergyK = (realnum)(fenergyWN*T1CM);
00282
00283 TransitionList::iterator tr = dBaseTrans[intNS].begin()+ipdBaseTrans[intNS][ipHi][ipLo];
00284 ASSERT(!(*tr).hasEmis());
00285 (*tr).AddLine2Stack();
00286 (*tr).Emis().Aul() = feinsteina;
00287 ASSERT( !isnan( (*tr).EnergyK() ) );
00288 fenergyWN = (realnum)((fenergyK)/T1CM);
00289 (*tr).EnergyWN() = fenergyWN;
00290 (*tr).WLAng() = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
00291
00292 (*tr).Emis().gf() = (realnum)GetGF((*tr).Emis().Aul(),(*tr).EnergyWN(), (*(*tr).Hi()).g());
00293
00294 if(DEBUGSTATE)
00295 {
00296 printf("The upper level is %ld \n",ipHi+1);
00297 printf("The lower level is %ld \n",ipLo+1);
00298 printf("The Einstein A is %E \n",(*tr).Emis().Aul());
00299 printf("The Energy of the transition is %E \n",(*tr).EnergyK());
00300 }
00301 }
00302
00303 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00304 {
00305 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00306 cdEXIT(EXIT_FAILURE);
00307 }
00308
00309 if(chLine[0]!='!')
00310 {
00311 fprintf( ioQQQ, " The number of radiative transitions in file %s is not correct.\n",chEFilename);
00312 cdEXIT(EXIT_FAILURE);
00313 }
00314
00315 long nCollPartners = -1;
00316
00317 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00318 {
00319 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00320 cdEXIT(EXIT_FAILURE);
00321 }
00322 else
00323 {
00324 nCollPartners = atoi(chLine);
00325 }
00326
00327 dBaseSpecies[intNS].lgLTE = ( nCollPartners < 0 );
00328
00329 if( nCollPartners > ipNCOLLIDER )
00330 {
00331 fprintf( ioQQQ, " The number of colliders is greater than what is expected in file %s.\n", chEFilename );
00332 cdEXIT(EXIT_FAILURE);
00333 }
00334
00335 FunctPtr func = new FunctLAMDA();
00336
00337 for( long ipPartner = 0; ipPartner < nCollPartners; ++ ipPartner )
00338 {
00339 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00340 {
00341 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00342 cdEXIT(EXIT_FAILURE);
00343 }
00344
00345 ASSERT( chLine[0] == '!' );
00346
00347 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00348 {
00349 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00350 cdEXIT(EXIT_FAILURE);
00351 }
00352
00353
00354
00355
00356
00357
00358 char *chCollName = strtok(chLine," ");
00359
00360 int intLColliderIndex = atoi(chCollName);
00361 int intCollIndex = -1;
00362
00363
00364
00365 if(intLColliderIndex == 1)
00366 {
00367 intCollIndex = ipH2;
00368 }
00369 else if(intLColliderIndex == 2)
00370 {
00371 intCollIndex = ipH2_PARA;
00372 }
00373 else if(intLColliderIndex == 3)
00374 {
00375 intCollIndex = ipH2_ORTHO;
00376 }
00377 else if(intLColliderIndex == 4)
00378 {
00379 intCollIndex = ipELECTRON;
00380 }
00381 else if(intLColliderIndex == 5)
00382 {
00383 intCollIndex = ipATOM_H;
00384 }
00385 else if(intLColliderIndex == 6)
00386 {
00387 intCollIndex = ipATOM_HE;
00388 }
00389 else
00390 {
00391
00392
00393 TotalInsanity();
00394 }
00395
00396 ASSERT( intCollIndex < ipNCOLLIDER );
00397
00398 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00399 {
00400 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00401 cdEXIT(EXIT_FAILURE);
00402 }
00403
00404 ASSERT( chLine[0] == '!' );
00405
00406 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00407 {
00408 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00409 cdEXIT(EXIT_FAILURE);
00410 }
00411
00412 long nCollTrans = atoi(chLine);
00413
00414 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00415 {
00416 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00417 cdEXIT(EXIT_FAILURE);
00418 }
00419
00420 ASSERT( chLine[0] == '!' );
00421
00422 if(read_whole_line( chLine , (int)sizeof(chLine) , ioLevData ) == NULL )
00423 {
00424 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEFilename);
00425 cdEXIT(EXIT_FAILURE);
00426 }
00427
00428 long intCollTemp = atoi(chLine);
00429
00430
00431 ReadCollisionRateTable( AtmolCollRateCoeff[intNS][intCollIndex],
00432 ioLevData, func, nMolLevs, intCollTemp, nCollTrans );
00433 }
00434 delete func;
00435
00436 fclose( ioLevData );
00437
00438 return;
00439 }
00440