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 "path.h"
00014 #include "trace.h"
00015 #include "lines.h"
00016 #include "lines_service.h"
00017 #include "ionbal.h"
00018 #include "struc.h"
00019 #include "geometry.h"
00020 #include "dynamics.h"
00021 #include "elementnames.h"
00022 #include "hyperfine.h"
00023 #include "atmdat.h"
00024
00025
00026
00027 #define INTBIG 2000000000
00028
00029
00030
00031
00032
00033
00034 long
00035 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
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126 #include "atomfeii.h"
00127
00128 void atmdat_readin(void)
00129 {
00130 long int i,
00131 iCase ,
00132 ion,
00133 ipDens ,
00134 ipISO ,
00135 ipTemp ,
00136 j,
00137 J,
00138 ig0,
00139 ig1,
00140 imax,
00141 nelem,
00142 nelec,
00143 magic1,
00144 magic2,
00145 mol ,
00146
00147 nUTA_Romas,
00148 nFeIonRomas[ipIRON],
00149
00150 nUTA_Behar ,
00151 nFeIonBehar[ipIRON],
00152
00153 nUTA_Gu06,
00154 nFeIonGu[ipIRON];
00155 double
00156 WLloRomas[ipIRON],
00157 WLhiRomas[ipIRON],
00158 WLloBehar[ipIRON],
00159 WLhiBehar[ipIRON],
00160 WLloGu[ipIRON],
00161 WLhiGu[ipIRON];
00162
00163 char cha;
00164 char chS2[3];
00165
00166 long ipZ;
00167 bool lgEOL;
00168
00169 FILE *ioDATA, *ioROMAS , *ioBEHAR=NULL;
00170 FILE *ioGU06=NULL;
00171 char chLine[FILENAME_PATH_LENGTH_2] ,
00172 chFilename[FILENAME_PATH_LENGTH_2];
00173
00174 static bool lgFirstCall = true;
00175
00176 DEBUG_ENTRY( "atmdat_readin()" );
00177
00178
00179 if( !lgFirstCall )
00180 {
00181
00182 bool lgTooBig = false;
00183 for( j=0; j < iterations.iter_malloc; j++ )
00184 {
00185 if( geometry.nend[j]>=struc.nzlim )
00186 lgTooBig = true;
00187 }
00188 if( lgTooBig )
00189 {
00190 fprintf(ioQQQ," This is the second or later calculation in a grid.\n");
00191 fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n");
00192 fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n");
00193 fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n");
00194 fprintf(ioQQQ," Sorry.\n");
00195 puts( "[Stop in atmdat_readin]" );
00196 cdEXIT(EXIT_FAILURE);
00197 }
00198
00199 DEBUG_EXIT( "atmdat_readin()" );
00200 return;
00201 }
00202
00203 lgFirstCall = false;
00204
00205
00206
00207 if( !mole.num_comole_calc )
00208 {
00209
00210 TotalInsanity();
00211 }
00212
00213
00214
00215 Badnell_rec_init();
00216
00217
00218
00219
00220 for( j=0; j < iterations.iter_malloc; j++ )
00221 {
00222 struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] );
00223 }
00224
00225
00226 ++struc.nzlim;
00227
00228 struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
00229
00230 struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double )));
00231
00232 struc.testr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00233
00234 struc.volstr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00235
00236 struc.drad_x_fillfac = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00237
00238 struc.histr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00239
00240 struc.hiistr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00241
00242 struc.ednstr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00243
00244 struc.o3str = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00245
00246 struc.pressure = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00247
00248 struc.pres_radiation_lines_curr = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00249
00250 struc.GasPressure = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00251
00252 struc.hden = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00253
00254 struc.DenParticles = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00255
00256 struc.DenMass = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00257
00258 struc.drad = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00259
00260 struc.depth = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00261
00262 struc.depth_last = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00263
00264 struc.drad_last = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00265
00266 struc.xLyman_depth = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00267
00268 struc.H2_molec = ((float**)MALLOC( (size_t)(N_H_MOLEC)*sizeof(float* )));
00269
00270 struc.CO_molec = ((float**)MALLOC( (size_t)(mole.num_comole_calc)*sizeof(float* )));
00271
00272 for(mol=0;mol<N_H_MOLEC;mol++)
00273 {
00274 struc.H2_molec[mol] = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00275 }
00276
00277 for(mol=0;mol<mole.num_comole_calc;mol++)
00278 {
00279 struc.CO_molec[mol] = ((float*)MALLOC( (size_t)(struc.nzlim)*sizeof(float )));
00280 }
00281
00282
00283 struc.gas_phase = (float **)MALLOC(sizeof(float *)*(unsigned)LIMELM );
00284
00285 for( ipZ=0; ipZ< LIMELM;++ipZ )
00286 {
00287 struc.gas_phase[ipZ] =
00288 (float*)MALLOC(sizeof(float)*(unsigned)(struc.nzlim) );
00289 }
00290
00291
00292 struc.xIonDense = (float ***)MALLOC(sizeof(float **)*(unsigned)(LIMELM+3) );
00293
00294 for( ipZ=0; ipZ< (LIMELM+3);++ipZ )
00295 {
00296
00297 struc.xIonDense[ipZ] =
00298 (float**)MALLOC(sizeof(float*)*(unsigned)(LIMELM+1) );
00299
00300 for( ion=0; ion < (LIMELM+1); ++ion )
00301 {
00302 struc.xIonDense[ipZ][ion] =
00303 (float*)MALLOC(sizeof(float)*(unsigned)(struc.nzlim) );
00304 }
00305 }
00306
00307
00308 for( i=0; i < struc.nzlim; i++ )
00309 {
00310 struc.testr[i] = 0.;
00311 struc.volstr[i] = 0.;
00312 struc.drad_x_fillfac[i] = 0.;
00313 struc.histr[i] = 0.;
00314 struc.hiistr[i] = 0.;
00315 struc.ednstr[i] = 0.;
00316 struc.o3str[i] = 0.;
00317 struc.heatstr[i] = 0.;
00318 struc.coolstr[i] = 0.;
00319 struc.pressure[i] = 0.;
00320 struc.pres_radiation_lines_curr[i] = 0.;
00321 struc.GasPressure[i] = 0.;
00322 struc.DenParticles[i] = 0.;
00323 struc.depth[i] = 0.;
00324 for(mol=0;mol<N_H_MOLEC;mol++)
00325 {
00326 struc.H2_molec[mol][i] = 0.;
00327 }
00328 for(mol=0;mol<mole.num_comole_calc;mol++)
00329 {
00330 struc.CO_molec[mol][i] = 0.;
00331 }
00332 }
00333
00334
00335 DynaCreateArrays( );
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 if( lgDataPathSet == true )
00346 {
00347
00348 strcpy( chFilename , chDataPath );
00349 strcat( chFilename , "level1.dat" );
00350 }
00351 else
00352 {
00353
00354 strcpy( chFilename , "level1.dat" );
00355 }
00356
00357 if( trace.lgTrace )
00358 fprintf( ioQQQ," atmdat_readin opening level1.dat:");
00359
00360 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00361 {
00362 fprintf( ioQQQ, " atmdat_readin could not open level1.dat\n" );
00363 if( lgDataPathSet == true )
00364 fprintf( ioQQQ, " even tried path\n" );
00365
00366 if( lgDataPathSet == true )
00367 {
00368 fprintf( ioQQQ, " atmdat_readin could not open level1.dat\n");
00369 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00370 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00371 }
00372
00373
00374 path_not_set();
00375
00376 puts( "[Stop in atmdat_readin]" );
00377 cdEXIT(EXIT_FAILURE);
00378 }
00379
00380
00381 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00382 {
00383 fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
00384 puts( "[Stop in atmdat_readin]" );
00385 cdEXIT(EXIT_FAILURE);
00386 }
00387
00388
00389
00390 nLevel1 = 0;
00391 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00392 {
00393
00394
00395 if( chLine[0] != '#')
00396 ++nLevel1;
00397 }
00398
00399
00400 TauLines = (EmLine *)MALLOC( (size_t)(nLevel1+1)*sizeof(EmLine ) );
00401
00402
00403 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
00404 {
00405 fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n");
00406 puts( "[Stop in atmdat_readin]" );
00407 cdEXIT(EXIT_FAILURE);
00408 }
00409
00410
00411 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00412 {
00413 fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n");
00414 puts( "[Stop in atmdat_readin]" );
00415 cdEXIT(EXIT_FAILURE);
00416 }
00417 i = 1;
00418
00419 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00420 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00421 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00422
00423
00424
00425 if( ( nelem != 5 ) || ( nelec != 12 ) || ( ion != 19 ) )
00426 {
00427 fprintf( ioQQQ,
00428 " atmdat_readin: the version of level1.dat is not the current version.\n" );
00429 fprintf( ioQQQ,
00430 " Please obtain the current version from the Cloudy web site.\n" );
00431 fprintf( ioQQQ,
00432 " I expected to find the number 04 11 5 and got %li %li %li instead.\n" ,
00433 nelem , nelec , ion );
00434 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00435 puts( "[Stop in atmdat_readin]" );
00436 cdEXIT(EXIT_FAILURE);
00437 }
00438
00439
00440
00441 i = 1;
00442 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00443 {
00444 if( i >= (nLevel1+1) )
00445 {
00446 fprintf(ioQQQ," version conflict.\n");
00447 }
00448
00449 if( chLine[0] != '#')
00450 {
00451
00452 EmLineJunk( &TauLines[i] );
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464 strncpy( chS2 , chLine , 2);
00465 chS2[2] = 0;
00466
00467 ipZ = 0;
00468 for( j=0; j<LIMELM; ++j)
00469 {
00470 if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 )
00471 {
00472
00473 ipZ = j + 1;
00474 break;
00475 }
00476 }
00477
00478
00479 if( ipZ == 0 )
00480 {
00481 fprintf( ioQQQ,
00482 " atmdat_readin could not identify chem symbol on this level 1line:\n");
00483 fprintf( ioQQQ, "%s\n", chLine );
00484 fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 );
00485 puts( "[Stop in atmdat_readin]" );
00486 cdEXIT(EXIT_FAILURE);
00487 }
00488
00489
00490 TauLines[i].nelem = (int)ipZ;
00491
00492
00493 j = 1;
00494 TauLines[i].IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00495 if( lgEOL )
00496 {
00497 fprintf( ioQQQ, " There should have been a number on this level1 line 1. Sorry.\n" );
00498 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00499 puts( "[Stop in atmdat_readin]" );
00500 cdEXIT(EXIT_FAILURE);
00501 }
00502
00503 TauLines[i].WLAng = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00504 if( lgEOL )
00505 {
00506 fprintf( ioQQQ, " There should have been a number on this level1 line 2. Sorry.\n" );
00507 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00508 puts( "[Stop in atmdat_readin]" );
00509 cdEXIT(EXIT_FAILURE);
00510 }
00511
00512 TauLines[i].EnergyWN = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00513 if( lgEOL )
00514 {
00515 fprintf( ioQQQ, " There should have been a number on this level1 line 3. Sorry.\n" );
00516 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00517 puts( "[Stop in atmdat_readin]" );
00518 cdEXIT(EXIT_FAILURE);
00519 }
00520
00521 TauLines[i].gLo = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00522 if( lgEOL )
00523 {
00524 fprintf( ioQQQ, " There should have been a number on this level1 line 4. Sorry.\n" );
00525 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00526 puts( "[Stop in atmdat_readin]" );
00527 cdEXIT(EXIT_FAILURE);
00528 }
00529
00530 TauLines[i].gHi = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00531 if( lgEOL )
00532 {
00533 fprintf( ioQQQ, " There should have been a number on this level1 line 5. Sorry.\n" );
00534 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00535 puts( "[Stop in atmdat_readin]" );
00536 cdEXIT(EXIT_FAILURE);
00537 }
00538
00539
00540
00541 TauLines[i].gf = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00542 if( lgEOL )
00543 {
00544 fprintf( ioQQQ, " There should have been a number on this level1 line 6. Sorry.\n" );
00545 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00546 puts( "[Stop in atmdat_readin]" );
00547 cdEXIT(EXIT_FAILURE);
00548 }
00549
00550 if( TauLines[i].gf < 0. )
00551 TauLines[i].gf = (float)pow(10.f,TauLines[i].gf);
00552
00553
00554 TauLines[i].Aul = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00555 if( lgEOL )
00556 {
00557 fprintf( ioQQQ, " There should have been a number on this level1 line 7. Sorry.\n" );
00558 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00559 puts( "[Stop in atmdat_readin]" );
00560 cdEXIT(EXIT_FAILURE);
00561 }
00562 if( TauLines[i].Aul < 0. )
00563 TauLines[i].Aul = (float)pow(10.f,TauLines[i].Aul);
00564
00565
00566 TauLines[i].iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00567 if( lgEOL )
00568 {
00569 fprintf( ioQQQ, " There should have been a number on this level1 line 8. Sorry.\n" );
00570 fprintf( ioQQQ, "string==%s==\n" ,chLine );
00571 puts( "[Stop in atmdat_readin]" );
00572 cdEXIT(EXIT_FAILURE);
00573 }
00574
00575
00576
00577 TauLines[i].cs = (float)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00578
00579
00580 ++i;
00581 }
00582 }
00583
00584
00585 fclose(ioDATA);
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597 if( nWindLine > 0 )
00598 {
00599
00600
00601
00602
00603 TauLine2 = ((EmLine *)MALLOC( (size_t)nWindLine*sizeof(EmLine )));
00604 cs1_flag_lev2 = ((float *)MALLOC( (size_t)nWindLine*sizeof(float )));
00605
00606
00607 for( i=0; i< nWindLine; ++i )
00608 {
00609
00610
00611 EmLineJunk( &TauLine2[i] );
00612 }
00613
00614
00615
00616 if( lgDataPathSet == true )
00617 {
00618
00619 strcpy( chFilename , chDataPath );
00620 strcat( chFilename , "level2.dat" );
00621 }
00622 else
00623 {
00624
00625 strcpy( chFilename , "level2.dat" );
00626 }
00627
00628 if( trace.lgTrace )
00629 fprintf( ioQQQ," atmdat_readin opening level2.dat:");
00630
00631 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00632 {
00633 fprintf( ioQQQ, " atmdat_readin could not open level2.dat\n" );
00634 if( lgDataPathSet == true )
00635 {
00636 fprintf( ioQQQ, " atmdat_readin could not open level2.dat, even tried path.\n");
00637 fprintf( ioQQQ, " path is *%s*\n",chDataPath );
00638 fprintf( ioQQQ, " final path is *%s*\n",chFilename );
00639 }
00640 puts( "[Stop in atmdat_readin]" );
00641 cdEXIT(EXIT_FAILURE);
00642 }
00643
00644
00645 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00646 {
00647 fprintf( ioQQQ, " level2.dat error getting magic number\n" );
00648 puts( "[Stop in atmdat_readin]" );
00649 cdEXIT(EXIT_FAILURE);
00650 }
00651
00652
00653 i = 1;
00654 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00655 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00656 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00657
00658
00659
00660
00661 if( lgEOL || ( nelem != 6 ) || ( nelec != 8 ) || ( ion != 10 ) )
00662 {
00663 fprintf( ioQQQ,
00664 " atmdat_readin: the version of level2.dat is not the current version.\n" );
00665 fprintf( ioQQQ,
00666 " I expected to find the number 06 08 10 and got %li %li %li instead.\n" ,
00667 nelem , nelec , ion );
00668 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00669 fprintf( ioQQQ, "Please obtain the correct version.\n");
00670 puts( "[Stop in atmdat_readin]" );
00671 cdEXIT(EXIT_FAILURE);
00672 }
00673
00674
00675 i = 0;
00676 while( i < nWindLine )
00677 {
00678 float tt[7];
00679 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00680 {
00681 fprintf( ioQQQ, " level2.dat error getting line %li\n", i );
00682 puts( "[Stop in atmdat_readin]" );
00683 cdEXIT(EXIT_FAILURE);
00684 }
00685
00686 if( chLine[0]!='#' )
00687 {
00688
00689 sscanf( chLine , "%f %f %f %f %f %f %f " ,
00690 &tt[0] ,
00691 &tt[1] ,
00692 &tt[2] ,
00693 &tt[3] ,
00694 &tt[4] ,
00695 &tt[5] ,
00696 &tt[6] );
00697
00698
00699 TauLine2[i].nelem = (int)tt[0];
00700 TauLine2[i].IonStg = (int)tt[1];
00701 TauLine2[i].gLo = tt[2];
00702 TauLine2[i].gHi = tt[3];
00703 TauLine2[i].gf = tt[4];
00704 TauLine2[i].EnergyWN = tt[5];
00705 cs1_flag_lev2[i] = tt[6];
00706 ++i;
00707 }
00708 }
00709
00710 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00711 {
00712 fprintf( ioQQQ, " level2.dat error getting last magic number\n" );
00713 puts( "[Stop in atmdat_readin]" );
00714 cdEXIT(EXIT_FAILURE);
00715 }
00716 sscanf( chLine , "%ld" , &magic2 );
00717 if( 999 != magic2 )
00718 {
00719 fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n",
00720 magic2 );
00721 puts( "[Stop in atmdat_readin]" );
00722 cdEXIT(EXIT_FAILURE);
00723 }
00724 fclose( ioDATA );
00725 if( trace.lgTrace )
00726 fprintf( ioQQQ," OK \n");
00727
00728
00729 }
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740 nUTA = 0;
00741 nUTA_Gu06 = 0;
00742 nUTA_Behar = 0;
00743 nUTA_Romas = 0;
00744
00745
00746
00747 if( ionbal.lgInnerShell_Gu06 )
00748 {
00749
00750
00751 if( lgDataPathSet == true )
00752 {
00753
00754 strcpy( chFilename , chDataPath );
00755 strcat( chFilename , "UTA_Gu06.dat" );
00756 }
00757 else
00758 {
00759
00760 strcpy( chFilename , "UTA_Gu06.dat" );
00761 }
00762
00763 if( trace.lgTrace )
00764 fprintf( ioQQQ," atmdat_readin opening UTA_Gu06.dat:");
00765
00766 if( ( ioGU06 = fopen( chFilename , "r" ) ) == NULL )
00767 {
00768 fprintf( ioQQQ, " atmdat_readin could not open UTA_Gu06.dat\n" );
00769 if( lgDataPathSet == true )
00770 fprintf( ioQQQ, " even tried path\n" );
00771
00772 if( lgDataPathSet == true )
00773 {
00774 fprintf( ioQQQ, " atmdat_readin could not open UTA_Gu06.dat\n");
00775 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00776 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00777 }
00778
00779 puts( "[Stop in atmdat_readin]" );
00780 cdEXIT(EXIT_FAILURE);
00781 }
00782
00783
00784 while( fgets( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
00785 {
00786
00787
00788 if( chLine[0] != '#')
00789 ++nUTA_Gu06;
00790 }
00791
00792 ASSERT( nUTA_Gu06 > 0 );
00793 --nUTA_Gu06;
00794
00795
00796 if( fseek( ioGU06 , 0 , SEEK_SET ) != 0 )
00797 {
00798 fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Gu06.dat.\n");
00799 puts( "[Stop in atmdat_readin]" );
00800 cdEXIT(EXIT_FAILURE);
00801 }
00802
00803
00804
00805
00806 while( fgets( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL )
00807 {
00808
00809
00810 if( chLine[0] != '#')
00811 break;
00812 }
00813
00814 j = 1;
00815 nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00816 nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00817 ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00818
00819
00820
00821 if( ( nelem != 2007 ) || ( nelec != 1 ) || ( ion != 23 ) )
00822 {
00823 fprintf( ioQQQ,
00824 " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" );
00825 fprintf( ioQQQ,
00826 " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" ,
00827 nelem , nelec , ion );
00828 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00829 puts( "[Stop in atmdat_readin]" );
00830 cdEXIT(EXIT_FAILURE);
00831 }
00832 }
00833 else
00834 {
00835
00836
00837
00838
00839
00840
00841
00842 if( lgDataPathSet == true )
00843 {
00844
00845 strcpy( chFilename , chDataPath );
00846 strcat( chFilename , "UTA_Behar.dat" );
00847 }
00848 else
00849 {
00850
00851 strcpy( chFilename , "UTA_Behar.dat" );
00852 }
00853
00854 if( trace.lgTrace )
00855 fprintf( ioQQQ," atmdat_readin opening UTA_Behar.dat:");
00856
00857 if( ( ioBEHAR = fopen( chFilename , "r" ) ) == NULL )
00858 {
00859 fprintf( ioQQQ, " atmdat_readin could not open UTA_Behar.dat\n" );
00860 if( lgDataPathSet == true )
00861 fprintf( ioQQQ, " even tried path\n" );
00862
00863 if( lgDataPathSet == true )
00864 {
00865 fprintf( ioQQQ, " atmdat_readin could not open UTA_Behar.dat\n");
00866 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00867 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00868 }
00869
00870 puts( "[Stop in atmdat_readin]" );
00871 cdEXIT(EXIT_FAILURE);
00872 }
00873
00874
00875
00876 if( fgets( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
00877 {
00878 fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n");
00879 puts( "[Stop in atmdat_readin]" );
00880 cdEXIT(EXIT_FAILURE);
00881 }
00882
00883 j = 1;
00884
00885 nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00886 nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00887 ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00888
00889
00890
00891
00892
00893 if( ( nelem != 2002 ) || ( nelec != 10 ) || ( ion != 22 ) )
00894 {
00895 fprintf( ioQQQ,
00896 " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" );
00897 fprintf( ioQQQ,
00898 " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" ,
00899 nelem , nelec , ion );
00900 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00901 puts( "[Stop in atmdat_readin]" );
00902 cdEXIT(EXIT_FAILURE);
00903 }
00904 }
00905
00906 if( trace.lgTrace )
00907 fprintf( ioQQQ," atmdat_readin UTA data files open OK\n");
00908
00909
00910
00911 nUTA_Behar = 0;
00912
00913
00914 for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
00915 {
00916 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00917 {
00918 ++nUTA_Behar;
00919 }
00920 }
00921
00922 if( !ionbal.lgInnerShell_Gu06 )
00923 {
00924
00925
00926 while( fgets( chLine , (int)sizeof(chLine) , ioBEHAR ) != NULL )
00927 {
00928
00929
00930 if( chLine[0] != '#')
00931 ++nUTA_Behar;
00932 }
00933
00934 if( fseek( ioBEHAR , 0 , SEEK_SET ) != 0 )
00935 {
00936 fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Behar.dat.\n");
00937 puts( "[Stop in atmdat_readin]" );
00938 cdEXIT(EXIT_FAILURE);
00939 }
00940
00941 if( fgets( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL )
00942 {
00943 fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n");
00944 puts( "[Stop in atmdat_readin]" );
00945 cdEXIT(EXIT_FAILURE);
00946 }
00947 }
00948
00949
00950
00951
00952 if( lgDataPathSet == true )
00953 {
00954
00955 strcpy( chFilename , chDataPath );
00956 strcat( chFilename , "UTA_Kisielius.dat" );
00957 }
00958 else
00959 {
00960
00961 strcpy( chFilename , "UTA_Kisielius.dat" );
00962 }
00963
00964 if( trace.lgTrace )
00965 fprintf( ioQQQ," atmdat_readin opening UTA_Kisielius.dat:");
00966 if( ( ioROMAS = fopen( chFilename , "r" ) ) == NULL )
00967 {
00968 fprintf( ioQQQ, " atmdat_readin could not open UTA_Kisielius.dat\n" );
00969 if( lgDataPathSet == true )
00970 fprintf( ioQQQ, " even tried path\n" );
00971
00972 if( lgDataPathSet == true )
00973 {
00974 fprintf( ioQQQ, " atmdat_readin could not open kisielius_uta.dat\n");
00975 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00976 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00977 }
00978
00979 puts( "[Stop in atmdat_readin]" );
00980 cdEXIT(EXIT_FAILURE);
00981 }
00982 nUTA_Romas = 0;
00983 while( fgets( chLine , (int)sizeof(chLine) , ioROMAS ) != NULL )
00984 {
00985
00986
00987 if( chLine[0] != '#')
00988 {
00989 ++nUTA_Romas;
00990 }
00991 }
00992
00993
00994
00995 UTALines = (EmLine *)MALLOC( (size_t)(nUTA_Behar+nUTA_Romas+nUTA_Gu06)*sizeof(EmLine) );
00996
00997
00998 for( ion=0; ion<ipIRON; ++ion )
00999 {
01000 nFeIonRomas[ion] = 0;
01001 WLloRomas[ion] = 1e10;
01002 WLhiRomas[ion] = 0.;
01003
01004 nFeIonBehar[ion] = 0;
01005 WLloBehar[ion] = 1e10;
01006 WLhiBehar[ion] = 0.;
01007
01008 nFeIonGu[ion] = 0;
01009 WLloGu[ion] = 1e10;
01010 WLhiGu[ion] = 0.;
01011 }
01012
01013
01014 i = 0;
01015
01016 if( fseek( ioROMAS , 0 , SEEK_SET ) != 0 )
01017 {
01018 fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Kisielius.dat.\n");
01019 puts( "[Stop in atmdat_readin]" );
01020 cdEXIT(EXIT_FAILURE);
01021 }
01022
01023 i = 0;
01024
01025
01026 for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO )
01027 {
01028 for( nelem=ipISO; nelem<LIMELM; ++nelem )
01029 {
01030 double f1;
01031 double sum_spon_auto;
01032 # define N 7
01033
01034 double alam[N]={ 1.859035 , 4.692138 , 10.324543 , 19.294364 ,
01035 40.401292 , 44.442754 , 72.959719 };
01036 double blam[N]={-9.7855968,-11.159739, -12.914931 , -14.987272,
01037 -18.853446 , -19.271781 ,-23.828388 };
01038 double clam[N]={ 8.2874628, 8.3043002, 8.3164038 , 8.3269937,
01039 8.4312895 , 8.3730893 , 8.493802 };
01040
01041
01042 double aA[N] = {1.9067 , 1.8715 , 2.1033, 2.2815 ,
01043 1.9511 , 1.9966 , 2.0852};
01044 double bA[N] = {8.4538 , 8.5528, 7.616, 6.9247 ,
01045 7.9479, 7.9777 , 7.8349 };
01046
01047 double aDep[N] = {29.54 , 12.07 , 24.23 , 7.35 , 7.37 , 11.14 , 11.99 };
01048 double bDep[N] = {-14.853 , -7.606 , -7.844 , 9.987 , 10.503 , 8.541 , 8.865 };
01049
01050
01051
01052 UTALines[i].nelem = nelem+1;
01053 ion = nelem + 1 - ipISO;
01054 UTALines[i].IonStg = ion;
01055
01056
01058 UTALines[i].gHi = 4.f;
01059 UTALines[i].gLo = 2.f;
01060
01061 f1 = alam[ipISO-2]*1e-4 + blam[ipISO-2]*1e-4*(nelem+1) +
01062 clam[ipISO-2]*1e-4*POW2(nelem+1.);
01063 f1 = (1./f1);
01064 UTALines[i].WLAng = (float)f1;
01065 UTALines[i].EnergyWN = (float)(1e8/f1);
01066 UTALines[i].EnergyErg = (float)(ERG1CM)*UTALines[i].EnergyWN;
01067
01068 f1 = aA[ipISO-2]*log((double)(nelem+1)) + bA[ipISO-2];
01069
01070 UTALines[i].Aul = (float)pow(10., f1 );
01071
01072 UTALines[i].gf =
01073 (float)(UTALines[i].Aul*UTALines[i].gHi*
01074 1.4992e-8*POW2(1e4/UTALines[i].EnergyWN));
01075
01076 UTALines[i].iRedisFun = 1;
01077
01078
01079 f1 = log((double)(nelem+1));
01080 f1 = POW2( f1 );
01081 sum_spon_auto = aDep[ipISO-2]*1e-2*f1 + bDep[ipISO-2]*1e-1;
01082
01083 if( ipISO == ipBE_LIKE )
01084 {
01085 sum_spon_auto = pow( E , sum_spon_auto );
01086 sum_spon_auto = pow(10., sum_spon_auto )*1e13;
01087 }
01088 else if( ipISO <= ipF_LIKE )
01089 {
01090 sum_spon_auto = pow(10., sum_spon_auto )*1e13;
01091 }
01092 else
01093 TotalInsanity();
01094
01095
01096
01097
01098
01099
01100 UTALines[i].cs = (float)MIN2(0., UTALines[i].Aul / sum_spon_auto - 1. );
01101
01102
01103 # if 0
01104 fprintf(ioQQQ," %2s%2li\t%.4f\t%.3e\t%.3e\t%.3e\n",
01105 elementnames.chElementSym[nelem],
01106 ion,
01107 UTALines[i].WLAng ,
01108 UTALines[i].Aul ,
01109 sum_spon_auto ,
01110 -UTALines[i].cs);
01111 # endif
01112
01113
01114 ++i;
01115 # undef N
01116 }
01117 }
01118
01119
01120 if( ionbal.lgInnerShell_Gu06 )
01121 {
01122 ioDATA = ioGU06;
01123 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01124 {
01125
01126 if( chLine[0] != '#')
01127 {
01128 long int i2;
01129 double WLAng, Aul, oscill , Aauto;
01130 EmLineJunk( &UTALines[i] );
01131
01132 sscanf( chLine ,"%4li%5li%8lf%13le%12le",
01133 &ion,&i2,&WLAng,&Aul,&Aauto);
01134 sscanf( &chLine[54] ,"%13le", &oscill);
01135
01136
01137 UTALines[i].nelem = ipIRON+1;
01138
01139
01140 ASSERT( ion >= 0 && ion < ipIRON );
01141
01142
01143 UTALines[i].IonStg = ion;
01144
01145
01146
01147
01148 UTALines[i].gHi = 1.f;
01149 UTALines[i].gLo = 1.f;
01150
01151
01152 UTALines[i].WLAng = (float)WLAng;
01153
01154
01155 ++nFeIonGu[ion-1];
01156 WLloGu[ion-1] = MIN2( WLAng , WLloGu[ion-1] );
01157 WLhiGu[ion-1] = MAX2( WLAng , WLloGu[ion-1] );
01158
01159
01160 UTALines[i].EnergyWN = (float)(1e8/WLAng);
01161 UTALines[i].EnergyErg = (float)(ERG1CM)*UTALines[i].EnergyWN;
01162
01163
01164
01165 UTALines[i].Aul = (float)Aul;
01166
01167
01168
01169
01170
01171 double frac_ioniz = Aauto/(Aul + Aauto);
01172 ASSERT( frac_ioniz >=0. && frac_ioniz <=1. );
01173 UTALines[i].cs = -(float)frac_ioniz;
01174
01175
01176 UTALines[i].gf = UTALines[i].gLo * (float)oscill;
01177
01178
01179
01180
01181
01182 UTALines[i].Aul = (float)eina(UTALines[i].gf,
01183 UTALines[i].EnergyWN,UTALines[i].gHi);
01184
01185 UTALines[i].iRedisFun = 1;
01186 {
01187
01188
01189 enum {DEBUG_LOC=false};
01190
01191 if( DEBUG_LOC && UTALines[i].nelem==ipIRON+1 && UTALines[i].IonStg==9)
01192 {
01193 fprintf(ioQQQ,"DEBUG Gu UTA\t%li\t%.3f\t%.3e\t%.3e\t%.3e\n",
01194 ion , WLAng, Aul ,
01195 eina(UTALines[i].gf,
01196 UTALines[i].EnergyWN,UTALines[i].gHi),
01197 Aauto );
01198 }
01199 }
01200
01201
01202 ++i;
01203 }
01204 }
01205 }
01206 else
01207 {
01208
01209 nUTA_Gu06 = 0;
01210
01211
01212
01213 ioDATA = ioBEHAR;
01214 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01215 {
01216
01217 if( chLine[0] != '#')
01218 {
01219 long int i1, i2, i3;
01220 double f1, f2, oscill;
01221 double frac_relax;
01222 EmLineJunk( &UTALines[i] );
01223
01224 sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
01225 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
01226
01227
01228 UTALines[i].nelem = ipIRON+1;
01229
01230
01231 ASSERT( i1 >= 0 && i1 < ipIRON );
01232
01233
01234 UTALines[i].IonStg = i1 + 1;
01235
01236
01237 if( i2>0 && i3>0 )
01238 {
01239 UTALines[i].gHi = (float)i2;
01240 UTALines[i].gLo = (float)i3;
01241 }
01242 else
01243 {
01244
01245
01246 UTALines[i].gHi = 1.f;
01247 UTALines[i].gLo = 1.f;
01248 }
01249
01250
01251 UTALines[i].WLAng = (float)f1;
01252
01253
01254 ++nFeIonBehar[i1];
01255 WLloBehar[i1] = MIN2( f1 , WLloBehar[i1] );
01256 WLhiBehar[i1] = MAX2( f1 , WLhiBehar[i1] );
01257
01258
01259 UTALines[i].EnergyWN = (float)(1e8/f1);
01260 UTALines[i].EnergyErg = (float)(ERG1CM)*UTALines[i].EnergyWN;
01261
01262
01263
01264 UTALines[i].Aul = (float)f2;
01265
01266
01267
01268 UTALines[i].gf = UTALines[i].gLo * (float)oscill;
01269
01270 UTALines[i].iRedisFun = 1;
01271
01272
01273
01274
01275 ASSERT( frac_relax >=0. && frac_relax <=1. );
01276 UTALines[i].cs = (float)(frac_relax-1.);
01277 # if 0
01278 fprintf(ioQQQ,"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n",
01279 nDataSet,
01280 i,
01281 elementnames.chElementSym[UTALines[i].nelem-1],
01282 UTALines[i].IonStg,
01283 UTALines[i].WLAng ,
01284 UTALines[i].Aul ,
01285 -UTALines[i].cs);
01286 # endif
01287
01288
01289 ++i;
01290 }
01291 }
01292 }
01293
01294
01295
01296
01297
01298
01299 ioDATA = ioROMAS;
01300 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01301 {
01302
01303 if( chLine[0] != '#')
01304 {
01305 long int i1, i2, i3;
01306 double f1, f2, oscill;
01307 double frac_relax;
01308 EmLineJunk( &UTALines[i] );
01309
01310 sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le",
01311 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill );
01312
01313
01314 UTALines[i].nelem = ipIRON+1;
01315
01316
01317 ASSERT( i1 >= 0 && i1 < ipIRON );
01318
01319
01320 UTALines[i].IonStg = i1 + 1;
01321
01322
01323 if( i2>0 && i3>0 )
01324 {
01325 UTALines[i].gHi = (float)i2;
01326 UTALines[i].gLo = (float)i3;
01327 }
01328 else
01329 {
01330
01331
01332 UTALines[i].gHi = 1.f;
01333 UTALines[i].gLo = 1.f;
01334 }
01335
01336
01337 UTALines[i].WLAng = (float)f1;
01338
01339
01340 ++nFeIonRomas[i1];
01341 WLloRomas[i1] = MIN2( f1 , WLloRomas[i1] );
01342 WLhiRomas[i1] = MAX2( f1 , WLhiRomas[i1] );
01343
01344 UTALines[i].EnergyWN = (float)(1e8/f1);
01345 UTALines[i].EnergyErg = (float)(ERG1CM)*UTALines[i].EnergyWN;
01346
01347
01348
01349 UTALines[i].Aul = (float)f2;
01350
01351
01352
01353 UTALines[i].gf = UTALines[i].gLo * (float)oscill;
01354
01355 UTALines[i].iRedisFun = 1;
01356
01357
01358
01359
01360 ASSERT( frac_relax >=0. && frac_relax <=1. );
01361 UTALines[i].cs = (float)(frac_relax-1.);
01362
01363
01364 ++i;
01365 }
01366 }
01367
01368
01369 ASSERT( i == (nUTA_Behar+nUTA_Romas+nUTA_Gu06) );
01370
01371 if( trace.lgTrace )
01372 {
01373 fprintf(ioQQQ , "summary of UTA data sets read in\nion numb WLlo WLhi\n");
01374 fprintf(ioQQQ , "Behar 01 total lines %li\n" , nUTA_Behar );
01375 for( ion=0; ion<ipIRON;++ion )
01376 {
01377 if( nFeIonBehar[ion] )
01378 {
01379 fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n",
01380 ion,
01381 nFeIonBehar[ion],
01382 WLloBehar[ion],
01383 WLhiBehar[ion]);
01384 }
01385 }
01386
01387 fprintf(ioQQQ , "Gu 06 total lines %li\n" , nUTA_Gu06 );
01388 for( ion=0; ion<ipIRON;++ion )
01389 {
01390 if( nFeIonGu[ion] )
01391 {
01392 fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n",
01393 ion,
01394 nFeIonGu[ion],
01395 WLloGu[ion],
01396 WLhiGu[ion]);
01397 }
01398 }
01399
01400 fprintf(ioQQQ , "Romas Kisielius 03 total lines %li\n", nUTA_Romas );
01401 for( ion=0; ion<ipIRON;++ion )
01402 {
01403 if( nFeIonRomas[ion] )
01404 {
01405 fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n",
01406 ion,
01407 nFeIonRomas[ion],
01408 WLloRomas[ion],
01409 WLhiRomas[ion] );
01410 }
01411 }
01412 }
01413
01414
01415
01416
01417
01418
01419
01420 if( ionbal.lgInnerShell_Kisielius )
01421 {
01422 for( i=0; i<nUTA_Behar+nUTA_Gu06; ++i )
01423 {
01424 if( UTALines[i].nelem-1 == ipIRON )
01425 {
01426
01427
01428
01429
01430 ASSERT( UTALines[i].IonStg-1< ipIRON );
01431 if( nFeIonRomas[UTALines[i].IonStg-1] )
01432 {
01433 UTALines[i].Aul = -1.;
01434 }
01435 }
01436 }
01437
01438 nUTA = nUTA_Behar + nUTA_Gu06 + nUTA_Romas;
01439 }
01440 else
01441 {
01442
01443
01444
01445 nUTA = nUTA_Behar + nUTA_Gu06;
01446 }
01447
01448
01449 fclose(ioROMAS);
01450 if( ionbal.lgInnerShell_Gu06 )
01451 fclose(ioGU06);
01452 else
01453 fclose(ioBEHAR);
01454
01455
01456
01457
01458 ASSERT( nCORotate > 0);
01459 C12O16Rotate = (EmLine *)MALLOC( (size_t)nCORotate*sizeof(EmLine) );
01460 C13O16Rotate = (EmLine *)MALLOC( (size_t)nCORotate*sizeof(EmLine) );
01461
01462
01463 lgCORotateMalloc = true;
01464
01465
01466 for( J=0; J< nCORotate; ++J )
01467 {
01468 EmLineJunk( &C12O16Rotate[J] );
01469 EmLineJunk( &C13O16Rotate[J] );
01470 }
01471
01472
01473
01474 HyperfineCreate();
01475
01476
01477 lines_setup();
01478
01479
01480
01481
01482
01483
01484 if( lgDataPathSet == true )
01485 {
01486
01487 strcpy( chFilename , chDataPath );
01488 strcat( chFilename , "mewe_gbar.dat" );
01489 }
01490 else
01491 {
01492
01493 strcpy( chFilename , "mewe_gbar.dat" );
01494 }
01495
01496 if( trace.lgTrace )
01497 fprintf( ioQQQ," atmdat_readin opening mewe_gbar.dat:");
01498
01499 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
01500 {
01501 fprintf( ioQQQ, " could not open mewe_gbar.dat\n" );
01502 if( lgDataPathSet == true )
01503 {
01504 fprintf( ioQQQ, " atmdat_readin could not open mewe_gbar.dat, even tried path.\n");
01505 fprintf( ioQQQ, " path is *%s*\n",chDataPath );
01506 fprintf( ioQQQ, " final path is *%s*\n",chFilename );
01507 }
01508 puts( "[Stop in atmdat_readin]" );
01509 cdEXIT(EXIT_FAILURE);
01510 }
01511
01512
01513 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01514 {
01515 fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" );
01516 puts( "[Stop in atmdat_readin]" );
01517 cdEXIT(EXIT_FAILURE);
01518 }
01519
01520 sscanf( chLine , "%ld" , &magic1 );
01521 if( magic1 != 9101 )
01522 {
01523 fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n",
01524 magic1 );
01525 puts( "[Stop in atmdat_readin]" );
01526 cdEXIT(EXIT_FAILURE);
01527 }
01528
01529
01530 for( i=1; i < 210; i++ )
01531 {
01532 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01533 {
01534 fprintf( ioQQQ, " mewe_gbar.dat error getting line %li\n", i );
01535 puts( "[Stop in atmdat_readin]" );
01536 cdEXIT(EXIT_FAILURE);
01537 }
01538
01539 sscanf( chLine , "%f %f %f %f " ,
01540 &MeweCoef.g[i][0] ,
01541 &MeweCoef.g[i][1] ,
01542 &MeweCoef.g[i][2] ,
01543 &MeweCoef.g[i][3] );
01544 }
01545
01546
01547 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01548 {
01549 fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" );
01550 puts( "[Stop in atmdat_readin]" );
01551 cdEXIT(EXIT_FAILURE);
01552 }
01553
01554 sscanf( chLine , "%ld" , &magic2 );
01555
01556 if( magic1 != magic2 )
01557 {
01558 fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n",
01559 magic2 );
01560 puts( "[Stop in atmdat_readin]" );
01561 cdEXIT(EXIT_FAILURE);
01562 }
01563
01564 fclose( ioDATA );
01565
01566 if( trace.lgTrace )
01567 fprintf( ioQQQ," OK \n");
01568
01569
01570
01571
01572 for( nelem=0; nelem < LIMELM; nelem++ )
01573 for( ion=0; ion < LIMELM; ion++ )
01574 Heavy.nsShells[nelem][ion] = LONG_MAX;
01575
01576
01577
01578
01579
01580 for( nelem=2; nelem < LIMELM; nelem++ )
01581 {
01582
01583 for( ion=0; ion <= nelem; ion++ )
01584 {
01585
01586 nelec = nelem - ion + 1;
01587
01588
01589
01590
01591
01592
01593
01594 atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1);
01595
01596 ASSERT( imax > 0 && imax <= 10 );
01597
01598
01599
01600 Heavy.nsShells[nelem][ion] = imax;
01601 }
01602 }
01603
01604
01605
01606
01607
01608
01609
01610 t_yield::Inst().init_yield();
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622 for( ipZ=0; ipZ<HS_NZ; ++ipZ )
01623 {
01624
01625 if( ipZ>1 && ipZ<5 ) continue;
01626
01627 for( iCase=0; iCase<2; ++iCase )
01628 {
01629
01630
01631 if( lgDataPathSet == true )
01632 {
01633
01634 strcpy( chFilename , chDataPath );
01635
01636 strcat( chFilename , "HS_e" );
01637 }
01638 else
01639 {
01640
01641 strcpy( chFilename , "e" );
01642 }
01643
01644
01645
01646
01647 chLine[0] = (char)(ipZ + '1');
01648
01649
01650 chLine[1] = '\0';
01651
01652
01653 strcat( chFilename , chLine );
01654
01655 if( iCase == 0 )
01656 {
01657
01658 strcat( chFilename , "a" );
01659 }
01660 else if(iCase == 1 )
01661 {
01662
01663 strcat( chFilename , "b" );
01664 }
01665
01666
01667 strcat( chFilename , ".dat");
01668
01669 if( trace.lgTrace )
01670 fprintf( ioQQQ," atmdat_readin opening Hummer Storey emission file %s:",chFilename);
01671
01672
01673 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
01674 {
01675 fprintf( ioQQQ, " could not open HS case b %s\n",chFilename);
01676 if( lgDataPathSet == true )
01677 {
01678 fprintf( ioQQQ, " atmdat_readin could not open even on path..\n");
01679 fprintf( ioQQQ, " path is *%s*\n",chDataPath );
01680 fprintf( ioQQQ, " final path is *%s*\n",chFilename );
01681 }
01682 puts( "[Stop in atmdat_readin]" );
01683 cdEXIT(EXIT_FAILURE);
01684 }
01685 if( trace.lgTrace )
01686 {
01687 fprintf( ioQQQ," OK\n");
01688 }
01689
01690
01691 i = fscanf( ioDATA, "%li %li ",
01692 &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] );
01693
01694
01695
01696 assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM );
01697 assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM);
01698
01699
01700 for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ )
01701 {
01702 for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ )
01703 {
01704 long int junk, junk2 , ne;
01705 fscanf( ioDATA, " %lf %li %lf %c %li %ld ",
01706 &atmdat.Density[iCase][ipZ][ipDens], &junk ,
01707 &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 ,
01708 &atmdat.ncut[iCase][ipZ] );
01709
01710 ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2;
01711 ASSERT( ne<=NLINEHS );
01712 for( j=0; j < ne; j++ )
01713 {
01714 fscanf( ioDATA, "%lf ", &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] );
01715 }
01716 }
01717 }
01718
01719
01720 fclose(ioDATA);
01721 # if 0
01722
01723 for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ )
01724 {
01725 fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]);
01726 }
01727 fprintf(ioQQQ,"\n");
01728 for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ )
01729 {
01730 fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]);
01731 }
01732 fprintf(ioQQQ,"\n");
01733 # endif
01734 }
01735 }
01736
01737
01738
01739
01740
01741
01742 FeIICreate( );
01743
01744 DEBUG_EXIT( "atmdat_readin()" );
01745 return;
01746 }
01747
01748 t_yield::t_yield()
01749 {
01750
01751
01752
01753 for( int nelem=0; nelem < LIMELM; nelem++ )
01754 {
01755 for( int ion=0; ion < LIMELM; ion++ )
01756 {
01757 for( int ns=0; ns < 7; ns++ )
01758 {
01759 n_elec_eject[nelem][ion][ns] = LONG_MAX;
01760 for( int nelec=0; nelec < 10; nelec++ )
01761 {
01762 frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX;
01763 }
01764 }
01765 }
01766 }
01767
01768 lgKillAuger = false;
01769 }
01770
01771 void t_yield::init_yield()
01772 {
01773 char chLine[FILENAME_PATH_LENGTH_2];
01774 char chFilename[FILENAME_PATH_LENGTH_2];
01775
01776 float temp[15];
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787 ASSERT( Heavy.nsShells[2][0] > 0 );
01788
01789
01790 n_elec_eject[ipHYDROGEN][0][0] = 1;
01791 n_elec_eject[ipHELIUM][0][0] = 1;
01792 n_elec_eject[ipHELIUM][1][0] = 1;
01793
01794 frac_elec_eject[ipHYDROGEN][0][0][0] = 1;
01795 frac_elec_eject[ipHELIUM][0][0][0] = 1;
01796 frac_elec_eject[ipHELIUM][1][0][0] = 1;
01797
01798
01799
01800 strcpy( chFilename, ( lgDataPathSet ? chDataPath : "" ) );
01801 strcat( chFilename, "mewe_nelectron.dat" );
01802
01803 if( trace.lgTrace )
01804 fprintf( ioQQQ, " init_yield opening %s:", chFilename );
01805
01806 FILE *ioDATA;
01807
01808
01809 if( ( ioDATA = fopen( chFilename, "r" ) ) == NULL )
01810 {
01811 fprintf( ioQQQ, " Could not open %s for reading\n", chFilename );
01812 fprintf( ioQQQ, " Is the path set correctly?\n" );
01813 puts( "[Stop in init_yield]" );
01814 cdEXIT(EXIT_FAILURE);
01815 }
01816
01817
01818
01819
01820
01821 for( int nelem=2; nelem < LIMELM; nelem++ )
01822 {
01823
01824 for( int ion=0; ion <= nelem; ion++ )
01825 {
01826 for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
01827 {
01828 char ch1 = '#';
01829
01830 while( ch1 == '#' || ch1 == '*' )
01831 {
01832 if( fgets( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
01833 {
01834 fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns );
01835 puts( "[Stop in init_yield]" );
01836 cdEXIT(EXIT_FAILURE);
01837 }
01838 ch1 = chLine[0];
01839 }
01840
01841 sscanf( chLine, "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",
01842 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
01843 &temp[5], &temp[6], &temp[7], &temp[8], &temp[9],
01844 &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] );
01845 n_elec_eject[nelem][ion][ns] = (long int)temp[3];
01846
01847 ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 );
01848
01849
01850 for( int j=0; j < 10; j++ )
01851 {
01852 frac_elec_eject[nelem][ion][ns][j] = temp[j+4];
01853 ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. );
01854 }
01855 }
01856 }
01857
01858
01859 }
01860
01861 fclose( ioDATA );
01862
01863
01864 if( lgKillAuger )
01865 do_kill_yield();
01866
01867 if( trace.lgTrace )
01868 fprintf( ioQQQ," OK \n");
01869
01870
01871
01872 strcpy( chFilename, ( lgDataPathSet ? chDataPath : "" ) );
01873 strcat( chFilename, "mewe_fluor.dat" );
01874
01875 if( trace.lgTrace )
01876 fprintf( ioQQQ, " init_yield opening %s:", chFilename );
01877
01878
01879 if( ( ioDATA = fopen( chFilename, "r" ) ) == NULL )
01880 {
01881 fprintf( ioQQQ, " Could not open %s for reading\n", chFilename );
01882 fprintf( ioQQQ, " Is the path set correctly?\n" );
01883 puts( "[Stop in init_yield]" );
01884 cdEXIT(EXIT_FAILURE);
01885 }
01886
01887
01888
01889
01890 do
01891 {
01892 if( fgets( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
01893 {
01894 fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 );
01895 puts( "[Stop in init_yield]" );
01896 cdEXIT(EXIT_FAILURE);
01897 }
01898 }
01899 while( chLine[0] == '#' );
01900
01901 bool lgEOL = false;
01902
01903 nfl_lines = 0;
01904 do
01905 {
01906 const int NKM = 10;
01907 int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 };
01908 int nAuger;
01909
01910 if( nfl_lines >= MEWE_FLUOR )
01911 TotalInsanity();
01912
01913
01914 sscanf( chLine, "%f %f %f %f %f %f %f",
01915 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4],
01916 &temp[5], &temp[6] );
01917
01918
01919 nfl_nelem[nfl_lines] = (int)temp[0]-1;
01920 ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM );
01921
01922
01923 nfl_ion[nfl_lines] = (int)temp[1]-1;
01924 ASSERT( nfl_ion[nfl_lines] >= 0 && nfl_ion[nfl_lines] <= nfl_nelem[nfl_lines]+1 );
01925
01926
01927 nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1];
01928 ASSERT( nfl_nshell[nfl_lines] >= 0 &&
01929
01930 nfl_nshell[nfl_lines] < Heavy.nsShells[nfl_nelem[nfl_lines]][nfl_ion[nfl_lines]]-1 );
01931
01932
01933 nAuger = (int)temp[3];
01934
01935 nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1;
01936
01937 ASSERT( nfl_ion_emit[nfl_lines] > 0 && nfl_ion_emit[nfl_lines] <= nfl_nelem[nfl_lines]+1);
01938
01939
01940 nfl_nLine[nfl_lines] = (int)temp[4];
01941
01942
01943 fl_energy[nfl_lines] = temp[5] / (float)EVRYD;
01944 ASSERT( fl_energy[nfl_lines] > 0. );
01945
01946
01947 fl_yield[nfl_lines] = temp[6];
01948
01949 ASSERT( fl_yield[nfl_lines] >= 0 );
01950
01951 ++nfl_lines;
01952
01953 do
01954 {
01955 if( fgets( chLine, (int)sizeof(chLine), ioDATA ) == NULL )
01956 lgEOL = true;
01957 }
01958 while( chLine[0]=='#' && !lgEOL );
01959 }
01960 while( !lgEOL );
01961
01962 fclose( ioDATA );
01963 }
01964
01965 void t_yield::do_kill_yield()
01966 {
01967
01968 ASSERT( lgKillAuger );
01969
01970 if( trace.lgTrace )
01971 fprintf( ioQQQ, " Auger data set zeroed.\n" );
01972 for( int nelem=0; nelem < LIMELM; nelem++ )
01973 {
01974 for( int nstage=0; nstage < nelem; nstage++ )
01975 {
01976 for( int nshell=0; nshell < 7; nshell++ )
01977 {
01978
01979 n_elec_eject[nelem][nstage][nshell] = 1;
01980 frac_elec_eject[nelem][nstage][nshell][0] = 1.;
01981 }
01982 }
01983 }
01984 }