cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*atmdat_readin read in some data files, but only if this is very first call, 00004 * called by Cloudy */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "taulines.h" 00008 #include "mewecoef.h" 00009 #include "iterations.h" 00010 #include "heavy.h" 00011 #include "mole.h" 00012 #include "yield.h" 00013 #include "trace.h" 00014 #include "lines.h" 00015 #include "lines_service.h" 00016 #include "ionbal.h" 00017 #include "struc.h" 00018 #include "geometry.h" 00019 #include "dynamics.h" 00020 #include "elementnames.h" 00021 #include "hyperfine.h" 00022 #include "atmdat.h" 00023 /* */ 00024 /* this was needed to get array to crash out of bounds if not set. 00025 * std limits on limits.h did not work with visual studio! */ 00026 #define INTBIG 2000000000 00027 00028 /* these are the individual pointers to the level 1 lines, they are set to 00029 * very large negative. 00030 * NB NB NB!! 00031 * these occur two times in the code!! 00032 * They are declared in taulines.h, and defined here */ 00033 long 00034 ipT1656=INTBIG, ipT9830=INTBIG, /*ipH21cm=INTBIG, ipHe3cm=INTBIG,*/ 00035 ipT8727=INTBIG, ipT1335=INTBIG, ipT1909=INTBIG, 00036 ipT977=INTBIG, ipT1550=INTBIG, ipT1548=INTBIG, ipT386=INTBIG, 00037 ipT310=INTBIG, ipc31175=INTBIG, ipT291=INTBIG, ipT280=INTBIG, 00038 ipT274=INTBIG, ipT270=INTBIG, ipT312=INTBIG, ipT610=INTBIG, 00039 ipT370=INTBIG, ipT157=INTBIG, ipT1085=INTBIG, 00040 ipT990=INTBIG, ipT1486=INTBIG, ipT765=INTBIG, ipT1243=INTBIG, 00041 ipT1239=INTBIG, ipT374g=INTBIG, ipT374x=INTBIG, ipT1200=INTBIG, 00042 ipT2140=INTBIG, ipT671=INTBIG, ipT315=INTBIG, ipT324=INTBIG, 00043 ipT333=INTBIG, ipT209=INTBIG, ipT122=INTBIG, ipT205=INTBIG, 00044 ipT57=INTBIG, ipT6300=INTBIG, ipT6363=INTBIG, ipT5577=INTBIG, 00045 ipT834=INTBIG, ipT1661=INTBIG, ipT1666=INTBIG, ipT835=INTBIG, 00046 ipT789=INTBIG, ipT630=INTBIG, ipT1304=INTBIG, ipSi10_606=INTBIG, 00047 ipT1039=INTBIG, ipT8446=INTBIG, ipT4368=INTBIG, ipTOI13=INTBIG, 00048 ipTOI11=INTBIG, ipTOI29=INTBIG, ipTOI46=INTBIG, ipTO1025=INTBIG, 00049 ipT304=INTBIG, ipT1214=INTBIG, ipT150=INTBIG, ipT146=INTBIG, 00050 ipT63=INTBIG, ipTO88=INTBIG, ipT52=INTBIG, ipT26=INTBIG, 00051 ipT1032=INTBIG, ipT1037=INTBIG, ipF0229=INTBIG, ipF0267=INTBIG, 00052 ipF444=INTBIG, ipF425=INTBIG, ipT770=INTBIG, ipT780=INTBIG, 00053 ipxNe0676=INTBIG, ipT895=INTBIG, ipT88=INTBIG, ipTNe13=INTBIG, 00054 ipTNe36=INTBIG, ipTNe16=INTBIG, ipTNe14=INTBIG, ipTNe24=INTBIG, 00055 ipT5895=INTBIG, ipfsNa373=INTBIG,ipfsNa490=INTBIG, ipfsNa421=INTBIG, 00056 ipxNa6143=INTBIG, ipxNa6862=INTBIG,ipxNa0746=INTBIG, ipMgI2853=INTBIG, 00057 ipMgI2026=INTBIG, ipT2796=INTBIG, ipT2804=INTBIG, 00058 ipT705=INTBIG, ipT4561=INTBIG, ipxMg51325=INTBIG, ipxMg52417=INTBIG, 00059 ipxMg52855=INTBIG,ipxMg71190=INTBIG, ipxMg72261=INTBIG, 00060 ipxMg72569=INTBIG,ipxMg08303=INTBIG, ipTMg610=INTBIG, ipTMg625=INTBIG, 00061 ipT58=INTBIG, ipTMg4=INTBIG, ipTMg14=INTBIG, ipTMg6=INTBIG, 00062 ipfsMg790=INTBIG, ipfsMg755=INTBIG, ipAlI3957=INTBIG, ipAlI3090=INTBIG, 00063 ipT1855=INTBIG, ipT1863=INTBIG, ipT2670=INTBIG, 00064 ipAl529=INTBIG, ipAl6366=INTBIG, ipAl6912=INTBIG, ipAl8575=INTBIG, 00065 ipAl8370=INTBIG, ipAl09204=INTBIG, ipT639=INTBIG, 00066 ipTAl550=INTBIG, ipTAl568=INTBIG, ipTAl48=INTBIG, ipSii2518=INTBIG, 00067 ipSii2215=INTBIG, ipT1808=INTBIG, 00068 ipT1207=INTBIG, ipT1895=INTBIG, ipT1394=INTBIG, ipT1403=INTBIG, 00069 ipT1527=INTBIG, ipT1305=INTBIG, ipT1260=INTBIG, ipSi619=INTBIG, 00070 ipSi10143=INTBIG, ipTSi499=INTBIG, ipTSi521=INTBIG, ipTSi41=INTBIG, 00071 ipTSi35=INTBIG, ipTSi25=INTBIG, ipTSi65=INTBIG, 00072 ipTSi3=INTBIG, ipTSi4=INTBIG, ipP0260=INTBIG, ipP0233=INTBIG, 00073 ipP0318=INTBIG, ipP713=INTBIG, ipP848=INTBIG, ipP817=INTBIG, 00074 ipP1027=INTBIG, ipP1018=INTBIG, ipT1256=INTBIG, ipT1194=INTBIG, 00075 ipTS1720=INTBIG, ipT1198=INTBIG, ipT786=INTBIG, 00076 ipT933=INTBIG, ipT944=INTBIG, ipfsS810=INTBIG, ipfsS912=INTBIG, 00077 ipfsS938=INTBIG, ipfsS1119=INTBIG, ipfsS1114=INTBIG, ipfsS1207=INTBIG, 00078 ipTSu418=INTBIG, ipTSu446=INTBIG, ipTSu30=INTBIG, ipTS19=INTBIG, 00079 ipTS34=INTBIG, ipTS11=INTBIG, ipfsCl214=INTBIG, ipfsCl233=INTBIG, 00080 ipCl04203=INTBIG, ipCl04117=INTBIG, ipCl973=INTBIG, ipCl1030=INTBIG, 00081 ipCl1092=INTBIG, ipT354=INTBIG, ipT389=INTBIG, ipT25=INTBIG, 00082 ipTAr7=INTBIG, ipTAr9=INTBIG, ipTAr22=INTBIG, ipTAr13=INTBIG, 00083 ipTAr8=INTBIG, ipAr06453=INTBIG, ipAr1055=INTBIG, ipAr1126=INTBIG, 00084 ipAr1178=INTBIG, ipKI7745=INTBIG, ipxK03462=INTBIG, ipxK04598=INTBIG, 00085 ipxK04154=INTBIG, ipxK06882=INTBIG, ipxK06557=INTBIG, 00086 ipxK07319=INTBIG, ipxK11425=INTBIG, ipCaI4228=INTBIG, ipT3934=INTBIG, 00087 ipT3969=INTBIG, ipT8498=INTBIG, ipT8542=INTBIG, 00088 ipT8662=INTBIG, ipT7291=INTBIG, ipT7324=INTBIG, ipTCa302=INTBIG, 00089 ipTCa345=INTBIG, ipTCa19=INTBIG, ipTCa3=INTBIG, ipTCa12=INTBIG, 00090 ipTCa4=INTBIG, ipCa0741=INTBIG, ipCa0761=INTBIG, ipCa08232=INTBIG, 00091 ipCa12333=INTBIG, ipSc05231=INTBIG, ipSc13264=INTBIG, 00092 ipTi06172=INTBIG, ipTi14212=INTBIG, ipVa07130=INTBIG, ipVa15172=INTBIG, 00093 ipCr08101=INTBIG, ipCr16141=INTBIG, ipxMn0979=INTBIG, 00094 ipxMn1712=INTBIG, ipFeI3884=INTBIG, ipFeI3729=INTBIG, ipFeI3457=INTBIG, 00095 ipFeI3021=INTBIG, ipFeI2966=INTBIG, ipTuv3=INTBIG, 00096 ipTr48=INTBIG, ipTFe16=INTBIG, ipTFe26=INTBIG, ipTFe34=INTBIG, 00097 ipTFe35=INTBIG, ipTFe46=INTBIG, ipTFe56=INTBIG, ipT1122=INTBIG, 00098 ipFe0795=INTBIG, ipFe0778=INTBIG, ipT245=INTBIG, ipT352=INTBIG, 00099 ipFe106375=INTBIG,ipT353=INTBIG, /*ipFe1310=INTBIG, ipFe1311=INTBIG,*/ 00100 ipT347=INTBIG, ipT192=INTBIG, ipT255=INTBIG, ipT11=INTBIG, 00101 ipT191=INTBIG, /*ipTFe07=INTBIG, ipTFe61=INTBIG,*/ ipFe18975=INTBIG, ipTFe23=INTBIG, 00102 ipTFe13=INTBIG, ipCo11527=INTBIG, ipxNi1242=INTBIG; 00103 long ipS4_1405=INTBIG,ipS4_1398=INTBIG,ipS4_1424=INTBIG,ipS4_1417=INTBIG,ipS4_1407=INTBIG, 00104 ipO4_1400=INTBIG,ipO4_1397=INTBIG,ipO4_1407=INTBIG,ipO4_1405=INTBIG,ipO4_1401=INTBIG, 00105 ipN3_1749=INTBIG,ipN3_1747=INTBIG,ipN3_1754=INTBIG,ipN3_1752=INTBIG,ipN3_1751=INTBIG, 00106 ipC2_2325=INTBIG,ipC2_2324=INTBIG,ipC2_2329=INTBIG,ipC2_2328=INTBIG,ipC2_2327=INTBIG, 00107 ipSi2_2334=INTBIG,ipSi2_2329=INTBIG,ipSi2_2350=INTBIG,ipSi2_2344=INTBIG,ipSi2_2336=INTBIG, 00108 ipFe22_247=INTBIG,ipFe22_217=INTBIG,ipFe22_348=INTBIG,ipFe22_292=INTBIG, 00109 ipFe22_253=INTBIG,ipFe22_846=INTBIG,ipTFe20_721=INTBIG,ipTFe20_578=INTBIG , ipZn04363=INTBIG, 00110 ipS12_520=INTBIG,ipS1_25m=INTBIG,ipS1_56m=INTBIG,ipCl1_11m=INTBIG, 00111 ipFe1_24m=INTBIG,ipFe1_35m=INTBIG,ipFe1_54m=INTBIG,ipFe1_111m=INTBIG, 00112 ipNi1_7m=INTBIG ,ipNi1_11m=INTBIG,ipSi1_130m=INTBIG,ipSi1_68m=INTBIG; 00113 00114 /* above are the individual pointers to the level 1 lines, they are set to 00115 * very large negative. 00116 * NB NB NB!! 00117 * these occur two times in the code!! 00118 * They are declared in taulines.h, and defined here */ 00119 00120 /* definition for whether level 2 lines are enabled, will be set to -1 00121 * with no level2 command */ 00122 /*long nWindLine = NWINDDIM;*/ 00123 /*realnum TauLine2[NWINDDIM][NTA];*/ 00124 /*realnum **TauLine2;*/ 00125 #include "atomfeii.h" 00126 00127 void atmdat_readin(void) 00128 { 00129 long int i, 00130 iCase , 00131 ion, 00132 ipDens , 00133 ipISO , 00134 ipTemp , 00135 j, 00136 J, 00137 ig0, 00138 ig1, 00139 imax, 00140 nelem, 00141 nelec, 00142 magic1, 00143 magic2, 00144 mol , 00145 00146 nUTA_Romas, 00147 nFeIonRomas[ipIRON], 00148 00149 nUTA_Behar , 00150 nFeIonBehar[ipIRON], 00151 00152 nUTA_Gu06, 00153 nFeIonGu[ipIRON]; 00154 double 00155 WLloRomas[ipIRON], 00156 WLhiRomas[ipIRON], 00157 WLloBehar[ipIRON], 00158 WLhiBehar[ipIRON], 00159 WLloGu[ipIRON], 00160 WLhiGu[ipIRON]; 00161 00162 char cha; 00163 char chS2[3]; 00164 00165 long ipZ; 00166 bool lgEOL; 00167 00168 FILE *ioDATA, *ioROMAS , *ioBEHAR=NULL; 00169 FILE *ioGU06=NULL; 00170 char chLine[FILENAME_PATH_LENGTH_2] , 00171 chFilename[FILENAME_PATH_LENGTH_2]; 00172 00173 static bool lgFirstCall = true; 00174 00175 DEBUG_ENTRY( "atmdat_readin()" ); 00176 00177 /* do nothing if not first call */ 00178 if( !lgFirstCall ) 00179 { 00180 /* do not do anything, but make sure that number of zones has not increased */ 00181 bool lgTooBig = false; 00182 for( j=0; j < iterations.iter_malloc; j++ ) 00183 { 00184 if( geometry.nend[j]>=struc.nzlim ) 00185 lgTooBig = true; 00186 } 00187 if( lgTooBig ) 00188 { 00189 fprintf(ioQQQ," This is the second or later calculation in a grid.\n"); 00190 fprintf(ioQQQ," The number of zones has been increased beyond what it was on the first calculation.\n"); 00191 fprintf(ioQQQ," This can\'t be done since space has already been allocated.\n"); 00192 fprintf(ioQQQ," Have the first calculation do the largest number of zones so that an increase is not needed.\n"); 00193 fprintf(ioQQQ," Sorry.\n"); 00194 cdEXIT(EXIT_FAILURE); 00195 } 00196 return; 00197 } 00198 00199 lgFirstCall = false; /* do not reevaluate again */ 00200 00201 /* make sure that molecules have been initialized - this will fail 00202 * if this routine is called before size of molecular network is known */ 00203 if( !mole.num_comole_calc ) 00204 { 00205 /* mole.num_comole_calc can't be zero */ 00206 TotalInsanity(); 00207 } 00208 00209 /* create space for the structure variables */ 00210 /* nzlim will be limit, and number allocated */ 00211 /* >>chng 01 jul 28, define this var, do all following mallocs */ 00212 for( j=0; j < iterations.iter_malloc; j++ ) 00213 { 00214 struc.nzlim = MAX2( struc.nzlim , geometry.nend[j] ); 00215 } 00216 00217 /* sloppy, but add one extra for safely */ 00218 ++struc.nzlim; 00219 00220 struc.coolstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double ))); 00221 00222 struc.heatstr = ((double*)MALLOC( (size_t)(struc.nzlim)*sizeof(double ))); 00223 00224 struc.testr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00225 00226 struc.volstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00227 00228 struc.drad_x_fillfac = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00229 00230 struc.histr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00231 00232 struc.hiistr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00233 00234 struc.ednstr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00235 00236 struc.o3str = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00237 00238 struc.pressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00239 00240 struc.windv = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00241 00242 struc.AccelTot = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00243 00244 struc.AccelGravity = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00245 00246 struc.pres_radiation_lines_curr = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00247 00248 struc.GasPressure = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00249 00250 struc.hden = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00251 00252 struc.DenParticles = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00253 00254 struc.DenMass = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00255 00256 struc.drad = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00257 00258 struc.depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00259 00260 struc.depth_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00261 00262 struc.drad_last = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00263 00264 struc.xLyman_depth = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00265 00266 struc.H2_molec = ((realnum**)MALLOC( (size_t)(N_H_MOLEC)*sizeof(realnum* ))); 00267 00268 struc.CO_molec = ((realnum**)MALLOC( (size_t)(mole.num_comole_calc)*sizeof(realnum* ))); 00269 00270 for(mol=0;mol<N_H_MOLEC;mol++) 00271 { 00272 struc.H2_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00273 } 00274 00275 for(mol=0;mol<mole.num_comole_calc;mol++) 00276 { 00277 struc.CO_molec[mol] = ((realnum*)MALLOC( (size_t)(struc.nzlim)*sizeof(realnum ))); 00278 } 00279 00280 /* create space for gas phase abundances array, first create space for the elements */ 00281 struc.gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)LIMELM ); 00282 /* last the zones */ 00283 for( ipZ=0; ipZ< LIMELM;++ipZ ) 00284 { 00285 struc.gas_phase[ipZ] = 00286 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) ); 00287 } 00288 00289 /* create space for struc.xIonDense array, first create space for the zones */ 00290 struc.xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(LIMELM+3) ); 00291 00292 for( ipZ=0; ipZ< (LIMELM+3);++ipZ ) 00293 { 00294 /* space for the ionization stage */ 00295 struc.xIonDense[ipZ] = 00296 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+1) ); 00297 /* now create diagonal of space for zones */ 00298 for( ion=0; ion < (LIMELM+1); ++ion ) 00299 { 00300 struc.xIonDense[ipZ][ion] = 00301 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) ); 00302 } 00303 } 00304 00305 /* some structure variables */ 00306 for( i=0; i < struc.nzlim; i++ ) 00307 { 00308 struc.testr[i] = 0.; 00309 struc.volstr[i] = 0.; 00310 struc.drad_x_fillfac[i] = 0.; 00311 struc.histr[i] = 0.; 00312 struc.hiistr[i] = 0.; 00313 struc.ednstr[i] = 0.; 00314 struc.o3str[i] = 0.; 00315 struc.heatstr[i] = 0.; 00316 struc.coolstr[i] = 0.; 00317 struc.pressure[i] = 0.; 00318 struc.pres_radiation_lines_curr[i] = 0.; 00319 struc.GasPressure[i] = 0.; 00320 struc.DenParticles[i] = 0.; 00321 struc.depth[i] = 0.; 00322 for(mol=0;mol<N_H_MOLEC;mol++) 00323 { 00324 struc.H2_molec[mol][i] = 0.; 00325 } 00326 for(mol=0;mol<mole.num_comole_calc;mol++) 00327 { 00328 struc.CO_molec[mol][i] = 0.; 00329 } 00330 } 00331 00332 /* allocate space for some arrays used by dynamics routines, and zero out vars */ 00333 DynaCreateArrays( ); 00334 00335 /************************************************************* 00336 * * 00337 * get the level 1 lines, with real collision data set * 00338 * * 00339 *************************************************************/ 00340 00341 if( trace.lgTrace ) 00342 fprintf( ioQQQ," atmdat_readin opening level1.dat:"); 00343 00344 ioDATA = open_data( "level1.dat", "r" ); 00345 00346 /* first line is a version number and does not count */ 00347 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00348 { 00349 fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n"); 00350 cdEXIT(EXIT_FAILURE); 00351 } 00352 00353 /* count how many lines are in the file, ignoring all lines 00354 * starting with '#' */ 00355 nLevel1 = 0; 00356 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00357 { 00358 /* we want to count the lines that do not start with # 00359 * since these contain data */ 00360 if( chLine[0] != '#') 00361 ++nLevel1; 00362 } 00363 00364 /* now malloc the TauLines array */ 00365 TauLines = (transition *)MALLOC( (size_t)(nLevel1+1)*sizeof(transition ) ); 00366 00367 for( i=0; i<nLevel1+1; i++ ) 00368 { 00369 TransitionJunk( &TauLines[i] ); 00370 TauLines[i].Hi = AddState2Stack(); 00371 TauLines[i].Lo = AddState2Stack(); 00372 TauLines[i].Emis = AddLine2Stack( true ); 00373 } 00374 00375 /* now rewind the file so we can read it a second time*/ 00376 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 ) 00377 { 00378 fprintf( ioQQQ, " atmdat_readin could not rewind level1.dat.\n"); 00379 cdEXIT(EXIT_FAILURE); 00380 } 00381 00382 /* check that magic number is ok */ 00383 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00384 { 00385 fprintf( ioQQQ, " atmdat_readin could not read first line of level1.dat.\n"); 00386 cdEXIT(EXIT_FAILURE); 00387 } 00388 i = 1; 00389 /* level 1 magic number */ 00390 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00391 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00392 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00393 00394 /* magic number 00395 * the following is the set of numbers that appear at the start of level1.dat */ 00396 if( ( nelem != 5 ) || ( nelec != 12 ) || ( ion != 19 ) ) 00397 { 00398 fprintf( ioQQQ, 00399 " atmdat_readin: the version of level1.dat is not the current version.\n" ); 00400 fprintf( ioQQQ, 00401 " Please obtain the current version from the Cloudy web site.\n" ); 00402 fprintf( ioQQQ, 00403 " I expected to find the number 04 11 5 and got %li %li %li instead.\n" , 00404 nelem , nelec , ion ); 00405 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00406 cdEXIT(EXIT_FAILURE); 00407 } 00408 00409 /* this starts at 1 not 0 since zero is reserved for the 00410 * dummy line - we counted number of lines above */ 00411 i = 1; 00412 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00413 { 00414 if( i >= (nLevel1+1) ) 00415 { 00416 fprintf(ioQQQ," version conflict.\n"); 00417 } 00418 /* only look at lines without '#' in first col */ 00419 if( chLine[0] != '#') 00420 { 00421 /* >>chng 00 apr 24, prevents crash in PresTotCurrent on call from ContSetIntensity, by PvH */ 00422 00423 00424 /*sscanf( chLine , 00425 "%2s%2li%9li %11f %2li %2li %10g %10g %2li\n", 00426 chS2,&i2,&i3,&f1,&i4,&i5,&f2,&f3,&i6 );*/ 00427 00428 /* first two cols of line has chem element symbol, 00429 * try to match it with the list of names 00430 * there are no H lines in the level 1 set so zero 00431 * is an invalid result */ 00432 00433 /* copy first two char into chS2 and null terminate */ 00434 strncpy( chS2 , chLine , 2); 00435 chS2[2] = 0; 00436 00437 ipZ = 0; 00438 for( j=0; j<LIMELM; ++j) 00439 { 00440 if( strcmp( elementnames.chElementSym[j], chS2 ) == 0 ) 00441 { 00442 /* ipZ is the actual atomic number starting from 1 */ 00443 ipZ = j + 1; 00444 break; 00445 } 00446 } 00447 00448 /* this happens if no valid chemical element in first two cols on this line */ 00449 if( ipZ == 0 ) 00450 { 00451 fprintf( ioQQQ, 00452 " atmdat_readin could not identify chem symbol on this level 1line:\n"); 00453 fprintf( ioQQQ, "%s\n", chLine ); 00454 fprintf( ioQQQ, "looking for this string==%2s==\n",chS2 ); 00455 cdEXIT(EXIT_FAILURE); 00456 } 00457 00458 /* now stuff them into the array, will convert over to proper units later */ 00459 TauLines[i].Hi->nelem = (int)ipZ; 00460 /* start the scan early on the line -- the element label will not be 00461 * picked up, first number will be the ion stage after the label, as in c 4 */ 00462 j = 1; 00463 TauLines[i].Hi->IonStg = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00464 if( lgEOL ) 00465 { 00466 fprintf( ioQQQ, " There should have been a number on this level1 line 1. Sorry.\n" ); 00467 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00468 cdEXIT(EXIT_FAILURE); 00469 } 00470 00471 TauLines[i].Lo->nelem = TauLines[i].Hi->nelem; 00472 TauLines[i].Lo->IonStg = TauLines[i].Hi->IonStg; 00473 00474 TauLines[i].WLAng = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00475 if( lgEOL ) 00476 { 00477 fprintf( ioQQQ, " There should have been a number on this level1 line 2. Sorry.\n" ); 00478 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00479 cdEXIT(EXIT_FAILURE); 00480 } 00481 /*TauLines[i].WLAng = (realnum)i3;*/ 00482 TauLines[i].EnergyWN = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00483 if( lgEOL ) 00484 { 00485 fprintf( ioQQQ, " There should have been a number on this level1 line 3. Sorry.\n" ); 00486 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00487 cdEXIT(EXIT_FAILURE); 00488 } 00489 /*TauLines[i].EnergyWN = (realnum)f1;*/ 00490 TauLines[i].Lo->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00491 if( lgEOL ) 00492 { 00493 fprintf( ioQQQ, " There should have been a number on this level1 line 4. Sorry.\n" ); 00494 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00495 cdEXIT(EXIT_FAILURE); 00496 } 00497 /*TauLines[i].Lo->g = (realnum)i4;*/ 00498 TauLines[i].Hi->g = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00499 if( lgEOL ) 00500 { 00501 fprintf( ioQQQ, " There should have been a number on this level1 line 5. Sorry.\n" ); 00502 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00503 cdEXIT(EXIT_FAILURE); 00504 } 00505 /*TauLines[i].Hi->g = (realnum)i5;*/ 00506 00507 /* gf is log if negative */ 00508 TauLines[i].Emis->gf = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00509 if( lgEOL ) 00510 { 00511 fprintf( ioQQQ, " There should have been a number on this level1 line 6. Sorry.\n" ); 00512 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00513 cdEXIT(EXIT_FAILURE); 00514 } 00515 /*TauLines[i].Emis->gf = (realnum)f2;*/ 00516 if( TauLines[i].Emis->gf < 0. ) 00517 TauLines[i].Emis->gf = (realnum)pow((realnum)10.f,TauLines[i].Emis->gf); 00518 00519 /* Emis.Aul is log if negative */ 00520 TauLines[i].Emis->Aul = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00521 if( lgEOL ) 00522 { 00523 fprintf( ioQQQ, " There should have been a number on this level1 line 7. Sorry.\n" ); 00524 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00525 cdEXIT(EXIT_FAILURE); 00526 } 00527 if( TauLines[i].Emis->Aul < 0. ) 00528 TauLines[i].Emis->Aul = (realnum)pow((realnum)10.f,TauLines[i].Emis->Aul); 00529 /*TauLines[i].Emis->Aul = f3;*/ 00530 00531 TauLines[i].Emis->iRedisFun = (int)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00532 if( lgEOL ) 00533 { 00534 fprintf( ioQQQ, " There should have been a number on this level1 line 8. Sorry.\n" ); 00535 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 00536 cdEXIT(EXIT_FAILURE); 00537 } 00538 /*TauLines[i].Emis->iRedisFun = i6;*/ 00539 00540 /* this is new in c94.01 and returns nothing (0.) for most lines */ 00541 TauLines[i].Coll.cs = (realnum)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00542 00543 /* finally increment i */ 00544 ++i; 00545 } 00546 } 00547 00548 /* now close the file */ 00549 fclose(ioDATA); 00550 00551 /************************************************************* 00552 * * 00553 * get the level 2 line, opacity project, data set * 00554 * * 00555 *************************************************************/ 00556 00557 /* nWindLine is initialized to the dimension of the vector when it is 00558 * initialized in the definition at the start of this file. 00559 * it is set to -1 with the "no level2" command, which 00560 * stops us from trying to establish this vector */ 00561 if( nWindLine > 0 ) 00562 { 00563 /* begin level2 level 2 wind block */ 00564 /* open file with level 2 line data */ 00565 00566 /* create the TauLine2 emline array */ 00567 TauLine2 = ((transition *)MALLOC( (size_t)nWindLine*sizeof(transition ))); 00568 cs1_flag_lev2 = ((realnum *)MALLOC( (size_t)nWindLine*sizeof(realnum ))); 00569 00570 /* first initialize entire array to dangerously large negative numbers */ 00571 for( i=0; i< nWindLine; ++i ) 00572 { 00573 /* >>chng 99 jul 16, from setting all t[] to flt_max, to call for 00574 * following, each member of structure set to own type of impossible value */ 00575 TransitionJunk( &TauLine2[i] ); 00576 00577 TauLine2[i].Hi = AddState2Stack(); 00578 TauLine2[i].Lo = AddState2Stack(); 00579 TauLine2[i].Emis = AddLine2Stack( true ); 00580 } 00581 00582 if( trace.lgTrace ) 00583 fprintf( ioQQQ," atmdat_readin opening level2.dat:"); 00584 00585 ioDATA = open_data( "level2.dat", "r" ); 00586 00587 /* get magic number off first line */ 00588 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00589 { 00590 fprintf( ioQQQ, " level2.dat error getting magic number\n" ); 00591 cdEXIT(EXIT_FAILURE); 00592 } 00593 00594 /* scan off three numbers, should be the yr mn dy of creation date */ 00595 i = 1; 00596 nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00597 nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00598 ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00599 00600 /* level2.dat magic number 00601 * the following is the set of numbers that appear at the start of level2.dat 01 05 29 00602 * >>chng 06 aug 10, chng to 06, 08, 10 */ 00603 if( lgEOL || ( nelem != 6 ) || ( nelec != 8 ) || ( ion != 10 ) ) 00604 { 00605 fprintf( ioQQQ, 00606 " atmdat_readin: the version of level2.dat is not the current version.\n" ); 00607 fprintf( ioQQQ, 00608 " I expected to find the number 06 08 10 and got %li %li %li instead.\n" , 00609 nelem , nelec , ion ); 00610 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00611 fprintf( ioQQQ, "Please obtain the correct version.\n"); 00612 cdEXIT(EXIT_FAILURE); 00613 } 00614 00615 /* now get the actual data */ 00616 i = 0; 00617 while( i < nWindLine ) 00618 { 00619 /* this must be double for sscanf to work below */ 00620 double tt[7]; 00621 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00622 { 00623 fprintf( ioQQQ, " level2.dat error getting line %li\n", i ); 00624 cdEXIT(EXIT_FAILURE); 00625 } 00626 /* option to skip any line starting with # */ 00627 if( chLine[0]!='#' ) 00628 { 00629 /*printf("string is %s\n",chLine );*/ 00630 sscanf( chLine , "%lf %lf %lf %lf %lf %lf %lf " , 00631 &tt[0] , 00632 &tt[1] , 00633 &tt[2] , 00634 &tt[3] , 00635 &tt[4] , 00636 &tt[5] , 00637 &tt[6] ); 00638 /* these are readjusted into their final form in the structure 00639 * in routine lines_setup*/ 00640 TauLine2[i].Hi->nelem = (int)tt[0]; 00641 TauLine2[i].Hi->IonStg = (int)tt[1]; 00642 TauLine2[i].Lo->g = (realnum)tt[2]; 00643 TauLine2[i].Hi->g = (realnum)tt[3]; 00644 TauLine2[i].Emis->gf = (realnum)tt[4]; 00645 TauLine2[i].EnergyWN = (realnum)tt[5]; 00646 cs1_flag_lev2[i] = (realnum)tt[6]; 00647 ++i; 00648 } 00649 } 00650 /* get magic number off last line */ 00651 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00652 { 00653 fprintf( ioQQQ, " level2.dat error getting last magic number\n" ); 00654 cdEXIT(EXIT_FAILURE); 00655 } 00656 sscanf( chLine , "%ld" , &magic2 ); 00657 if( 999 != magic2 ) 00658 { 00659 fprintf( ioQQQ, " level2.dat ends will wrong magic number=%ld \n", 00660 magic2 ); 00661 cdEXIT(EXIT_FAILURE); 00662 } 00663 fclose( ioDATA ); 00664 if( trace.lgTrace ) 00665 fprintf( ioQQQ," OK \n"); 00666 00667 /* end of level2 block*/ 00668 } 00669 00670 /************************* 00671 * 00672 * get the UTA line sets 00673 * >>chng 06 nov 26, use Gu et al. 2006 UTA data by default 00674 * with option to use older Behar 00675 * 00676 *************************/ 00677 00678 /* these will count number of UTA lines of each type */ 00679 nUTA = 0; 00680 nUTA_Gu06 = 0; 00681 nUTA_Behar = 0; 00682 nUTA_Romas = 0; 00683 /* this is option to use older Behar data rather than new Gu data 00684 * variable lgInnerShell_Gu06 is true by default, set false with 00685 * SET UTA BEHAR command */ 00686 if( ionbal.lgInnerShell_Gu06 ) 00687 { 00688 if( trace.lgTrace ) 00689 fprintf( ioQQQ," atmdat_readin opening UTA_Gu06.dat:"); 00690 00691 ioGU06 = open_data( "UTA_Gu06.dat", "r" ); 00692 00693 /* count how many non-comment lines there are in Gu06 file */ 00694 while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL ) 00695 { 00696 /* we want to count the lines that do not start with # 00697 * since these contain data */ 00698 if( chLine[0] != '#') 00699 ++nUTA_Gu06; 00700 } 00701 /* we counted a magic number as a real line, back off by one */ 00702 ASSERT( nUTA_Gu06 > 0 ); 00703 --nUTA_Gu06; 00704 00705 /* now rewind the file so we can read it a second time*/ 00706 if( fseek( ioGU06 , 0 , SEEK_SET ) != 0 ) 00707 { 00708 fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Gu06.dat.\n"); 00709 cdEXIT(EXIT_FAILURE); 00710 } 00711 00712 /* get up to magic number in Gu 06 file - there is a large header of 00713 * comments with the first data the magic number */ 00714 /* count how many non-comment lines there are in Gu06 file */ 00715 while( read_whole_line( chLine , (int)sizeof(chLine) , ioGU06 ) != NULL ) 00716 { 00717 /* we want to break on first line that does not start with # 00718 * since that contains data */ 00719 if( chLine[0] != '#') 00720 break; 00721 } 00722 /* now get Gu magic number */ 00723 j = 1; 00724 nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00725 nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00726 ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00727 00728 /* is magic number correct? 00729 * the following is the set of numbers that appear at the start of UTA_Gu06.dat 2006 11 17 */ 00730 if( ( nelem != 2007 ) || ( nelec != 1 ) || ( ion != 23 ) ) 00731 { 00732 fprintf( ioQQQ, 00733 " atmdat_readin: the version of UTA_Gu06.dat is not the current version.\n" ); 00734 fprintf( ioQQQ, 00735 " I expected to find the number 2007 1 23 and got %li %li %li instead.\n" , 00736 nelem , nelec , ion ); 00737 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00738 cdEXIT(EXIT_FAILURE); 00739 } 00740 } 00741 else 00742 { 00743 /* this branch, get the Behar 01 data */ 00744 /* >>refer Fe inner shell UTA Behar, E., Sako, M, Kahn, S.M., 00745 * >>refercon 2001, ApJ, 563, 497-504 00746 * the fits were based on the paper 00747 * >>refer Fe inner shell UTA Behar, E., & Netzer, H., 2002, ApJ, 570, 165-170 */ 00748 00749 if( trace.lgTrace ) 00750 fprintf( ioQQQ," atmdat_readin opening UTA_Behar.dat:"); 00751 00752 ioBEHAR = open_data( "UTA_Behar.dat", "r" ); 00753 00754 /* count number of lines in Behar file 00755 * first line is a version number and does not count */ 00756 if( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL ) 00757 { 00758 fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n"); 00759 cdEXIT(EXIT_FAILURE); 00760 } 00761 00762 j = 1; 00763 /* UTA_Behar.dat magic number */ 00764 nelem = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00765 nelec = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00766 ion = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00767 00768 /* magic number 00769 * the following is the set of numbers that appear at the start of UTA_Behar.dat 2002 8 19 */ 00770 /* >>chng 02 oct 22, had incorrect stat weights for all stages of ionization except atom. 00771 * no longer used stored values when stat weight is zero */ 00772 if( ( nelem != 2002 ) || ( nelec != 10 ) || ( ion != 22 ) ) 00773 { 00774 fprintf( ioQQQ, 00775 " atmdat_readin: the version of UTA_Behar.dat is not the current version.\n" ); 00776 fprintf( ioQQQ, 00777 " I expected to find the number 2002 10 22 and got %li %li %li instead.\n" , 00778 nelem , nelec , ion ); 00779 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00780 cdEXIT(EXIT_FAILURE); 00781 } 00782 } 00783 00784 if( trace.lgTrace ) 00785 fprintf( ioQQQ," atmdat_readin UTA data files open OK\n"); 00786 00787 /* count how many lines are in the file, ignoring all lines 00788 * starting with '#' */ 00789 nUTA_Behar = 0; 00790 00791 /* first fill in Behar & Netzer interpolation - first pass just count number of lines */ 00792 for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO ) 00793 { 00794 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00795 { 00796 ++nUTA_Behar; 00797 } 00798 } 00799 00800 if( !ionbal.lgInnerShell_Gu06 ) 00801 { 00802 /* count number of lines in Behar file if we are going to use it 00803 * wee always use above inner shell UTA from B&N */ 00804 while( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) != NULL ) 00805 { 00806 /* we want to count the lines that do not start with # 00807 * since these contain data */ 00808 if( chLine[0] != '#') 00809 ++nUTA_Behar; 00810 } 00811 /* now rewind the file so we can read it a second time*/ 00812 if( fseek( ioBEHAR , 0 , SEEK_SET ) != 0 ) 00813 { 00814 fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Behar.dat.\n"); 00815 cdEXIT(EXIT_FAILURE); 00816 } 00817 /* first line is a version number and does not count */ 00818 if( read_whole_line( chLine , (int)sizeof(chLine) , ioBEHAR ) == NULL ) 00819 { 00820 fprintf( ioQQQ, " atmdat_readin could not read first line of UTA_Behar.dat.\n"); 00821 cdEXIT(EXIT_FAILURE); 00822 } 00823 } 00824 00825 /* now read Romas Kisielius data, taken from 00826 * >>refer Fe UTA Kisielius, R., Hibbert, A., Ferland, G.J., Foord, M.E., Rose, S.J., 00827 * >>refercon van Hoof, P.A.M. & Keenan, F.P. 2003, MNRAS, 344, 696 */ 00828 00829 if( trace.lgTrace ) 00830 fprintf( ioQQQ," atmdat_readin opening UTA_Kisielius.dat:"); 00831 00832 ioROMAS = open_data( "UTA_Kisielius.dat", "r" ); 00833 00834 nUTA_Romas = 0; 00835 while( read_whole_line( chLine , (int)sizeof(chLine) , ioROMAS ) != NULL ) 00836 { 00837 /* we want to count the lines that do not start with # 00838 * since these contain data */ 00839 if( chLine[0] != '#') 00840 { 00841 ++nUTA_Romas; 00842 } 00843 } 00844 00845 /* now malloc the UTA lines array with enough space for both 00846 * data sets */ 00847 UTALines = (transition *)MALLOC( (size_t)(nUTA_Behar+nUTA_Romas+nUTA_Gu06)*sizeof(transition) ); 00848 00849 for( i=0; i<nUTA_Behar+nUTA_Romas+nUTA_Gu06; i++ ) 00850 { 00851 TransitionJunk( &UTALines[i] ); 00852 UTALines[i].Hi = AddState2Stack(); 00853 UTALines[i].Lo = AddState2Stack(); 00854 UTALines[i].Emis = AddLine2Stack( true ); 00855 } 00856 00857 /* these will be used to keep track of what is in each data set */ 00858 for( ion=0; ion<ipIRON; ++ion ) 00859 { 00860 nFeIonRomas[ion] = 0; 00861 WLloRomas[ion] = 1e10; 00862 WLhiRomas[ion] = 0.; 00863 00864 nFeIonBehar[ion] = 0; 00865 WLloBehar[ion] = 1e10; 00866 WLhiBehar[ion] = 0.; 00867 00868 nFeIonGu[ion] = 0; 00869 WLloGu[ion] = 1e10; 00870 WLhiGu[ion] = 0.; 00871 } 00872 00873 /* this will count number of lines, must equal sum of UTAs at end */ 00874 i = 0; 00875 /* now rewind the file so we can read it a second time*/ 00876 if( fseek( ioROMAS , 0 , SEEK_SET ) != 0 ) 00877 { 00878 fprintf( ioQQQ, " atmdat_readin could not rewind UTA_Kisielius.dat.\n"); 00879 cdEXIT(EXIT_FAILURE); 00880 } 00881 00882 i = 0; 00883 /* first fill in Behar & Netzer interpolation - now do real thing, 00884 * these are Nick Abel's fits to the BN data */ 00885 for( ipISO=ipLI_LIKE; ipISO<=ipF_LIKE; ++ipISO ) 00886 { 00887 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00888 { 00889 double f1; 00890 double sum_spon_auto; 00891 # define N 7 00892 00893 double alam[N]={ 1.859035 , 4.692138 , 10.324543 , 19.294364 , 00894 40.401292 , 44.442754 , 72.959719 }; 00895 double blam[N]={-9.7855968,-11.159739, -12.914931 , -14.987272, 00896 -18.853446 , -19.271781 ,-23.828388 }; 00897 double clam[N]={ 8.2874628, 8.3043002, 8.3164038 , 8.3269937, 00898 8.4312895 , 8.3730893 , 8.493802 }; 00899 /* >>chng 05 mar 07, by NA, update fits to better form */ 00900 00901 double aA[N] = {1.9067 , 1.8715 , 2.1033, 2.2815 , 00902 1.9511 , 1.9966 , 2.0852}; 00903 double bA[N] = {8.4538 , 8.5528, 7.616, 6.9247 , 00904 7.9479, 7.9777 , 7.8349 }; 00905 00906 double aDep[N] = {29.54 , 12.07 , 24.23 , 7.35 , 7.37 , 11.14 , 11.99 }; 00907 double bDep[N] = {-14.853 , -7.606 , -7.844 , 9.987 , 10.503 , 8.541 , 8.865 }; 00908 00909 /* NB - these are plus one is because 1 will be subtracted when used */ 00910 /* derive element and ion */ 00911 UTALines[i].Hi->nelem = nelem+1; 00912 ion = nelem + 1 - ipISO; 00913 UTALines[i].Hi->IonStg = ion; 00914 00915 /* these were the statistical weights */ 00917 UTALines[i].Hi->g = 4.f; 00918 UTALines[i].Lo->g = 2.f; 00919 00920 f1 = alam[ipISO-2]*1e-4 + blam[ipISO-2]*1e-4*(nelem+1) + 00921 clam[ipISO-2]*1e-4*POW2(nelem+1.); 00922 f1 = (1./f1); 00923 UTALines[i].WLAng = (realnum)f1; 00924 UTALines[i].EnergyWN = (realnum)(1e8/f1); 00925 UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN; 00926 00927 f1 = aA[ipISO-2]*log((double)(nelem+1)) + bA[ipISO-2]; 00928 00929 UTALines[i].Emis->Aul = (realnum)pow(10., f1 ); 00930 /* generate Emis.gf from A */ 00931 UTALines[i].Emis->gf = 00932 (realnum)(UTALines[i].Emis->Aul*UTALines[i].Hi->g* 00933 1.4992e-8*POW2(1e4/UTALines[i].EnergyWN)); 00934 00935 UTALines[i].Emis->iRedisFun = 1; 00936 00937 /* find sum of spontaneous relax and autoionization As */ 00938 f1 = log((double)(nelem+1)); 00939 f1 = POW2( f1 ); 00940 sum_spon_auto = aDep[ipISO-2]*1e-2*f1 + bDep[ipISO-2]*1e-1; 00941 00942 if( ipISO == ipBE_LIKE ) 00943 { 00944 sum_spon_auto = exp( sum_spon_auto ); 00945 sum_spon_auto = pow(10., sum_spon_auto )*1e13; 00946 } 00947 else if( ipISO <= ipF_LIKE ) 00948 { 00949 sum_spon_auto = pow(10., sum_spon_auto )*1e13; 00950 } 00951 else 00952 TotalInsanity(); 00953 00954 /* last number is negative branching ratio for autoionization, 00955 * which we will store as negative of cs so not misinterpreted as cs */ 00956 /* >>chng 03 feb 18, change to negative of br */ 00957 00958 /*ASSERT( UTALines[i].Emis->Aul / sum_spon_auto <= 1. && UTALines[i].Aul / sum_spon_auto >= 0. );*/ 00959 UTALines[i].Coll.cs = (realnum)MIN2(0., UTALines[i].Emis->Aul / sum_spon_auto - 1. ); 00960 00961 /* option to print UTAs */ 00962 # if 0 00963 fprintf(ioQQQ," %2s%2li\t%.4f\t%.3e\t%.3e\t%.3e\n", 00964 elementnames.chElementSym[nelem], 00965 ion, 00966 UTALines[i].WLAng , 00967 UTALines[i].Emis->Aul , 00968 sum_spon_auto , 00969 -UTALines[i].cs); 00970 # endif 00971 00972 /* finally increment i */ 00973 ++i; 00974 # undef N 00975 } 00976 } 00977 00978 /* next read in the Gu file */ 00979 if( ionbal.lgInnerShell_Gu06 ) 00980 { 00981 ioDATA = ioGU06; 00982 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00983 { 00984 /* only look at lines without '#' in first col */ 00985 if( chLine[0] != '#') 00986 { 00987 long int i2; 00988 double WLAng, Aul, oscill , Aauto; 00989 00990 sscanf( chLine ,"%4li%5li%8lf%13le%12le", 00991 &ion,&i2,&WLAng,&Aul,&Aauto); 00992 sscanf( &chLine[54] ,"%13le", &oscill); 00993 00994 /* all these are iron, first number was ion stage with 0 the atom */ 00995 UTALines[i].Hi->nelem = ipIRON+1; 00996 00997 /* now do stage of ionization */ 00998 ASSERT( ion >= 0 && ion < ipIRON ); 00999 /* the ion stage - 1 is atom - this data file has different 01000 * format from other two - others are 0 for atom */ 01001 UTALines[i].Hi->IonStg = ion; 01002 01003 /* these were the statistical weights 01004 * not included in the original data files, and not needed, 01005 * so use unity */ 01006 UTALines[i].Hi->g = 1.f; 01007 UTALines[i].Lo->g = 1.f; 01008 01009 /* wavelength in Angstroms */ 01010 UTALines[i].WLAng = (realnum)WLAng; 01011 01012 /* keep track of what we have */ 01013 ++nFeIonGu[ion-1]; 01014 WLloGu[ion-1] = MIN2( WLAng , WLloGu[ion-1] ); 01015 WLhiGu[ion-1] = MAX2( WLAng , WLloGu[ion-1] ); 01016 01017 /* energy in wavenumbers */ 01018 UTALines[i].EnergyWN = (realnum)(1e8/WLAng); 01019 UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN; 01020 01021 /* this is true spontaneous rate for doubly excited state to inner 01022 * shell UTA, and is used for pumping, and also relaxing back to inner shell */ 01023 UTALines[i].Emis->Aul = (realnum)Aul; 01024 01025 /* cs is branching ratio for autoionization, 01026 * which we will store as negative cs so not misinterpreted as real cs 01027 * this multiplies the direct pump rate to get ionization rate 01028 * rate evaluated in conv_base.cpp*/ 01029 double frac_ioniz = Aauto/(Aul + Aauto); 01030 ASSERT( frac_ioniz >=0. && frac_ioniz <=1. ); 01031 UTALines[i].Coll.cs = -(realnum)frac_ioniz; 01032 01033 /* save gf scanned from line */ 01034 UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill; 01035 01036 /* statistical weights are not given in the data file - 01037 * assume unity and convert absorption oscillator strength 01038 * into an effective Aul, which is correct for unit statistical 01039 * weight */ 01040 UTALines[i].Emis->Aul = (realnum)eina(UTALines[i].Emis->gf, 01041 UTALines[i].EnergyWN,UTALines[i].Hi->g); 01042 01043 UTALines[i].Emis->iRedisFun = 1; 01044 { 01045 /* DEBUG OTS rates - turn this on to debug line, continuum or both rates */ 01046 /*@-redef@*/ 01047 enum {DEBUG_LOC=false}; 01048 /*@+redef@*/ 01049 if( DEBUG_LOC && UTALines[i].Hi->nelem==ipIRON+1 && UTALines[i].Hi->IonStg==9) 01050 { 01051 fprintf(ioQQQ,"DEBUG Gu UTA\t%li\t%.3f\t%.3e\t%.3e\t%.3e\n", 01052 ion , WLAng,Aul , 01053 eina(UTALines[i].Emis->gf, 01054 UTALines[i].EnergyWN,UTALines[i].Hi->g), 01055 Aauto ); 01056 } 01057 } 01058 01059 /* finally increment i */ 01060 ++i; 01061 } 01062 } 01063 } 01064 else 01065 { 01066 /* ionbal.lgInnerShell_Gu06 is false, do not use this data */ 01067 nUTA_Gu06 = 0; 01068 01069 /* the Behar 01 data */ 01070 /* read in the Behar data */ 01071 ioDATA = ioBEHAR; 01072 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 01073 { 01074 /* only look at lines without '#' in first col */ 01075 if( chLine[0] != '#') 01076 { 01077 long int i1, i2, i3; 01078 double f1, f2, oscill; 01079 double frac_relax; 01080 01081 sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le", 01082 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill ); 01083 01084 /* all these are iron, first number was ion stage with 0 the atom */ 01085 UTALines[i].Hi->nelem = ipIRON+1; 01086 01087 /* now do stage of ionization */ 01088 ASSERT( i1 >= 0 && i1 < ipIRON ); 01089 /* NB - the plus one is because 1 will be subtracted when used, 01090 * in the original data file i1 is 0 for the atom */ 01091 UTALines[i].Hi->IonStg = i1 + 1; 01092 01093 /* these were the statistical weights */ 01094 if( i2>0 && i3>0 ) 01095 { 01096 UTALines[i].Hi->g = (realnum)i2; 01097 UTALines[i].Lo->g = (realnum)i3; 01098 } 01099 else 01100 { 01101 /* nearly all stages of ionization do not include the stat weights, 01102 * so use unity */ 01103 UTALines[i].Hi->g = 1.f; 01104 UTALines[i].Lo->g = 1.f; 01105 } 01106 01107 /* wavelength in Angstroms */ 01108 UTALines[i].WLAng = (realnum)f1; 01109 01110 /* keep track of what we have */ 01111 ++nFeIonBehar[i1]; 01112 WLloBehar[i1] = MIN2( f1 , WLloBehar[i1] ); 01113 WLhiBehar[i1] = MAX2( f1 , WLhiBehar[i1] ); 01114 01115 /* energy in wavenumbers */ 01116 UTALines[i].EnergyWN = (realnum)(1e8/f1); 01117 UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN; 01118 01119 /* this is true spontaneous rate for doubly excited state to inner shell UTA, 01120 * and is used for pumping, and also relaxing back to inner shell */ 01121 UTALines[i].Emis->Aul = (realnum)f2; 01122 01123 /* >>chng 02 oct 22, use absorption oscillator strength in data file 01124 * since don't have stat weights for most lines */ 01125 UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill; 01126 01127 UTALines[i].Emis->iRedisFun = 1; 01128 01129 /* last number is branching ratio for autoionizing, 01130 * which we will store in negative in cs so not misinterpreted as real cs*/ 01131 /* >>chng 03 feb 18, cs to negative of branching ratio */ 01132 ASSERT( frac_relax >=0. && frac_relax <=1. ); 01133 UTALines[i].Coll.cs = (realnum)(frac_relax-1.); 01134 # if 0 01135 fprintf(ioQQQ,"data set %li line %li %2s%2li\t%.4f\t%.3e\t%.3e\n", 01136 nDataSet, 01137 i, 01138 elementnames.chElementSym[UTALines[i].Hi->nelem-1], 01139 UTALines[i].Hi->IonStg, 01140 UTALines[i].WLAng , 01141 UTALines[i].Emis->Aul , 01142 -UTALines[i].cs); 01143 # endif 01144 01145 /* finally increment the counter i */ 01146 ++i; 01147 } 01148 } 01149 } 01150 01151 /* last read in the Romas Kisielius data last since will use it in preference to other 01152 * two sets 01153 *>>refer Fe UTA Kisielius, R.; Hibbert, A.; Ferland, G. J.; Foord, M. E.; Rose, S. J.; 01154 *>>refercon van Hoof, P. A. M.; Keenan, F. P. 2003, MNRAS, 344, 696 01155 */ 01156 ioDATA = ioROMAS; 01157 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 01158 { 01159 /* only look at lines without '#' in first col */ 01160 if( chLine[0] != '#') 01161 { 01162 long int i1, i2, i3; 01163 double f1, f2, oscill; 01164 double frac_relax; 01165 01166 sscanf( chLine ,"%li\t%li\t%li\t%lf\t%le\t%le\t%le", 01167 &i1,&i2,&i3,&f1,&f2,&frac_relax,&oscill ); 01168 01169 /* all these are iron, first number was ion stage with 0 the atom */ 01170 UTALines[i].Hi->nelem = ipIRON+1; 01171 01172 /* now do stage of ionization */ 01173 ASSERT( i1 >= 0 && i1 < ipIRON ); 01174 /* NB - the plus one is because 1 will be subtracted when used, 01175 * in the original data file i1 is 0 for the atom */ 01176 UTALines[i].Hi->IonStg = i1 + 1; 01177 01178 /* these were the statistical weights */ 01179 if( i2>0 && i3>0 ) 01180 { 01181 UTALines[i].Hi->g = (realnum)i2; 01182 UTALines[i].Lo->g = (realnum)i3; 01183 } 01184 else 01185 { 01186 /* nearly all stages of ionization do not include the stat weights, 01187 * so use unity */ 01188 UTALines[i].Hi->g = 1.f; 01189 UTALines[i].Lo->g = 1.f; 01190 } 01191 01192 /*TauLines[i].Hi->IonStg = i2;*/ 01193 UTALines[i].WLAng = (realnum)f1; 01194 01195 /* keep track of what we have */ 01196 ++nFeIonRomas[i1]; 01197 WLloRomas[i1] = MIN2( f1 , WLloRomas[i1] ); 01198 WLhiRomas[i1] = MAX2( f1 , WLhiRomas[i1] ); 01199 01200 UTALines[i].EnergyWN = (realnum)(1e8/f1); 01201 UTALines[i].EnergyErg = (realnum)(ERG1CM)*UTALines[i].EnergyWN; 01202 01203 /* this is true spontaneous rate for doubly excited state to inner shell, 01204 * and is used for pumping, and also relaxing back to inner shell */ 01205 UTALines[i].Emis->Aul = (realnum)f2; 01206 /* generate gf from A */ 01207 /* >>chng 02 oct 22, use absorption oscillator strength in data file 01208 * since don't have stat weights for most lines */ 01209 UTALines[i].Emis->gf = UTALines[i].Lo->g * (realnum)oscill; 01210 01211 UTALines[i].Emis->iRedisFun = 1; 01212 01213 /* last number is branching ratio for autoionizing, 01214 * which we will store in negative in cs so not misinterpreted as real cs*/ 01215 /* >>chng 03 feb 18, cs to negative of br */ 01216 ASSERT( frac_relax >=0. && frac_relax <=1. ); 01217 UTALines[i].Coll.cs = (realnum)(frac_relax-1.); 01218 01219 /* finally increment i */ 01220 ++i; 01221 } 01222 } 01223 01224 /* these have to agree */ 01225 ASSERT( i == (nUTA_Behar+nUTA_Romas+nUTA_Gu06) ); 01226 01227 if( trace.lgTrace ) 01228 { 01229 fprintf(ioQQQ , "summary of UTA data sets read in\nion numb WLlo WLhi\n"); 01230 fprintf(ioQQQ , "Behar 01 total lines %li\n" , nUTA_Behar ); 01231 for( ion=0; ion<ipIRON;++ion ) 01232 { 01233 if( nFeIonBehar[ion] ) 01234 { 01235 fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n", 01236 ion, 01237 nFeIonBehar[ion], 01238 WLloBehar[ion], 01239 WLhiBehar[ion]); 01240 } 01241 } 01242 01243 fprintf(ioQQQ , "Gu 06 total lines %li\n" , nUTA_Gu06 ); 01244 for( ion=0; ion<ipIRON;++ion ) 01245 { 01246 if( nFeIonGu[ion] ) 01247 { 01248 fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n", 01249 ion, 01250 nFeIonGu[ion], 01251 WLloGu[ion], 01252 WLhiGu[ion]); 01253 } 01254 } 01255 01256 fprintf(ioQQQ , "Romas Kisielius 03 total lines %li\n", nUTA_Romas ); 01257 for( ion=0; ion<ipIRON;++ion ) 01258 { 01259 if( nFeIonRomas[ion] ) 01260 { 01261 fprintf(ioQQQ, "%3li %6li %7.3f %7.3f \n", 01262 ion, 01263 nFeIonRomas[ion], 01264 WLloRomas[ion], 01265 WLhiRomas[ion] ); 01266 } 01267 } 01268 } 01269 01270 /* this is default option to use the data Romas Kisielius computed 01271 * to replace any existing data . */ 01272 01273 /* here there are two options - we can add the Romas Kisielius data to the 01274 * first block of UTA data, but zero out any overlapping data, 01275 * or simply ignore the Romas Kisielius data */ 01276 if( ionbal.lgInnerShell_Kisielius ) 01277 { 01278 for( i=0; i<nUTA_Behar+nUTA_Gu06; ++i ) 01279 { 01280 if( UTALines[i].Hi->nelem-1 == ipIRON ) 01281 { 01282 /* nFeIonRomas counts how many of the Romas Kisielius lines 01283 * fall in this stage of ionization - if zero then 01284 * none, so keep Behar data, if positive, then 01285 * there are new data, and zero out Behar data */ 01286 ASSERT( UTALines[i].Hi->IonStg-1< ipIRON ); 01287 if( nFeIonRomas[UTALines[i].Hi->IonStg-1] ) 01288 { 01289 UTALines[i].Emis->Aul = -1.; 01290 } 01291 } 01292 } 01293 /* the total number of UTAs, a global variable */ 01294 nUTA = nUTA_Behar + nUTA_Gu06 + nUTA_Romas; 01295 } 01296 else 01297 { 01298 /* the Kisielius data come last, and nUTA is the global variable 01299 * that counts the number of UTA lines. Setting nUTA to only the 01300 * first two amounts to ignoring the Kisielius data */ 01301 nUTA = nUTA_Behar + nUTA_Gu06; 01302 } 01303 01304 /* now close the files */ 01305 fclose(ioROMAS); 01306 if( ionbal.lgInnerShell_Gu06 ) 01307 fclose(ioGU06); 01308 else 01309 fclose(ioBEHAR); 01310 01311 /* end UTA */ 01312 01313 /* allocate space for the carbon monoxide CO rotation levels */ 01314 ASSERT( nCORotate > 0); 01315 C12O16Rotate = (transition *)MALLOC( (size_t)nCORotate*sizeof(transition) ); 01316 C13O16Rotate = (transition *)MALLOC( (size_t)nCORotate*sizeof(transition) ); 01317 01318 /* this says we cannot change the number of levels again, since space is allocated */ 01319 lgCORotateMalloc = true; 01320 01321 /* next initialize entire array to dangerously large negative numbers */ 01322 for( J=0; J< nCORotate; ++J ) 01323 { 01324 TransitionJunk( &C12O16Rotate[J] ); 01325 C12O16Rotate[J].Hi = AddState2Stack(); 01326 C12O16Rotate[J].Lo = AddState2Stack(); 01327 C12O16Rotate[J].Emis = AddLine2Stack( true ); 01328 01329 TransitionJunk( &C13O16Rotate[J] ); 01330 C13O16Rotate[J].Hi = AddState2Stack(); 01331 C13O16Rotate[J].Lo = AddState2Stack(); 01332 C13O16Rotate[J].Emis = AddLine2Stack( true ); 01333 } 01334 01335 /* read in data for the set of hyperfine structure lines, and allocate 01336 * space for the transition HFLines[nHFLines] structure */ 01337 HyperfineCreate(); 01338 01339 /* this is Humeshkar's work on generic atoms and molecules. */ 01340 //Nemala_Start(); 01341 01342 /* initialize the large block of level 1 real lines, and OP level 2 lines */ 01343 lines_setup(); 01344 01345 /* mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat mewe_gbar.dat ========*/ 01346 /* read in g-bar data taken from 01347 *>>refer all cs-gbar Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985, 01348 *>>refercon A&AS, 62, 197 */ 01349 /* open file with Mewe coefficients */ 01350 01351 if( trace.lgTrace ) 01352 fprintf( ioQQQ," atmdat_readin opening mewe_gbar.dat:"); 01353 01354 ioDATA = open_data( "mewe_gbar.dat", "r" ); 01355 01356 /* get magic number off first line */ 01357 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01358 { 01359 fprintf( ioQQQ, " mewe_gbar.dat error getting magic number\n" ); 01360 cdEXIT(EXIT_FAILURE); 01361 } 01362 /* check the number */ 01363 sscanf( chLine , "%ld" , &magic1 ); 01364 if( magic1 != 9101 ) 01365 { 01366 fprintf( ioQQQ, " mewe_gbar.dat starts with wrong magic number=%ld \n", 01367 magic1 ); 01368 cdEXIT(EXIT_FAILURE); 01369 } 01370 01371 /* now get the actual data, indices are correct for c, in Fort went from 2 to 210 */ 01372 for( i=1; i < 210; i++ ) 01373 { 01374 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01375 { 01376 fprintf( ioQQQ, " mewe_gbar.dat error getting line %li\n", i ); 01377 cdEXIT(EXIT_FAILURE); 01378 } 01379 /*printf("%s\n",chLine);*/ 01380 double help[4]; 01381 sscanf( chLine, "%lf %lf %lf %lf ", &help[0], &help[1], &help[2], &help[3] ); 01382 for( int l=0; l < 4; ++l ) 01383 MeweCoef.g[i][l] = (realnum)help[l]; 01384 } 01385 01386 /* get magic number off last line */ 01387 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 01388 { 01389 fprintf( ioQQQ, " mewe_gbar.dat error getting last magic number\n" ); 01390 cdEXIT(EXIT_FAILURE); 01391 } 01392 01393 sscanf( chLine , "%ld" , &magic2 ); 01394 01395 if( magic1 != magic2 ) 01396 { 01397 fprintf( ioQQQ, " mewe_gbar.dat ends will wrong magic number=%ld \n", 01398 magic2 ); 01399 cdEXIT(EXIT_FAILURE); 01400 } 01401 01402 fclose( ioDATA ); 01403 01404 if( trace.lgTrace ) 01405 fprintf( ioQQQ," OK \n"); 01406 01407 /* This is what remains of the t_yield initialization 01408 * this should not be in the constructor of t_yield 01409 * since it initializes a different struct */ 01410 for( nelem=0; nelem < LIMELM; nelem++ ) 01411 for( ion=0; ion < LIMELM; ion++ ) 01412 Heavy.nsShells[nelem][ion] = LONG_MAX; 01413 01414 /* now read in all auger yields 01415 * will do elements from li on up, 01416 * skip any line starting with # 01417 * this loop goes from lithium to Zn */ 01418 for( nelem=2; nelem < LIMELM; nelem++ ) 01419 { 01420 /* nelem is on the shifted C scale, so 2 is Li */ 01421 for( ion=0; ion <= nelem; ion++ ) 01422 { 01423 /* number of bound electrons, = atomic number for neutral */ 01424 nelec = nelem - ion + 1; 01425 /* one of dima's routines to determine the number of electrons 01426 * for this species, nelem +1 to shift to physical number */ 01427 /* subroutine atmdat_outer_shell(iz,in,imax,ig0,ig1) 01428 * iz - atomic number from 1 to 30 (integer) 01429 * in - number of electrons from 1 to iz (integer) 01430 * Output: imax - number of the outer shell 01431 */ 01432 atmdat_outer_shell(nelem+1,nelec,&imax,&ig0,&ig1); 01433 01434 ASSERT( imax > 0 && imax <= 10 ); 01435 01436 /* nsShells[nelem][ion] is outer shell number for ion with nelec electrons 01437 * on physics scale, with K shell being 1 */ 01438 Heavy.nsShells[nelem][ion] = imax; 01439 } 01440 } 01441 01442 /************************************************************* 01443 * * 01444 * get the Auger electron yield data set * 01445 * * 01446 *************************************************************/ 01447 01448 t_yield::Inst().init_yield(); 01449 01450 /**************************************************************** 01451 * * 01452 * get the Hummer and Storey model case A, B results, these are * 01453 * the two data files e1b.dat and e2b.dat, for H and He * 01454 * * 01455 ****************************************************************/ 01456 01457 /* now get Hummer and Storey case b data, loop over H, He */ 01458 /* >>chng 01 aug 08, generalized this to both case A and B, and all elements 01459 * up to HS_NZ, now 8 */ 01460 for( ipZ=0; ipZ<HS_NZ; ++ipZ ) 01461 { 01462 /* don't do the minor elements, Li-B */ 01463 if( ipZ>1 && ipZ<5 ) continue; 01464 01465 for( iCase=0; iCase<2; ++iCase ) 01466 { 01467 /* open file with case B data */ 01468 /* >>chng 01 aug 08, add HS_ to start of file names to indicate Hummer Storey origin */ 01469 /* create file name for this charge 01470 * first character after e is charge number for this element, 01471 * then follows whether this is case A or case B data */ 01472 sprintf( chFilename, "HS_e%ld%c.dat", ipZ+1, ( iCase == 0 ) ? 'a' : 'b' ); 01473 01474 if( trace.lgTrace ) 01475 fprintf( ioQQQ," atmdat_readin opening Hummer Storey emission file %s:",chFilename); 01476 01477 /* open the file */ 01478 ioDATA = open_data( chFilename, "r" ); 01479 01480 if( trace.lgTrace ) 01481 { 01482 fprintf( ioQQQ," OK\n"); 01483 } 01484 01485 /* read in the number of temperatures and densities*/ 01486 i = fscanf( ioDATA, "%li %li ", 01487 &atmdat.ntemp[iCase][ipZ], &atmdat.nDensity[iCase][ipZ] ); 01488 01489 /* check that ntemp and nDensity are below NHSDIM, 01490 * set to 15 in atmdat_HS_caseB.h */ 01491 assert (atmdat.ntemp[iCase][ipZ] >0 && atmdat.ntemp[iCase][ipZ] <= NHSDIM ); 01492 assert (atmdat.nDensity[iCase][ipZ] > 0 && atmdat.nDensity[iCase][ipZ] <= NHSDIM); 01493 01494 /* loop reading in line emissivities for all temperatures*/ 01495 for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][ipZ]; ipTemp++ ) 01496 { 01497 for( ipDens=0; ipDens < atmdat.nDensity[iCase][ipZ]; ipDens++ ) 01498 { 01499 long int junk, junk2 , ne; 01500 fscanf( ioDATA, " %lf %li %lf %c %li %ld ", 01501 &atmdat.Density[iCase][ipZ][ipDens], &junk , 01502 &atmdat.ElecTemp[iCase][ipZ][ipTemp], &cha , &junk2 , 01503 &atmdat.ncut[iCase][ipZ] ); 01504 01505 ne = atmdat.ncut[iCase][ipZ]*(atmdat.ncut[iCase][ipZ] - 1)/2; 01506 ASSERT( ne<=NLINEHS ); 01507 for( j=0; j < ne; j++ ) 01508 { 01509 fscanf( ioDATA, "%lf ", &atmdat.Emiss[iCase][ipZ][ipTemp][ipDens][j] ); 01510 } 01511 } 01512 } 01513 01514 /*this is end of read-in loop */ 01515 fclose(ioDATA); 01516 # if 0 01517 /* print list of densities and temperatures */ 01518 for( ipDens=0; ipDens<atmdat.nDensity[iCase][ipZ]; ipDens++ ) 01519 { 01520 fprintf(ioQQQ," %e,", atmdat.Density[iCase][ipZ][ipDens]); 01521 } 01522 fprintf(ioQQQ,"\n"); 01523 for( ipTemp=0; ipTemp<atmdat.ntemp[iCase][ipZ]; ipTemp++ ) 01524 { 01525 fprintf(ioQQQ," %e,", atmdat.ElecTemp[iCase][ipZ][ipTemp]); 01526 } 01527 fprintf(ioQQQ,"\n"); 01528 # endif 01529 } 01530 } 01531 01532 /* Verner's model FeII atom 01533 * do not call if Netzer model used, or it Fe(2) is zero 01534 * exception is when code is searching for first soln */ 01535 /* read the atomic data from files */ 01536 /* convert line arrays to internal form needed by code */ 01537 FeIICreate(); 01538 return; 01539 } 01540 01541 t_yield::t_yield() 01542 { 01543 /* preset two arrays that will hold auger data 01544 * set to very large values so that code will blow 01545 * if not set properly below */ 01546 for( int nelem=0; nelem < LIMELM; nelem++ ) 01547 { 01548 for( int ion=0; ion < LIMELM; ion++ ) 01549 { 01550 for( int ns=0; ns < 7; ns++ ) 01551 { 01552 n_elec_eject[nelem][ion][ns] = LONG_MAX; 01553 for( int nelec=0; nelec < 10; nelec++ ) 01554 { 01555 frac_elec_eject[nelem][ion][ns][nelec] = FLT_MAX; 01556 } 01557 } 01558 } 01559 } 01560 01561 lgKillAuger = false; 01562 } 01563 01564 void t_yield::init_yield() 01565 { 01566 char chLine[FILENAME_PATH_LENGTH_2]; 01567 const char* chFilename; 01568 01569 /* following is double for sscanf to work */ 01570 double temp[15]; 01571 01572 DEBUG_ENTRY( "init_yield()" ); 01573 01574 /************************************************************* 01575 * * 01576 * get the Auger electron yield data set * 01577 * * 01578 *************************************************************/ 01579 01580 /* NB NB -- This test of Heavy.nsShells remains here to assure 01581 * that it contains meaningful values since needed below, once 01582 * t_Heavy is a Singleton, it can be removed !!! */ 01583 ASSERT( Heavy.nsShells[2][0] > 0 ); 01584 01585 /* hydrogen and helium will not be done below, so set yields here*/ 01586 n_elec_eject[ipHYDROGEN][0][0] = 1; 01587 n_elec_eject[ipHELIUM][0][0] = 1; 01588 n_elec_eject[ipHELIUM][1][0] = 1; 01589 01590 frac_elec_eject[ipHYDROGEN][0][0][0] = 1; 01591 frac_elec_eject[ipHELIUM][0][0][0] = 1; 01592 frac_elec_eject[ipHELIUM][1][0][0] = 1; 01593 01594 /* open file auger.dat, yield data file that came from 01595 * >>refer all auger Kaastra, J.S., & Mewe, R. 1993, A&AS, 97, 443-482 */ 01596 chFilename = "mewe_nelectron.dat"; 01597 01598 if( trace.lgTrace ) 01599 fprintf( ioQQQ, " init_yield opening %s:", chFilename ); 01600 01601 FILE *ioDATA; 01602 01603 /* open the file */ 01604 ioDATA = open_data( chFilename, "r" ); 01605 01606 /* now read in all auger yields 01607 * will do elements from li on up, 01608 * skip any line starting with # 01609 * this loop goes from lithium to Zn */ 01610 for( int nelem=2; nelem < LIMELM; nelem++ ) 01611 { 01612 /* nelem is on the shifted C scale, so 2 is Li */ 01613 for( int ion=0; ion <= nelem; ion++ ) 01614 { 01615 for( int ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 01616 { 01617 char ch1 = '#'; 01618 /* the * is the old comment char, accept it, but really want # */ 01619 while( ch1 == '#' || ch1 == '*' ) 01620 { 01621 if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL ) 01622 { 01623 fprintf( ioQQQ, " %s error getting line %i\n", chFilename, ns ); 01624 cdEXIT(EXIT_FAILURE); 01625 } 01626 ch1 = chLine[0]; 01627 } 01628 /*printf("string is %s\n",chLine );*/ 01629 sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", 01630 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4], 01631 &temp[5], &temp[6], &temp[7], &temp[8], &temp[9], 01632 &temp[10],&temp[11],&temp[12],&temp[13],&temp[14] ); 01633 n_elec_eject[nelem][ion][ns] = (long int)temp[3]; 01634 01635 ASSERT( n_elec_eject[nelem][ion][ns] >= 0 && n_elec_eject[nelem][ion][ns] < 11 ); 01636 01637 /* can pop off up to 10 auger electrons, these are the probabilities*/ 01638 for( int j=0; j < 10; j++ ) 01639 { 01640 frac_elec_eject[nelem][ion][ns][j] = (realnum)temp[j+4]; 01641 ASSERT( frac_elec_eject[nelem][ion][ns][j] >= 0. ); 01642 } 01643 } 01644 } 01645 /* activate this print statement to get yield for k-shell of all atoms */ 01646 /*fprintf(ioQQQ,"yyyield\t%li\t%.4e\n", nelem+1, frac_elec_eject[nelem][0][0][0] );*/ 01647 } 01648 01649 fclose( ioDATA ); 01650 01651 if( trace.lgTrace ) 01652 { 01653 /* this is set with no auger command */ 01654 if( lgKillAuger ) 01655 fprintf( ioQQQ, " Auger yields will be killed.\n"); 01656 fprintf( ioQQQ, " OK\n"); 01657 } 01658 01659 /* open file mewe_fluor.dat, yield data file that came from 01660 * >>refer all auger Kaastra, J.S., & Mewe, R. 1993, 97, 443-482 */ 01661 chFilename = "mewe_fluor.dat"; 01662 01663 if( trace.lgTrace ) 01664 fprintf( ioQQQ, " init_yield opening %s:", chFilename ); 01665 01666 /* open the file */ 01667 ioDATA = open_data( chFilename, "r" ); 01668 01669 /* now read in all auger yields 01670 * will do elements from li on up, 01671 * skip any line starting with # */ 01672 do 01673 { 01674 if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL ) 01675 { 01676 fprintf( ioQQQ, " %s error getting line %i\n", chFilename, 0 ); 01677 cdEXIT(EXIT_FAILURE); 01678 } 01679 } 01680 while( chLine[0] == '#' ); 01681 01682 bool lgEOL = false; 01683 01684 nfl_lines = 0; 01685 do 01686 { 01687 const int NKM = 10; 01688 int nDima[NKM] = { 0, 1, 2, 2, 3, 4, 4, 5, 5, 6 }; 01689 int nAuger; 01690 01691 if( nfl_lines >= MEWE_FLUOR ) 01692 TotalInsanity(); 01693 01694 /*printf("string is %s\n",chLine );*/ 01695 sscanf( chLine, "%lf %lf %lf %lf %lf %lf %lf", 01696 &temp[0], &temp[1], &temp[2], &temp[3], &temp[4], 01697 &temp[5], &temp[6] ); 01698 01699 /* the atomic number, C is 5 */ 01700 nfl_nelem[nfl_lines] = (int)temp[0]-1; 01701 ASSERT( nfl_nelem[nfl_lines] >= 0 && nfl_nelem[nfl_lines] < LIMELM ); 01702 01703 /* the ion stage for target, atom is 0 */ 01704 nfl_ion[nfl_lines] = (int)temp[1]-1; 01705 ASSERT( nfl_ion[nfl_lines] >= 0 && nfl_ion[nfl_lines] <= nfl_nelem[nfl_lines]+1 ); 01706 01707 /* the target's shell */ 01708 nfl_nshell[nfl_lines] = nDima[(long)temp[2]-1]; 01709 ASSERT( nfl_nshell[nfl_lines] >= 0 && 01710 /* nsShells is shell number, where K is 1 */ 01711 nfl_nshell[nfl_lines] < Heavy.nsShells[nfl_nelem[nfl_lines]][nfl_ion[nfl_lines]]-1 ); 01712 01713 /* this is the number of Auger electrons ejected */ 01714 nAuger = (int)temp[3]; 01715 /* so this is the spectrum of the photons that are emitted */ 01716 nfl_ion_emit[nfl_lines] = nfl_ion[nfl_lines] + nAuger + 1; 01717 /* must be gt 0 since at least photoelectron comes off */ 01718 ASSERT( nfl_ion_emit[nfl_lines] > 0 && nfl_ion_emit[nfl_lines] <= nfl_nelem[nfl_lines]+1); 01719 01720 /* this is the type of line as defined in their paper */ 01721 nfl_nLine[nfl_lines] = (int)temp[4]; 01722 01723 /* energy in Ryd */ 01724 fl_energy[nfl_lines] = (realnum)temp[5] / (realnum)EVRYD; 01725 ASSERT( fl_energy[nfl_lines] > 0. ); 01726 01727 /* fluor yield */ 01728 fl_yield[nfl_lines] = (realnum)temp[6]; 01729 /* NB cannot assert <=1 since data file has yields around 1.3 - 1.4 */ 01730 ASSERT( fl_yield[nfl_lines] >= 0 ); 01731 01732 ++nfl_lines; 01733 01734 do 01735 { 01736 if( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) == NULL ) 01737 lgEOL = true; 01738 } 01739 while( chLine[0]=='#' && !lgEOL ); 01740 } 01741 while( !lgEOL ); 01742 01743 fclose( ioDATA ); 01744 }