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