00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "physconst.h"
00007 #include "taulines.h"
00008 #include "mewecoef.h"
00009 #include "iterations.h"
00010 #include "heavy.h"
00011 #include "mole.h"
00012 #include "yield.h"
00013 #include "trace.h"
00014 #include "lines.h"
00015 #include "lines_service.h"
00016 #include "ionbal.h"
00017 #include "struc.h"
00018 #include "geometry.h"
00019 #include "dense.h"
00020 #include "dynamics.h"
00021 #include "elementnames.h"
00022 #include "hyperfine.h"
00023 #include "atmdat.h"
00024 #include "iso.h"
00025
00026
00027
00028 const long INTBIG = 2000000000;
00029
00030
00031
00032
00033
00034
00035 long ipT1656=INTBIG, ipT9830=INTBIG,
00036 ipT8727=INTBIG, ipT1335=INTBIG, ipT1909=INTBIG,
00037 ipT977=INTBIG, ipT1550=INTBIG, ipT1548=INTBIG, ipT386=INTBIG,
00038 ipT310=INTBIG, ipc31175=INTBIG, ipT291=INTBIG, ipT280=INTBIG,
00039 ipT274=INTBIG, ipT270=INTBIG, ipT312=INTBIG, ipT610=INTBIG,
00040 ipT370=INTBIG, ipT157=INTBIG, ipT1085=INTBIG,
00041 ipT990=INTBIG, ipT1486=INTBIG, ipT765=INTBIG, ipT1243=INTBIG,
00042 ipT1239=INTBIG, ipT374g=INTBIG, ipT374x=INTBIG, ipT1200=INTBIG,
00043 ipT2140=INTBIG, ipT671=INTBIG, ipT315=INTBIG, ipT324=INTBIG,
00044 ipT333=INTBIG, ipT209=INTBIG, ipT122=INTBIG, ipT205=INTBIG,
00045 ipT57=INTBIG, ipT6300=INTBIG, ipT6363=INTBIG, ipT5577=INTBIG,
00046 ipT834=INTBIG, ipT1661=INTBIG, ipT1666=INTBIG, ipT835=INTBIG,
00047 ipT789=INTBIG, ipT630=INTBIG, ipT1304=INTBIG, ipSi10_606=INTBIG,
00048 ipT1039=INTBIG, ipT8446=INTBIG, ipT4368=INTBIG, ipTOI13=INTBIG,
00049 ipTOI11=INTBIG, ipTOI29=INTBIG, ipTOI46=INTBIG, ipTO1025=INTBIG,
00050 ipT304=INTBIG, ipT1214=INTBIG, ipT150=INTBIG, ipT146=INTBIG,
00051 ipT63=INTBIG, ipTO88=INTBIG, ipT52=INTBIG, ipT26=INTBIG,
00052 ipT1032=INTBIG, ipT1037=INTBIG, ipF0229=INTBIG, ipF0267=INTBIG,
00053 ipF444=INTBIG, ipF425=INTBIG, ipT770=INTBIG, ipT780=INTBIG,
00054 ipxNe0676=INTBIG, ipT895=INTBIG, ipT88=INTBIG, ipTNe13=INTBIG,
00055 ipTNe36=INTBIG, ipTNe16=INTBIG, ipTNe14=INTBIG, ipTNe24=INTBIG,
00056 ipT5895=INTBIG, ipfsNa373=INTBIG,ipfsNa490=INTBIG, ipfsNa421=INTBIG,
00057 ipxNa6143=INTBIG, ipxNa6862=INTBIG,ipxNa0746=INTBIG, ipMgI2853=INTBIG,
00058 ipMgI2026=INTBIG, ipT2796=INTBIG, ipT2804=INTBIG,
00059 ipT705=INTBIG, ipT4561=INTBIG, ipxMg51325=INTBIG, ipxMg52417=INTBIG,
00060 ipxMg52855=INTBIG,ipxMg71190=INTBIG, ipxMg72261=INTBIG,
00061 ipxMg72569=INTBIG,ipxMg08303=INTBIG, ipTMg610=INTBIG, ipTMg625=INTBIG,
00062 ipT58=INTBIG, ipTMg4=INTBIG, ipTMg14=INTBIG, ipTMg6=INTBIG,
00063 ipfsMg790=INTBIG, ipfsMg755=INTBIG, ipAlI3957=INTBIG, ipAlI3090=INTBIG,
00064 ipT1855=INTBIG, ipT1863=INTBIG, ipT2670=INTBIG,
00065 ipAl529=INTBIG, ipAl6366=INTBIG, ipAl6912=INTBIG, ipAl8575=INTBIG,
00066 ipAl8370=INTBIG, ipAl09204=INTBIG, ipT639=INTBIG,
00067 ipTAl550=INTBIG, ipTAl568=INTBIG, ipTAl48=INTBIG, ipSii2518=INTBIG,
00068 ipSii2215=INTBIG, ipT1808=INTBIG,
00069 ipT1207=INTBIG, ipT1895=INTBIG, ipT1394=INTBIG, ipT1403=INTBIG,
00070 ipT1527=INTBIG, ipT1305=INTBIG, ipT1260=INTBIG, ipSi619=INTBIG,
00071 ipSi10143=INTBIG, ipTSi499=INTBIG, ipTSi521=INTBIG, ipTSi41=INTBIG,
00072 ipTSi35=INTBIG, ipTSi25=INTBIG, ipTSi65=INTBIG,
00073 ipTSi3=INTBIG, ipTSi4=INTBIG, ipP0260=INTBIG, ipP0233=INTBIG,
00074 ipP0318=INTBIG, ipP713=INTBIG, ipP848=INTBIG, ipP817=INTBIG,
00075 ipP1027=INTBIG, ipP1018=INTBIG, ipT1256=INTBIG, ipT1194=INTBIG,
00076 ipTS1720=INTBIG, ipT1198=INTBIG, ipT786=INTBIG,
00077 ipT933=INTBIG, ipT944=INTBIG, ipfsS810=INTBIG, ipfsS912=INTBIG,
00078 ipfsS938=INTBIG, ipfsS1119=INTBIG, ipfsS1114=INTBIG, ipfsS1207=INTBIG,
00079 ipTSu418=INTBIG, ipTSu446=INTBIG, ipTSu30=INTBIG, ipTS19=INTBIG,
00080 ipTS34=INTBIG, ipTS11=INTBIG, ipfsCl214=INTBIG, ipfsCl233=INTBIG,
00081 ipCl04203=INTBIG, ipCl04117=INTBIG, ipCl973=INTBIG, ipCl1030=INTBIG,
00082 ipCl1092=INTBIG, ipT354=INTBIG, ipT389=INTBIG, ipT25=INTBIG,
00083 ipTAr7=INTBIG, ipTAr9=INTBIG, ipTAr22=INTBIG, ipTAr13=INTBIG,
00084 ipTAr8=INTBIG, ipAr06453=INTBIG, ipAr1055=INTBIG, ipAr1126=INTBIG,
00085 ipAr1178=INTBIG, ipKI7745=INTBIG, ipxK03462=INTBIG, ipxK04598=INTBIG,
00086 ipxK04154=INTBIG, ipxK06882=INTBIG, ipxK06557=INTBIG,
00087 ipxK07319=INTBIG, ipxK11425=INTBIG, ipCaI4228=INTBIG, ipT3934=INTBIG,
00088 ipT3969=INTBIG, ipT8498=INTBIG, ipT8542=INTBIG,
00089 ipT8662=INTBIG, ipT7291=INTBIG, ipT7324=INTBIG, ipTCa302=INTBIG,
00090 ipTCa345=INTBIG, ipTCa19=INTBIG, ipTCa3=INTBIG, ipTCa12=INTBIG,
00091 ipTCa4=INTBIG, ipCa0741=INTBIG, ipCa0761=INTBIG, ipCa08232=INTBIG,
00092 ipCa12333=INTBIG, ipSc05231=INTBIG, ipSc13264=INTBIG,
00093 ipTi06172=INTBIG, ipTi14212=INTBIG, ipVa07130=INTBIG, ipVa15172=INTBIG,
00094 ipCr08101=INTBIG, ipCr16141=INTBIG, ipxMn0979=INTBIG,
00095 ipxMn1712=INTBIG, ipFeI3884=INTBIG, ipFeI3729=INTBIG, ipFeI3457=INTBIG,
00096 ipFeI3021=INTBIG, ipFeI2966=INTBIG, ipTuv3=INTBIG,
00097 ipTr48=INTBIG, ipTFe16=INTBIG, ipTFe26=INTBIG, ipTFe34=INTBIG,
00098 ipTFe35=INTBIG, ipTFe46=INTBIG, ipTFe56=INTBIG, ipT1122=INTBIG,
00099 ipFe0795=INTBIG, ipFe0778=INTBIG, ipT245=INTBIG, ipT352=INTBIG,
00100 ipFe106375=INTBIG,ipT353=INTBIG,
00101 ipT347=INTBIG, ipT192=INTBIG, ipT255=INTBIG, ipT11=INTBIG,
00102 ipT191=INTBIG, ipFe18975=INTBIG, ipTFe23=INTBIG,
00103 ipTFe13=INTBIG, ipCo11527=INTBIG, ipxNi1242=INTBIG;
00104 long ipS4_1405=INTBIG,ipS4_1398=INTBIG,ipS4_1424=INTBIG,ipS4_1417=INTBIG,ipS4_1407=INTBIG,
00105 ipO4_1400=INTBIG,ipO4_1397=INTBIG,ipO4_1407=INTBIG,ipO4_1405=INTBIG,ipO4_1401=INTBIG,
00106 ipN3_1749=INTBIG,ipN3_1747=INTBIG,ipN3_1754=INTBIG,ipN3_1752=INTBIG,ipN3_1751=INTBIG,
00107 ipC2_2325=INTBIG,ipC2_2324=INTBIG,ipC2_2329=INTBIG,ipC2_2328=INTBIG,ipC2_2327=INTBIG,
00108 ipSi2_2334=INTBIG,ipSi2_2329=INTBIG,ipSi2_2350=INTBIG,ipSi2_2344=INTBIG,ipSi2_2336=INTBIG,
00109 ipFe22_247=INTBIG,ipFe22_217=INTBIG,ipFe22_348=INTBIG,ipFe22_292=INTBIG,
00110 ipFe22_253=INTBIG,ipFe22_846=INTBIG,ipTFe20_721=INTBIG,ipTFe20_578=INTBIG , ipZn04363=INTBIG,
00111 ipS12_520=INTBIG,ipS1_25m=INTBIG,ipS1_56m=INTBIG,ipCl1_11m=INTBIG,
00112 ipFe1_24m=INTBIG,ipFe1_35m=INTBIG,ipFe1_54m=INTBIG,ipFe1_111m=INTBIG,
00113 ipNi1_7m=INTBIG ,ipNi1_11m=INTBIG,ipSi1_130m=INTBIG,ipSi1_68m=INTBIG,
00114 ipNI_pumpDirect[NI_NDP]={INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG, INTBIG},
00115 ipNI_pumpIndirect=INTBIG, ipFe17_17=INTBIG;
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 #include "atomfeii.h"
00129
00130
00131 typedef enum { IS_NONE, IS_K_SHELL, IS_L1_SHELL, IS_L2_SHELL, IS_TOP } exc_type;
00132
00133 struct t_BadnellLevel
00134 {
00135 string config;
00136 int irsl;
00137 int S;
00138 int L;
00139 int g;
00140 realnum energy;
00141 bool lgAutoIonizing;
00142 exc_type WhichShell;
00143 t_BadnellLevel() : irsl(0), S(0), L(0), g(0), energy(0.f), lgAutoIonizing(false), WhichShell(IS_NONE) {}
00144 };
00145
00146
00147 STATIC void read_SH98_He1_cross_sections(void);
00148
00149 STATIC void read_Helike_cross_sections(void);
00150
00151
00152 STATIC void ReadBadnellAIData(const string& fnam,
00153 long nelem,
00154 long ion,
00155 vector<transition>& UTA,
00156 bitset<IS_TOP> Skip);
00157
00158 inline transition& InitTransition(transition& t);
00159 inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl);
00160
00161
00162
00163 const realnum gf_cutoff = 1.e-5f;
00164
00165 void atmdat_readin(void)
00166 {
00167 long int i,
00168 iCase ,
00169 ion,
00170 ipDens ,
00171 ipISO ,
00172 ipTemp ,
00173 j,
00174 ig0,
00175 ig1,
00176 imax,
00177 nelem,
00178 nelec,
00179 magic1,
00180 magic2,
00181 mol;
00182
00183 char cha;
00184 char chS2[3];
00185
00186 long ipZ;
00187 bool lgEOL;
00188
00189 FILE *ioDATA;
00190 char chLine[FILENAME_PATH_LENGTH_2] ,
00191 chFilename[FILENAME_PATH_LENGTH_2];
00192
00193 static bool lgFirstCall = true;
00194
00195 DEBUG_ENTRY( "atmdat_readin()" );
00196
00197
00198 if( !lgFirstCall )
00199 {
00200
00201 bool lgTooBig = false;
00202 for( j=0; j < iterations.iter_malloc; j++ )
00203 {
00204 if( geometry.nend[j]>=struc.nzlim )
00205 lgTooBig = true;
00206 }
00207 if( lgTooBig )
00208 {
00209 fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
00210 fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
00211 fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
00212 fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
00213 fprintf(ioQQQ," Sorry.\n");
00214 cdEXIT(EXIT_FAILURE);
00215 }
00216 return;
00217 }
00218
00219 lgFirstCall = false;
00220
00221
00222
00223 if( !mole.num_comole_calc )
00224 {
00225
00226 TotalInsanity();
00227 }
00228
00229
00230
00231
00232 for( j=0; j < iterations.iter_malloc; j++ )
00233 {
00234 struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
00235 }
00236
00237
00238 ++struc.nzlim;
00239
00240 struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
00241
00242 struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
00243
00244 struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00245
00246 struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00247
00248 struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00249
00250 struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00251
00252 struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00253
00254 struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00255
00256 struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00257
00258 struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00259
00260 struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00261
00262 struc.AccelTotalOutward = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00263
00264 struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00265
00266 struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00267
00268 struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00269
00270 struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00271
00272 struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00273
00274 struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00275
00276 struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00277
00278 struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00279
00280 struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00281
00282 struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00283
00284 struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00285
00286 struc.H2_molec = ((realnum**)MALLOC( (size_t)(N_H_MOLEC)*sizeof(realnum* )));
00287
00288 struc.CO_molec = ((realnum**)MALLOC( (size_t)(mole.num_comole_calc)*sizeof(realnum* )));
00289
00290 for(mol=0;mol<N_H_MOLEC;mol++)
00291 {
00292 struc.H2_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00293 }
00294
00295 for(mol=0;mol<mole.num_comole_calc;mol++)
00296 {
00297 struc.CO_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum )));
00298 }
00299
00300
00301 struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)LIMELM );
00302
00303 for( ipZ=0; ipZ< LIMELM;++ipZ )
00304 {
00305 struc.gas_phase[ipZ] =
00306 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
00307 }
00308
00309
00310 struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM+3) );
00311
00312 for( ipZ=0; ipZ< (LIMELM+3);++ipZ )
00313 {
00314
00315 struc.xIonDense[ipZ] =
00316 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+1) );
00317
00318 for( ion=0; ion < (LIMELM+1); ++ion )
00319 {
00320 struc.xIonDense[ipZ][ion] =
00321 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
00322 }
00323 }
00324
00325 struc.StatesElemNEW = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(LIMELM) );
00326 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00327 {
00328 if( dense.lgElmtOn[nelem] )
00329 {
00330 struc.StatesElemNEW[nelem] = (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(nelem+1) );
00331 for( long ion=0; ion<nelem+1; ion++ )
00332 {
00333 long ipISO = nelem-ion;
00334 if( ipISO < NISO )
00335 {
00336 struc.StatesElemNEW[nelem][ion] = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)iso.numLevels_max[ipISO][nelem]);
00337 for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
00338 {
00339 struc.StatesElemNEW[nelem][ion][level] =
00340 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
00341 }
00342 }
00343 else
00344 {
00345 fixit();
00346 struc.StatesElemNEW[nelem][ion] = NULL;
00347 }
00348 }
00349 }
00350 }
00351
00352
00353 for( i=0; i < struc.nzlim; i++ )
00354 {
00355 struc.testr[i] = 0.;
00356 struc.volstr[i] = 0.;
00357 struc.drad_x_fillfac[i] = 0.;
00358 struc.histr[i] = 0.;
00359 struc.hiistr[i] = 0.;
00360 struc.ednstr[i] = 0.;
00361 struc.o3str[i] = 0.;
00362 struc.heatstr[i] = 0.;
00363 struc.coolstr[i] = 0.;
00364 struc.pressure[i] = 0.;
00365 struc.pres_radiation_lines_curr[i] = 0.;
00366 struc.GasPressure[i] = 0.;
00367 struc.DenParticles[i] = 0.;
00368 struc.depth[i] = 0.;
00369 for(mol=0;mol<N_H_MOLEC;mol++)
00370 {
00371 struc.H2_molec[mol][i] = 0.;
00372 }
00373 for(mol=0;mol<mole.num_comole_calc;mol++)
00374 {
00375 struc.CO_molec[mol][i] = 0.;
00376 }
00377 }
00378
00379
00380 DynaCreateArrays( );
00381
00382
00383
00384
00385
00386
00387
00388 if( trace.lgTrace )
00389 fprintf( ioQQQ," atmdat_readin reading level1.dat\n");
00390
00391 ioDATA = open_data( "level1.dat", "r" );
00392
00393
00394 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00395 {
00396 fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
00397 cdEXIT(EXIT_FAILURE);
00398 }
00399
00400
00401
00402 nLevel1 = 0;
00403 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00404 {
00405
00406
00407 if( chLine[0] != '#')
00408 ++nLevel1;
00409 }
00410
00411
00412 TauLines = (transition *)MALLOC( (size_t)(nLevel1+1)*sizeof(transition ) );
00413
00414 for( i=0; i<nLevel1+1; i++ )
00415 {
00416 TauLines[i].Junk();
00417 TauLines[i].Hi = AddState2Stack();
00418 TauLines[i].Lo = AddState2Stack();
00419 TauLines[i].Emis = AddLine2Stack( true );
00420 }
00421
00422
00423 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
00424 {
00425 fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
00426 cdEXIT(EXIT_FAILURE);
00427 }
00428
00429
00430 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00431 {
00432 fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
00433 cdEXIT(EXIT_FAILURE);
00434 }
00435 i = 1;
00436
00437 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00438 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00439 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00440
00441
00442
00443 if( ( nelem != 10 ) || ( nelec != 3 ) || ( ion != 3 ) )
00444 {
00445 fprintf( ioQQQ,
00446 " atmdat_readin: the version of level1.dat is not the current version.\n" );
00447 fprintf( ioQQQ,
00448 " Please obtain the current version from the Cloudy web site.\n" );
00449 fprintf( ioQQQ,
00450 " I expected to find the number 09 11 18 and got %2.2li %2.2li %2.2li instead.\n" ,
00451 nelem , nelec , ion );
00452 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00453 cdEXIT(EXIT_FAILURE);
00454 }
00455
00456
00457
00458 i = 1;
00459 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00460 {
00461 if( i >= (nLevel1+1) )
00462 TotalInsanity();
00463
00464
00465 if( chLine[0] != '#')
00466 {
00467
00468
00469
00470
00471
00472
00473 strncpy( chS2 , chLine , 2);
00474 chS2[2] = 0;
00475
00476 ipZ = 0;
00477 for( j=0; j<LIMELM; ++j)
00478 {
00479 if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
00480 {
00481
00482 ipZ = j + 1;
00483 break;
00484 }
00485 }
00486
00487
00488 if( ipZ == 0 )
00489 {
00490 fprintf( ioQQQ,
00491 " atmdat_readin could not identify chem symbol on this level 1line:\n");
00492 fprintf( ioQQQ, "%s\n", chLine );
00493 fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
00494 cdEXIT(EXIT_FAILURE);
00495 }
00496
00497
00498 TauLines[i].Hi->nelem = (int)ipZ;
00499
00500
00501 j = 1;
00502 TauLines[i].Hi->IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00503 if( lgEOL )
00504 {
00505 fprintf( ioQQQ, " There should have been a number on this level1 line 1. Sorry.\n" );
00506 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00507 cdEXIT(EXIT_FAILURE);
00508 }
00509
00510 TauLines[i].Lo->nelem = TauLines[i].Hi->nelem;
00511 TauLines[i].Lo->IonStg = TauLines[i].Hi->IonStg;
00512
00513 TauLines[i].WLAng = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00514 if( lgEOL )
00515 {
00516 fprintf( ioQQQ, " There should have been a number on this level1 line 2. Sorry.\n" );
00517 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00518 cdEXIT(EXIT_FAILURE);
00519 }
00520
00521 TauLines[i].EnergyWN = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00522 if( lgEOL )
00523 {
00524 fprintf( ioQQQ, " There should have been a number on this level1 line 3. Sorry.\n" );
00525 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00526 cdEXIT(EXIT_FAILURE);
00527 }
00528
00529 TauLines[i].Lo->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00530 if( lgEOL )
00531 {
00532 fprintf( ioQQQ, " There should have been a number on this level1 line 4. Sorry.\n" );
00533 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00534 cdEXIT(EXIT_FAILURE);
00535 }
00536
00537 TauLines[i].Hi->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00538 if( lgEOL )
00539 {
00540 fprintf( ioQQQ, " There should have been a number on this level1 line 5. Sorry.\n" );
00541 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00542 cdEXIT(EXIT_FAILURE);
00543 }
00544
00545
00546 TauLines[i].Emis->gf = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00547 if( lgEOL )
00548 {
00549 fprintf( ioQQQ, " There should have been a number on this level1 line 6. Sorry.\n" );
00550 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00551 cdEXIT(EXIT_FAILURE);
00552 }
00553
00554 if( TauLines[i].Emis->gf < 0. )
00555 TauLines[i].Emis->gf = (realnum)pow((realnum)10.f,TauLines[i].Emis->gf);
00556
00557
00558 TauLines[i].Emis->Aul = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00559 if( lgEOL )
00560 {
00561 fprintf( ioQQQ, " There should have been a number on this level1 line 7. Sorry.\n" );
00562 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00563 cdEXIT(EXIT_FAILURE);
00564 }
00565 if( TauLines[i].Emis->Aul < 0. )
00566 TauLines[i].Emis->Aul = (realnum)pow((realnum)10.f,TauLines[i].Emis->Aul);
00567
00568 TauLines[i].Emis->iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00569 if( lgEOL )
00570 {
00571 fprintf( ioQQQ, " There should have been a number on this level1 line 8. Sorry.\n" );
00572 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00573 cdEXIT(EXIT_FAILURE);
00574 }
00575
00576
00577 TauLines[i].Coll.col_str = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00578
00579
00580 ++i;
00581 }
00582 }
00583
00584
00585 fclose(ioDATA);
00586 if( trace.lgTrace )
00587 fprintf( ioQQQ, " reading level1.dat OK\n" );
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 if( nWindLine > 0 )
00601 {
00602
00603
00604
00605
00606 TauLine2 = ((transition *)MALLOC( (size_t)nWindLine*sizeof(transition )));
00607 cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum )));
00608
00609
00610 for( i=0; i< nWindLine; ++i )
00611 {
00612
00613
00614 TauLine2[i].Junk();
00615
00616 TauLine2[i].Hi = AddState2Stack();
00617 TauLine2[i].Lo = AddState2Stack();
00618 TauLine2[i].Emis = AddLine2Stack( true );
00619 }
00620
00621 if( trace.lgTrace )
00622 fprintf( ioQQQ," atmdat_readin reading level2.dat\n");
00623
00624 ioDATA = open_data( "level2.dat", "r" );
00625
00626
00627 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00628 {
00629 fprintf( ioQQQ, " level2.dat error getting magic number\n" );
00630 cdEXIT(EXIT_FAILURE);
00631 }
00632
00633
00634 i = 1;
00635 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00636 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00637 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00638
00639
00640
00641 if( lgEOL || ( nelem != 9 ) || ( nelec != 11 ) || ( ion != 18 ) )
00642 {
00643 fprintf( ioQQQ,
00644 " atmdat_readin: the version of level2.dat is not the current version.\n" );
00645 fprintf( ioQQQ,
00646 " I expected to find the number 09 11 18 and got %2.2li %2.2li %2.2li instead.\n" ,
00647 nelem , nelec , ion );
00648 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00649 fprintf( ioQQQ, "Please obtain the correct version.\n");
00650 cdEXIT(EXIT_FAILURE);
00651 }
00652
00653
00654 i = 0;
00655 while( i < nWindLine )
00656 {
00657
00658 double tt[7];
00659 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00660 {
00661 fprintf( ioQQQ, " level2.dat error getting line %li\n", i );
00662 cdEXIT(EXIT_FAILURE);
00663 }
00664
00665 if( chLine[0]!='#' )
00666 {
00667
00668 sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " ,
00669 &tt[0] ,
00670 &tt[1] ,
00671 &tt[2] ,
00672 &tt[3] ,
00673 &tt[4] ,
00674 &tt[5] ,
00675 &tt[6] );
00676
00677
00678 TauLine2[i].Hi->nelem = (int)tt[0];
00679 TauLine2[i].Hi->IonStg = (int)tt[1];
00680 TauLine2[i].Lo->g = (realnum)tt[2];
00681 TauLine2[i].Hi->g = (realnum)tt[3];
00682 TauLine2[i].Emis->gf = (realnum)tt[4];
00683 TauLine2[i].EnergyWN = (realnum)tt[5];
00684 cs1_flag_lev2[i] = (realnum)tt[6];
00685 ++i;
00686 }
00687 }
00688
00689 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00690 {
00691 fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
00692 cdEXIT(EXIT_FAILURE);
00693 }
00694 sscanf( chLine , "%ld" , &magic2 );
00695 if( 999 != magic2 )
00696 {
00697 fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n",
00698 magic2 );
00699 cdEXIT(EXIT_FAILURE);
00700 }
00701 fclose( ioDATA );
00702 if( trace.lgTrace )
00703 fprintf( ioQQQ," reading level2.dat OK\n");
00704
00705
00706 }
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717 UTALines.reserve( 4000 );
00718
00719
00720 const char* chElmSymLC[] =
00721 { "h", "he", "li", "be", "b", "c", "n", "o", "f", "ne", "na", "mg", "al", "si", "p",
00722 "s", "cl", "ar", "k", "ca", "sc", "ti", "v", "cr", "mn", "fe", "co", "ni", "cu", "zn" };
00723
00724
00725 for( ipISO=ipLI_LIKE; ipISO < ipAL_LIKE; ++ipISO )
00726 {
00727 for( nelem=ipISO; nelem < LIMELM; ++nelem )
00728 {
00729
00730 long ion = nelem - ipISO;
00731 bitset<IS_TOP> Skip;
00732 if( ipISO < ipNA_LIKE )
00733 {
00734
00735 ostringstream oss;
00736
00737 oss << "nrb00#" << chElmSymLC[ipISO-1] << "_";
00738
00739 oss << chElmSymLC[nelem] << ion+1 << "ic1-2.dat";
00740
00741 ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
00742 }
00743 else
00744 {
00745
00746
00747 ostringstream oss;
00748 oss << "nrb00#" << chElmSymLC[ipISO-1] << "_";
00749 oss << chElmSymLC[nelem] << ion+1 << "ic1-3.dat";
00750 ReadBadnellAIData( oss.str(), nelem, ion, UTALines, Skip );
00751
00752
00753 if( ionbal.lgInnerShell_Kisielius && nelem == ipIRON &&
00754 ipISO >= ipNA_LIKE && ipISO <= ipAL_LIKE )
00755 Skip[IS_L2_SHELL] = 1;
00756
00757 ostringstream oss2;
00758 oss2 << "nrb00#" << chElmSymLC[ipISO-1] << "_";
00759 oss2 << chElmSymLC[nelem] << ion+1 << "ic2-3.dat";
00760 ReadBadnellAIData( oss2.str(), nelem, ion, UTALines, Skip );
00761 }
00762 }
00763 }
00764
00765
00766 const realnum StatWeightGroundLevelIron[] =
00767 { 9.f, 10.f, 9.f, 6.f, 1.f, 4.f, 5.f, 4.f, 1.f, 4.f, 5.f, 4.f, 1.f,
00768 2.f, 1.f, 2.f, 1.f, 4.f, 5.f, 4.f, 1.f, 2.f, 1.f, 2.f, 1.f, 2.f };
00769
00770
00771 transition BlankLine;
00772 BlankLine.Junk();
00773
00774
00775 if( ionbal.lgInnerShell_Gu06 )
00776 {
00777
00778
00779
00780 if( trace.lgTrace )
00781 fprintf( ioQQQ," atmdat_readin reading UTA_Gu06.dat\n");
00782
00783 FILE *ioGU06 = open_data( "UTA_Gu06.dat", "r" );
00784
00785 nelem = 0;
00786
00787
00788 while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
00789 {
00790
00791
00792 if( chLine[0] != '#')
00793 break;
00794 }
00795
00796 sscanf( chLine, "%li %li %li", &nelem, &nelec, &ion );
00797
00798
00799
00800 if( nelem != 2007 || nelec != 1 || ion != 23 )
00801 {
00802 fprintf( ioQQQ,
00803 " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
00804 fprintf( ioQQQ,
00805 " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
00806 nelem , nelec , ion );
00807 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00808 cdEXIT(EXIT_FAILURE);
00809 }
00810
00811 while( read_whole_line( chLine, (int)sizeof(chLine), ioGU06 ) != NULL )
00812 {
00813 if( chLine[0] != '#' )
00814 {
00815 long int i2;
00816 double WLAng, Aul, oscill, Aauto;
00817
00818 sscanf( chLine, "%4li%5li%8lf%13lf%12lf",
00819 &ion, &i2, &WLAng, &Aul, &Aauto );
00820 sscanf( &chLine[54], "%13lf", &oscill );
00821
00822
00823
00824
00825 ipISO = ipIRON - ion + 1;
00826 int ipThres = ionbal.lgInnerShell_Kisielius ? ipAL_LIKE : ipMG_LIKE;
00827 if( ipISO <= ipThres )
00828 continue;
00829
00830 UTALines.push_back( InitTransition( BlankLine ) );
00831
00832
00833 UTALines.back().Hi->nelem = ipIRON+1;
00834
00835
00836 ASSERT( ion > 0 && ion <= ipIRON );
00837
00838
00839 UTALines.back().Hi->IonStg = ion;
00840
00841
00842
00843 if( strstr_s( chLine, "(J=1/2)" ) != NULL )
00844 UTALines.back().Hi->g = 2.f;
00845 else if( strstr_s( chLine, "(J=1)" ) != NULL )
00846 UTALines.back().Hi->g = 3.f;
00847 else if( strstr_s( chLine, "(J=3/2)" ) != NULL )
00848 UTALines.back().Hi->g = 4.f;
00849 else if( strstr_s( chLine, "(J=2)" ) != NULL )
00850 UTALines.back().Hi->g = 5.f;
00851 else if( strstr_s( chLine, "(J=5/2)" ) != NULL )
00852 UTALines.back().Hi->g = 6.f;
00853 else if( strstr_s( chLine, "(J=3)" ) != NULL )
00854 UTALines.back().Hi->g = 7.f;
00855 else if( strstr_s( chLine, "(J=7/2)" ) != NULL )
00856 UTALines.back().Hi->g = 8.f;
00857 else if( strstr_s( chLine, "(J=4)" ) != NULL )
00858 UTALines.back().Hi->g = 9.f;
00859 else if( strstr_s( chLine, "(J=9/2)" ) != NULL )
00860 UTALines.back().Hi->g = 10.f;
00861 else if( strstr_s( chLine, "(J=5)" ) != NULL )
00862 UTALines.back().Hi->g = 11.f;
00863 else if( strstr_s( chLine, "(J=11/2)" ) != NULL )
00864 UTALines.back().Hi->g = 12.f;
00865 else
00866 TotalInsanity();
00867 UTALines.back().Lo->g = StatWeightGroundLevelIron[ion-1];
00868
00869
00870 UTALines.back().WLAng = (realnum)WLAng;
00871 UTALines.back().EnergyWN = (realnum)(1e8/WLAng);
00872 UTALines.back().EnergyErg = (realnum)(ERG1CM)*UTALines.back().EnergyWN;
00873
00874
00875 double frac_ioniz = Aauto/(Aul + Aauto);
00876 ASSERT( frac_ioniz >= 0. && frac_ioniz <= 1. );
00877 UTALines.back().Emis->AutoIonizFrac = (realnum)frac_ioniz;
00878
00879
00880 UTALines.back().Emis->gf = UTALines.back().Lo->g * (realnum)oscill;
00881
00882
00883
00884 UTALines.back().Emis->Aul =
00885 (realnum)eina( UTALines.back().Emis->gf,
00886 UTALines.back().EnergyWN,
00887 UTALines.back().Hi->g );
00888
00889 ASSERT( fp_equal_tol( (realnum)Aul, UTALines.back().Emis->Aul, 1.e-3f*(realnum)Aul ) );
00890
00891 UTALines.back().Emis->iRedisFun = ipPRD;
00892
00893
00894 if( UTALines.back().Emis->gf < gf_cutoff )
00895 UTALines.pop_back();
00896 }
00897 }
00898
00899 fclose( ioGU06 );
00900
00901 if( trace.lgTrace )
00902 fprintf( ioQQQ, " reading UTA_Gu06.dat OK\n" );
00903 }
00904 else
00905 {
00906
00907
00908 if( trace.lgTrace )
00909 fprintf( ioQQQ," atmdat_readin reading UTA_Behar.dat\n");
00910
00911 FILE *ioBEHAR = open_data( "UTA_Behar.dat", "r" );
00912
00913
00914 nelem = 0;
00915 if( read_whole_line( chLine, (int)sizeof(chLine), ioBEHAR ) != NULL )
00916 sscanf( chLine, "%li %li %li", &nelem, &nelec, &ion );
00917
00918
00919
00920 if( nelem != 2002 || nelec != 10 || ion != 22 )
00921 {
00922 fprintf( ioQQQ,
00923 " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" );
00924 fprintf( ioQQQ,
00925 " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
00926 nelem , nelec , ion );
00927 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00928 cdEXIT(EXIT_FAILURE);
00929 }
00930
00931
00932 while( read_whole_line( chLine, (int)sizeof(chLine), ioBEHAR ) != NULL )
00933 {
00934 if( chLine[0] != '#' )
00935 {
00936 long int i1, i2, i3;
00937 double f1, f2, oscill;
00938 double frac_relax;
00939
00940 sscanf( chLine ,"%li\t%li\t%li\t%lf\t%lf\t%lf\t%lf",
00941 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
00942
00943
00944
00945
00946 ipISO = ipIRON - i1;
00947 int ipThres = ionbal.lgInnerShell_Kisielius ? ipAL_LIKE : ipMG_LIKE;
00948 if( ipISO <= ipThres )
00949 continue;
00950
00951 UTALines.push_back( InitTransition( BlankLine ) );
00952
00953
00954 UTALines.back().Hi->nelem = ipIRON+1;
00955
00956
00957 ASSERT( i1 >= 0 && i1 < ipIRON );
00958
00959
00960 UTALines.back().Hi->IonStg = i1 + 1;
00961
00962
00963 UTALines.back().WLAng = (realnum)f1;
00964 UTALines.back().EnergyWN = 1.e8f/(realnum)f1;
00965 UTALines.back().EnergyErg = (realnum)ERG1CM*UTALines.back().EnergyWN;
00966
00967 UTALines.back().Lo->g = StatWeightGroundLevelIron[i1];
00968
00969
00970 double flu_test = GetGF( f2, UTALines.back().EnergyWN, 1. )/UTALines.back().Lo->g;
00971
00972 UTALines.back().Hi->g = (realnum)nint( oscill/flu_test );
00973
00974 ASSERT( i2 == 0 || fp_equal( UTALines.back().Hi->g, (realnum)i2 ) );
00975 ASSERT( i3 == 0 || fp_equal( UTALines.back().Lo->g, (realnum)i3 ) );
00976
00977
00978
00979 UTALines.back().Emis->Aul = (realnum)f2;
00980
00981
00982
00983 UTALines.back().Emis->gf = UTALines.back().Lo->g * (realnum)oscill;
00984
00985 UTALines.back().Emis->iRedisFun = ipPRD;
00986
00987
00988 ASSERT( frac_relax >= 0.f && frac_relax <= 1.f );
00989 UTALines.back().Emis->AutoIonizFrac = 1.f-(realnum)frac_relax;
00990
00991
00992 if( UTALines.back().Emis->gf < gf_cutoff )
00993 UTALines.pop_back();
00994 }
00995 }
00996
00997 fclose( ioBEHAR );
00998
00999 if( trace.lgTrace )
01000 fprintf( ioQQQ, " reading UTA_Behar.dat OK\n" );
01001 }
01002
01003 if( ionbal.lgInnerShell_Kisielius )
01004 {
01005
01006
01007
01008 if( trace.lgTrace )
01009 fprintf( ioQQQ," atmdat_readin reading UTA_Kisielius.dat\n");
01010
01011 FILE *ioROMAS = open_data( "UTA_Kisielius.dat", "r" );
01012
01013 while( read_whole_line( chLine, (int)sizeof(chLine), ioROMAS ) != NULL )
01014 {
01015
01016 if( chLine[0] != '#' )
01017 {
01018 long int i1, i2, i3;
01019 double f1, f2, oscill;
01020 double frac_relax;
01021
01022 sscanf( chLine, "%li\t%li\t%li\t%lf\t%lf\t%lf\t%lf",
01023 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
01024
01025
01026
01027 if( i2 == StatWeightGroundLevelIron[i1] )
01028 {
01029 UTALines.push_back( InitTransition( BlankLine ) );
01030
01031
01032 UTALines.back().Hi->nelem = ipIRON+1;
01033
01034
01035 ASSERT( i1 >= 0 && i1 < ipIRON );
01036
01037
01038 UTALines.back().Hi->IonStg = i1 + 1;
01039
01040
01041 UTALines.back().Hi->g = (realnum)i3;
01042 UTALines.back().Lo->g = (realnum)i2;
01043
01044 UTALines.back().WLAng = (realnum)f1;
01045 UTALines.back().EnergyWN = 1.e8f/(realnum)f1;
01046 UTALines.back().EnergyErg = (realnum)ERG1CM*UTALines.back().EnergyWN;
01047
01048
01049
01050 UTALines.back().Emis->gf = UTALines.back().Lo->g * (realnum)oscill;
01051 UTALines.back().Emis->Aul =
01052 (realnum)eina( UTALines.back().Emis->gf,
01053 UTALines.back().EnergyWN,
01054 UTALines.back().Hi->g );
01055
01056 ASSERT( fp_equal_tol( (realnum)f2, UTALines.back().Emis->Aul, 1.e-3f*(realnum)f2 ) );
01057
01058 UTALines.back().Emis->iRedisFun = ipPRD;
01059
01060
01061 ASSERT( frac_relax >= 0.f && frac_relax <= 1.f );
01062 UTALines.back().Emis->AutoIonizFrac = 1.f-(realnum)frac_relax;
01063
01064
01065 if( UTALines.back().Emis->gf < gf_cutoff )
01066 UTALines.pop_back();
01067 }
01068 }
01069 }
01070
01071 fclose( ioROMAS );
01072
01073 if( trace.lgTrace )
01074 fprintf( ioQQQ, " reading UTA_Kisielius.dat OK\n" );
01075 }
01076
01077 nUTA = (long)UTALines.size();
01078
01079
01080 if( false )
01081 {
01082 for( i=0; i < nUTA; ++i )
01083 dprintf( ioQQQ, "%5ld %s %2ld wavl %7.3f glo %2g gup %2g Aul %.2e gf %.2e ai branch %.3f\n",
01084 i,
01085 elementnames.chElementSym[UTALines[i].Hi->nelem-1],
01086 UTALines[i].Hi->IonStg,
01087 UTALines[i].WLAng,
01088 UTALines[i].Lo->g,
01089 UTALines[i].Hi->g,
01090 UTALines[i].Emis->Aul,
01091 UTALines[i].Emis->gf,
01092 UTALines[i].Emis->AutoIonizFrac );
01093 cdEXIT(EXIT_SUCCESS);
01094 }
01095
01096
01097
01098
01099
01100 HyperfineCreate();
01101
01102
01103 if( atmdat.lgLamdaOn || atmdat.lgChiantiOn )
01104 database_readin();
01105 else
01106 nSpecies = 0;
01107
01108
01109 lines_setup();
01110
01111
01112
01113
01114
01115
01116
01117 if( trace.lgTrace )
01118 fprintf( ioQQQ," atmdat_readin reading mewe_gbar.dat\n");
01119
01120 ioDATA = open_data( "mewe_gbar.dat", "r" );
01121
01122
01123 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01124 {
01125 fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" );
01126 cdEXIT(EXIT_FAILURE);
01127 }
01128
01129 sscanf( chLine , "%ld" , &magic1 );
01130 if( magic1 != 9101 )
01131 {
01132 fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n",
01133 magic1 );
01134 cdEXIT(EXIT_FAILURE);
01135 }
01136
01137
01138 for( i=1; i < 210; i++ )
01139 {
01140 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01141 {
01142 fprintf( ioQQQ, " mewe_gbar.dat error getting line %li\n", i );
01143 cdEXIT(EXIT_FAILURE);
01144 }
01145
01146 double help[4];
01147 sscanf( chLine, "%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] );
01148 for( int l=0; l < 4; ++l )
01149 MeweCoef.g[i][l] = (realnum)help[l];
01150 }
01151
01152
01153 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01154 {
01155 fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
01156 cdEXIT(EXIT_FAILURE);
01157 }
01158
01159 sscanf( chLine , "%ld" , &magic2 );
01160
01161 if( magic1 != magic2 )
01162 {
01163 fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n",
01164 magic2 );
01165 cdEXIT(EXIT_FAILURE);
01166 }
01167
01168 fclose( ioDATA );
01169
01170 if( trace.lgTrace )
01171 fprintf( ioQQQ," reading mewe_gbar.dat OK \n");
01172
01173
01174
01175
01176 for( nelem=0; nelem < LIMELM; nelem++ )
01177 for( ion=0; ion < LIMELM; ion++ )
01178 Heavy.nsShells[nelem][ion] = LONG_MAX;
01179
01180
01181
01182
01183
01184 for( nelem=2; nelem < LIMELM; nelem++ )
01185 {
01186
01187 for( ion=0; ion <= nelem; ion++ )
01188 {
01189
01190 nelec = nelem - ion + 1;
01191
01192
01193
01194
01195
01196
01197
01198 atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);
01199
01200 ASSERT( imax > 0 && imax <= 10 );
01201
01202
01203
01204 Heavy.nsShells[nelem][ion] = imax;
01205 }
01206 }
01207
01208
01209
01210
01211
01212
01213
01214 t_yield::Inst().init_yield();
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226 for( ipZ=0; ipZ<HS_NZ; ++ipZ )
01227 {
01228
01229 if( ipZ>1 && ipZ<5 ) continue;
01230
01231 for( iCase=0; iCase<2; ++iCase )
01232 {
01233
01234
01235
01236
01237
01238 sprintf( chFilename, "HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ? 'a' : 'b' );
01239
01240 if( trace.lgTrace )
01241 fprintf( ioQQQ," atmdat_readin reading Hummer Storey emission file %s\n",chFilename );
01242
01243
01244 ioDATA = open_data( chFilename, "r" );
01245
01246
01247 i = fscanf( ioDATA, "%li %li ",
01248 &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );
01249
01250
01251
01252 assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
01253 assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM);
01254
01255
01256 for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
01257 {
01258 for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
01259 {
01260 long int junk, junk2 , ne;
01261 i = fscanf( ioDATA, " %lf %li %lf %c %li %ld ",
01262 &atmdat.Density[iCase][ipZ][ipDens], &junk ,
01263 &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 ,
01264 &atmdat.ncut[iCase][ipZ] );
01265
01266 ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
01267 ASSERT( ne<=NLINEHS );
01268 for( j=0; j < ne; j++ )
01269 {
01270 i = fscanf( ioDATA, "%lf ",
01271 &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
01272 }
01273 }
01274 }
01275
01276
01277 fclose(ioDATA);
01278 if( trace.lgTrace )
01279 fprintf( ioQQQ," reading %s OK\n", chFilename );
01280
01281 # if 0
01282
01283 for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
01284 {
01285 fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]);
01286 }
01287 fprintf(ioQQQ,"\n");
01288 for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
01289 {
01290 fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]);
01291 }
01292 fprintf(ioQQQ,"\n");
01293 # endif
01294 }
01295 }
01296
01297
01298 read_SH98_He1_cross_sections();
01299
01300 read_Helike_cross_sections();
01301
01302
01303
01304
01305
01306
01307 FeIICreate();
01308 return;
01309 }
01310
01311 STATIC void read_SH98_He1_cross_sections(void)
01312 {
01313 DEBUG_ENTRY( "read_SH98_He1_cross_sections()" );
01314
01315 FILE *ioDATA;
01316 bool lgEOL;
01317
01318 char chPath[FILENAME_PATH_LENGTH_2],
01319 chDirectory[FILENAME_PATH_LENGTH_2],
01320 chLine[FILENAME_PATH_LENGTH_2];
01321
01322 const int ipNUM_FILES = 10;
01323
01324 char chFileNames[ipNUM_FILES][10] = {
01325 "p0202.3se",
01326 "p0202.3po",
01327 "p0202.3ge",
01328 "p0202.3fo",
01329 "p0202.3de",
01330 "p0202.1se",
01331 "p0202.1po",
01332 "p0202.1ge",
01333 "p0202.1fo",
01334 "p0202.1de" };
01335
01336 HS_He1_Xsectn = ((double****)MALLOC( (size_t)(26)*sizeof(double***)));
01337 HS_He1_Energy = ((double****)MALLOC( (size_t)(26)*sizeof(double***)));
01338
01339
01340 HS_He1_Xsectn[0] = NULL;
01341 HS_He1_Energy[0] = NULL;
01342
01343 for( long in = 1; in<=25; in++ )
01344 {
01345
01346 HS_He1_Xsectn[in] = ((double***)MALLOC( (size_t)(MIN2(5,in))*sizeof(double**)));
01347 HS_He1_Energy[in] = ((double***)MALLOC( (size_t)(MIN2(5,in))*sizeof(double**)));
01348 for( long il = 0; il<MIN2(5,in); il++ )
01349 {
01350
01351 HS_He1_Xsectn[in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
01352 HS_He1_Energy[in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
01353 HS_He1_Xsectn[in][il][0] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
01354 HS_He1_Energy[in][il][0] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
01355 HS_He1_Xsectn[in][il][1] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
01356 HS_He1_Energy[in][il][1] = ((double*)MALLOC( (size_t)(NUM_HS98_DATA_POINTS)*sizeof(double)));
01357 }
01358 }
01359
01360 # ifdef _WIN32
01361 strcpy( chDirectory, "sh98_he1\\pi\\" );
01362 # else
01363 strcpy( chDirectory, "sh98_he1/pi/" );
01364 # endif
01365
01366
01367
01368 for( long ipFile=0; ipFile<ipNUM_FILES; ipFile++ )
01369 {
01370 long S, L, index, N=0;
01371 long UNUSED P;
01372
01373 strcpy( chPath, chDirectory );
01374 strcat( chPath, chFileNames[ipFile] );
01375 ioDATA = open_data( chPath, "r" );
01376
01377 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01378 {
01379 long i=1, s;
01380 long i1, i2, i3;
01381 long numDataPoints;
01382
01383
01384 i1 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01385 i2 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01386 i3 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01387 if( i1==0 && i2==0 && i3==0 )
01388 break;
01389
01390
01391 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01392 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01393
01394 i=1;
01395
01396 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01397 S = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01398 L = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01399 P = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01400 index = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01401
01402
01403 ASSERT( index >= 1 );
01404
01405
01406
01407
01408 if( L==0 && S==3 )
01409 N = L + index + 1;
01410 else
01411 N = L + index;
01412
01413
01414 ASSERT( N<=25 );
01415
01416 if( S==1 )
01417 s=0;
01418 else if( S==3 )
01419 s=1;
01420 else
01421 TotalInsanity();
01422
01423 i=1;
01424
01425 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01426
01427 FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01428 numDataPoints = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01429
01430 ASSERT( numDataPoints == NUM_HS98_DATA_POINTS );
01431
01432
01433 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01434
01435
01436
01437
01438 for( long k=0; k<NUM_HS98_DATA_POINTS; k++ )
01439 {
01440 i=1;
01441 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01442 HS_He1_Energy[N][L][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01443 HS_He1_Xsectn[N][L][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01444 }
01445 }
01446
01447
01448 ASSERT( N==25 );
01449
01450 fclose( ioDATA );
01451 }
01452
01453 return;
01454 }
01455
01456 STATIC void read_Helike_cross_sections(void)
01457 {
01458 DEBUG_ENTRY( "read_Helike_cross_sections()" );
01459
01460 FILE *ioDATA;
01461 bool lgEOL;
01462
01463 char chLine[FILENAME_PATH_LENGTH_2];
01464 char chFileName[23] = "helike_pcs_topbase.dat";
01465
01466
01467 const int MaxN = 10;
01468 long last_i1=0;
01469
01470
01471
01472
01473 OP_Helike_Xsectn = ((double*****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(double****)));
01474 OP_Helike_Energy = ((double*****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(double****)));
01475 OP_Helike_NumPts = ((long****)MALLOC( (size_t)(ipCALCIUM+1)*sizeof(long***)));
01476
01477
01478 OP_Helike_Xsectn[ipHYDROGEN] = NULL;
01479 OP_Helike_Energy[ipHYDROGEN] = NULL;
01480 OP_Helike_NumPts[ipHYDROGEN] = 0;
01481 OP_Helike_Xsectn[ipHELIUM] = NULL;
01482 OP_Helike_Energy[ipHELIUM] = NULL;
01483 OP_Helike_NumPts[ipHELIUM] = 0;
01484
01485 for( long nelem = ipLITHIUM; nelem<= ipCALCIUM; nelem++ )
01486 {
01487
01488 OP_Helike_Xsectn[nelem] = ((double****)MALLOC( (size_t)(MaxN+1)*sizeof(double***)));
01489 OP_Helike_Energy[nelem] = ((double****)MALLOC( (size_t)(MaxN+1)*sizeof(double***)));
01490 OP_Helike_NumPts[nelem] = ((long***)MALLOC( (size_t)(MaxN+1)*sizeof(long**)));
01491
01492 for( long in = 1; in<=MaxN; in++ )
01493 {
01494
01495 OP_Helike_Xsectn[nelem][in] = ((double***)MALLOC( (size_t)(in)*sizeof(double**)));
01496 OP_Helike_Energy[nelem][in] = ((double***)MALLOC( (size_t)(in)*sizeof(double**)));
01497 OP_Helike_NumPts[nelem][in] = ((long**)MALLOC( (size_t)(in)*sizeof(long*)));
01498 for( long il = 0; il<in; il++ )
01499 {
01500
01501 OP_Helike_Xsectn[nelem][in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
01502 OP_Helike_Energy[nelem][in][il] = ((double**)MALLOC( (size_t)(2)*sizeof(double*)));
01503 OP_Helike_NumPts[nelem][in][il] = ((long*)MALLOC( (size_t)(2)*sizeof(long)));
01504
01505 OP_Helike_Xsectn[nelem][in][il][0] = NULL;
01506 OP_Helike_Energy[nelem][in][il][0] = NULL;
01507 OP_Helike_NumPts[nelem][in][il][0] = 0;
01508 OP_Helike_Xsectn[nelem][in][il][1] = NULL;
01509 OP_Helike_Energy[nelem][in][il][1] = NULL;
01510 OP_Helike_NumPts[nelem][in][il][1] = 0;
01511 }
01512 }
01513 }
01514
01515 ioDATA = open_data( chFileName, "r" );
01516
01517
01518
01519
01520
01521
01522
01523
01524 for( long i=0; i<3; i++)
01525 {
01526 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01527 {
01528 fprintf( ioQQQ,"PROBLEM corruption in TOPbase Helike pcs datafile.\nSorry\n" );
01529 cdEXIT(EXIT_FAILURE);
01530 }
01531 }
01532
01533 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01534 {
01535 long i=1;
01536 long n, l, s;
01537 long i1, i2, i3, i4, i5, i7;
01538 double i6;
01539
01540 i1 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01541 i2 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01542 i3 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01543 i4 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01544 i5 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01545 i6 = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01546 i7 = (long)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01547
01548 if( lgEOL )
01549 {
01550 fprintf( ioQQQ,"PROBLEM corruption in TOPbase Helike pcs datafile.\nSorry\n" );
01551 cdEXIT(EXIT_FAILURE);
01552 }
01553
01554
01555 if( i1==i2 && i1==i3 && i1==i4 && i1==i5 && i1==i7 && i1==-1 )
01556 break;
01557
01558
01559
01560 ASSERT( i1>0 && i1==(last_i1 + 1) && i1<=795 );
01561 last_i1 = i1;
01562
01563 ASSERT( (i2-1)<=ipCALCIUM && (i2-1)>=ipLITHIUM );
01564
01565 ASSERT( i3==2 );
01566
01567 ASSERT( i4>=100 && i4<400 );
01568 if( i4 >= 300 )
01569 s=1;
01570 else
01571 s=0;
01572 l = (i4 - (2*s+1)*100)/10;
01573
01574 ASSERT( l<=2 );
01575
01576 ASSERT( i5>=1 && i5<=10 );
01577 if( s==1 && l==0 )
01578 n = i5 + 1;
01579 else
01580 n = i5 + l;
01581 ASSERT( l<=MaxN );
01582
01583
01584 ASSERT( i6 < 0. );
01585
01586 OP_Helike_NumPts[i2-1][n][l][s] = i7;
01587 if( i7==0 )
01588 continue;
01589
01590 ASSERT( i7 > 0 );
01591 ASSERT( OP_Helike_Xsectn[i2-1][n][l][s]==NULL );
01592 ASSERT( OP_Helike_Energy[i2-1][n][l][s]==NULL );
01593 OP_Helike_Xsectn[i2-1][n][l][s] = ((double*)MALLOC( (size_t)(i7)*sizeof(double)));
01594 OP_Helike_Energy[i2-1][n][l][s] = ((double*)MALLOC( (size_t)(i7)*sizeof(double)));
01595
01596
01597
01598 for( long k=0; k<i7; k++ )
01599 {
01600 i=1;
01601 read_whole_line( chLine , (int)sizeof(chLine) , ioDATA );
01602 OP_Helike_Energy[i2-1][n][l][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01603 OP_Helike_Xsectn[i2-1][n][l][s][k] = (double)FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01604
01605
01606 if( k > 0 )
01607 {
01608 ASSERT( OP_Helike_Energy[i2-1][n][l][s][k] > OP_Helike_Energy[i2-1][n][l][s][k-1] );
01609 ASSERT( OP_Helike_Xsectn[i2-1][n][l][s][k] < OP_Helike_Xsectn[i2-1][n][l][s][k-1] );
01610 }
01611
01612
01613 FFmtRead( chLine, &i, sizeof(chLine), &lgEOL );
01614 ASSERT( lgEOL );
01615 }
01616 }
01617
01618 fclose( ioDATA );
01619
01620 return;
01621 }
01622
01623 t_yield::t_yield()
01624 {
01625
01626
01627
01628 for( int nelem=0; nelem < LIMELM; nelem++ )
01629 {
01630 for( int ion=0; ion < LIMELM; ion++ )
01631 {
01632 for( int ns=0; ns < 7; ns++ )
01633 {
01634 n_elec_eject[nelem][ion][ns] = LONG_MAX;
01635 for( int nelec=0; nelec < 10; nelec++ )
01636 {
01637 frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
01638 }
01639 }
01640 }
01641 }
01642
01643 lgKillAuger = false;
01644 }
01645
01646 void t_yield::init_yield()
01647 {
01648 char chLine[FILENAME_PATH_LENGTH_2];
01649 const char* chFilename;
01650
01651
01652 double temp[15];
01653
01654 DEBUG_ENTRY( "init_yield()" );
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665 ASSERT( Heavy.nsShells[2][0] > 0 );
01666
01667
01668 n_elec_eject[ipHYDROGEN][0][0] = 1;
01669 n_elec_eject[ipHELIUM][0][0] = 1;
01670 n_elec_eject[ipHELIUM][1][0] = 1;
01671
01672 frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
01673 frac_elec_eject[ipHELIUM][0][0][0] = 1;
01674 frac_elec_eject[ipHELIUM][1][0][0] = 1;
01675
01676
01677
01678 chFilename = "mewe_nelectron.dat";
01679
01680 if( trace.lgTrace )
01681 fprintf( ioQQQ, " init_yield reading %s\n", chFilename );
01682
01683 FILE *ioDATA;
01684
01685
01686 ioDATA = open_data( chFilename, "r" );
01687
01688
01689
01690
01691
01692 for( int nelem=2; nelem < LIMELM; nelem++ )
01693 {
01694
01695 for( int ion=0; ion <= nelem; ion++ )
01696 {
01697 for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
01698 {
01699 char ch1 = '#';
01700
01701 while( ch1 == '#' || ch1 == '*' )
01702 {
01703 if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
01704 {
01705 fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns );
01706 cdEXIT(EXIT_FAILURE);
01707 }
01708 ch1 = chLine[0];
01709 }
01710
01711 sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
01712 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
01713 &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
01714 &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] );
01715 n_elec_eject[nelem][ion][ns] = (long int)temp[3];
01716
01717 ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 );
01718
01719
01720 for( int j=0; j < 10; j++ )
01721 {
01722 frac_elec_eject[nelem][ion][ns][j] = (realnum)temp[j+4];
01723 ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. );
01724 }
01725 }
01726 }
01727
01728
01729 }
01730
01731 fclose( ioDATA );
01732
01733 if( trace.lgTrace )
01734 {
01735
01736 if( lgKillAuger )
01737 fprintf( ioQQQ, " Auger yields will be killed.\n");
01738 fprintf( ioQQQ, " reading %s OK\n", chFilename );
01739 }
01740
01741
01742
01743 chFilename = "mewe_fluor.dat";
01744
01745 if( trace.lgTrace )
01746 fprintf( ioQQQ, " init_yield reading %s\n", chFilename );
01747
01748
01749 ioDATA = open_data( chFilename, "r" );
01750
01751
01752
01753
01754 do
01755 {
01756 if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
01757 {
01758 fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 );
01759 cdEXIT(EXIT_FAILURE);
01760 }
01761 }
01762 while( chLine[0] == '#' );
01763
01764 bool lgEOL = false;
01765
01766 nfl_lines = 0;
01767 do
01768 {
01769 const int NKM = 10;
01770 int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
01771 int nAuger;
01772
01773 if( nfl_lines >= MEWE_FLUOR )
01774 TotalInsanity();
01775
01776
01777 sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf",
01778 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
01779 &temp[5], &temp[6] );
01780
01781
01782 nfl_nelem[nfl_lines] = (int)temp[0]-1;
01783 ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM );
01784
01785
01786 nfl_ion[nfl_lines] = (int)temp[1]-1;
01787 ASSERT( nfl_ion[nfl_lines] >= 0 && nfl_ion[nfl_lines] <= nfl_nelem[nfl_lines]+1 );
01788
01789
01790 nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1];
01791 ASSERT( nfl_nshell[nfl_lines] >= 0 &&
01792
01793 nfl_nshell[nfl_lines] < Heavy.nsShells[nfl_nelem[nfl_lines]][nfl_ion[nfl_lines]]-1 );
01794
01795
01796 nAuger = (int)temp[3];
01797
01798 nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1;
01799
01800 ASSERT( nfl_ion_emit[nfl_lines] > 0 && nfl_ion_emit[nfl_lines] <= nfl_nelem[nfl_lines]+1);
01801
01802
01803 nfl_nLine[nfl_lines] = (int)temp[4];
01804
01805
01806 fl_energy[nfl_lines] = (realnum)temp[5] / (realnum)EVRYD;
01807 ASSERT( fl_energy[nfl_lines] > 0. );
01808
01809
01810 fl_yield[nfl_lines] = (realnum)temp[6];
01811
01812 ASSERT( fl_yield[nfl_lines] >= 0 );
01813
01814 ++nfl_lines;
01815
01816 do
01817 {
01818 if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
01819 lgEOL = true;
01820 }
01821 while( chLine[0]=='#' && !lgEOL );
01822 }
01823 while( !lgEOL );
01824
01825 fclose( ioDATA );
01826 if( trace.lgTrace )
01827 fprintf( ioQQQ, " reading %s OK\n", chFilename );
01828 }
01829
01830
01831 STATIC void ReadBadnellAIData(const string& fnam,
01832 long nelem,
01833 long ion,
01834 vector<transition>& UTA,
01835 bitset<IS_TOP> Skip)
01836 {
01837 DEBUG_ENTRY( "ReadBadnellAIData()" );
01838
01839 if( trace.lgTrace )
01840 fprintf( ioQQQ," ReadBadnellAIData reading %s\n", fnam.c_str() );
01841
01842 fstream ioDATA;
01843 open_data( ioDATA, fnam.c_str(), mode_r );
01844
01845 string line;
01846 getline( ioDATA, line );
01847 ASSERT( line.substr(0,4) == "SEQ=" );
01848 getline( ioDATA, line );
01849 getline( ioDATA, line );
01850
01851 ASSERT( line.substr(3,21) == "PARENT LEVEL INDEXING" );
01852 int nParent;
01853 istringstream iss( line.substr(65,4) );
01854 iss >> nParent;
01855
01856 int nMulti = (nParent+5)/6;
01857 for( int i=0; i < nParent+5; ++i )
01858 getline( ioDATA, line );
01859
01860
01861 ASSERT( line.substr(3,26) == "IC RESOLVED LEVEL INDEXING" );
01862 int nLevel;
01863 istringstream iss2( line.substr(63,6) );
01864 iss2 >> nLevel;
01865
01866 for( int i=0; i < 3; ++i )
01867 getline( ioDATA, line );
01868
01869
01870 vector<t_BadnellLevel> level( nLevel );
01871 for( int i=0; i < nLevel; ++i )
01872 {
01873 getline( ioDATA, line );
01874 istringstream iss3( line );
01875 int indx, irsl;
01876 iss3 >> indx >> irsl;
01877 level[indx-1].irsl = irsl;
01878 level[indx-1].config = line.substr(16,20);
01879 istringstream iss4( line.substr(37,1) );
01880 iss4 >> level[indx-1].S;
01881 istringstream iss5( line.substr(39,1) );
01882 iss5 >> level[indx-1].L;
01883 istringstream iss6( line.substr(41,4) );
01884 double J;
01885 iss6 >> J;
01886 level[indx-1].g = nint(2.*J + 1.);
01887 istringstream iss7( line.substr(46,11) );
01888 iss7 >> level[indx-1].energy;
01889
01890
01891 level[indx-1].lgAutoIonizing = ( line[57] == '*' );
01892 if( level[indx-1].lgAutoIonizing )
01893 {
01894 if( level[indx-1].config.find( "1S1" ) != string::npos )
01895 level[indx-1].WhichShell = IS_K_SHELL;
01896 else if( level[indx-1].config.find( "2S1" ) != string::npos )
01897 level[indx-1].WhichShell = IS_L1_SHELL;
01898 else if( level[indx-1].config.find( "2P5" ) != string::npos )
01899 level[indx-1].WhichShell = IS_L2_SHELL;
01900 else
01901 TotalInsanity();
01902 }
01903 else
01904 {
01905 level[indx-1].WhichShell = IS_NONE;
01906 }
01907 }
01908
01909
01910
01911 while( getline( ioDATA, line ) )
01912 {
01913 if( line.find( "IRSL IRSL" ) != string::npos )
01914 break;
01915 }
01916
01917 for( int i=0; i < nMulti-1; ++i )
01918 getline( ioDATA, line );
01919
01920
01921 transition BlankLine;
01922 BlankLine.Junk();
01923
01924
01925 while( getline( ioDATA, line ) )
01926 {
01927
01928 if( line.size() < 10 )
01929 break;
01930
01931
01932
01933 if( line.size() < 50 )
01934 continue;
01935
01936
01937 line[19] = ' ';
01938
01939 int irsl_lo, irsl_hi, dum;
01940 double edif, Bij, Rji, Aai;
01941 istringstream iss8( line );
01942
01943
01944
01945
01946
01947
01948 iss8 >> irsl_lo >> irsl_hi >> dum >> dum >> edif >> Bij >> Rji >> Aai;
01949
01950 int ind_lo = irsl2ind( level, irsl_lo );
01951 int ind_hi = irsl2ind( level, irsl_hi );
01952 ASSERT( level[ind_hi].lgAutoIonizing );
01953
01954 for( int i=0; i < nMulti-1; ++i )
01955 getline( ioDATA, line );
01956
01957
01958 if( ind_lo == 0 && !Skip[level[ind_hi].WhichShell] )
01959 {
01960 UTA.push_back( InitTransition( BlankLine ) );
01961
01962
01963 UTA.back().Hi->nelem = nelem+1;
01964 UTA.back().Hi->IonStg = ion+1;
01965
01966 UTA.back().Hi->g = (realnum)level[ind_hi].g;
01967 UTA.back().Lo->g = (realnum)level[ind_lo].g;
01968
01969 double WavNum = edif*RYD_INF;
01970
01971
01972 UTA.back().WLAng = (realnum)(1e8/WavNum);
01973 UTA.back().EnergyWN = (realnum)WavNum;
01974 UTA.back().EnergyErg = (realnum)(ERG1CM)*UTA.back().EnergyWN;
01975
01976
01977 double frac_ioniz = Aai/(Rji + Aai);
01978 ASSERT( frac_ioniz >= 0. && frac_ioniz <= 1. );
01979 UTA.back().Emis->AutoIonizFrac = (realnum)frac_ioniz;
01980
01981
01982
01983
01984
01985 UTA.back().Emis->Aul = (realnum)(Bij*UTA.back().Lo->g/UTA.back().Hi->g);
01986 UTA.back().Emis->gf =
01987 (realnum)GetGF( UTA.back().Emis->Aul, UTA.back().EnergyWN, UTA.back().Hi->g );
01988
01989 UTA.back().Emis->iRedisFun = ipPRD;
01990
01991
01992 if( UTA.back().Emis->gf < gf_cutoff )
01993 UTA.pop_back();
01994 }
01995 }
01996
01997
01998 getline( ioDATA, line );
01999 ASSERT( line.substr(3,7) == "NRSLMX=" );
02000
02001 ioDATA.close();
02002
02003 if( trace.lgTrace )
02004 fprintf( ioQQQ, " reading %s OK\n", fnam.c_str() );
02005 }
02006
02007 inline transition& InitTransition(transition& t)
02008 {
02009 t.Hi = AddState2Stack();
02010 t.Lo = AddState2Stack();
02011 t.Emis = AddLine2Stack( true );
02012 return t;
02013 }
02014
02015 inline int irsl2ind(vector<t_BadnellLevel>& level, int irsl)
02016 {
02017 for( unsigned int i=0; i < level.size(); ++i )
02018 {
02019 if( level[i].irsl == irsl )
02020 return (int)i;
02021 }
02022
02023 TotalInsanity();
02024 }