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 /*MeanInc increment mean ionization fractions and temperatures over computed structure, in radius_increment */ 00004 /*MeanZero zero mean of ionization fractions array */ 00005 /*MeanIonRadius do radius mean ionization fractions or temperature over radius for any element */ 00006 /*MeanIonVolume do volume mean ionization fractions or temperature over volume for any element */ 00007 /*aver compute average of various quantities over the computed geometry 00008 * called by startenditer to initialize, radinc to increment, and prtfinal for final results */ 00009 #include "cddefines.h" 00010 #include "physconst.h" 00011 #include "radius.h" 00012 #include "dense.h" 00013 #include "hyperfine.h" 00014 #include "magnetic.h" 00015 #include "hmi.h" 00016 #include "phycon.h" 00017 #include "geometry.h" 00018 #include "mean.h" 00019 00020 /*MeanInc increment mean ionization fractions and temperatures over computed structure, in radius_increment */ 00021 void MeanInc(void) 00022 { 00023 long int ion, 00024 nelem; 00025 double dc, 00026 dcv; 00027 00028 /* this routine is called by radius_increment during the calculation to update the 00029 * sums that will become the vol and rad weighted mean ionizations */ 00030 00031 DEBUG_ENTRY( "MeanInc()" ); 00032 00033 for( nelem=0; nelem < LIMELM; nelem++ ) 00034 { 00035 /* this is mean ionization */ 00036 if( nelem==ipHYDROGEN ) 00037 { 00038 /* >>chng 04 jun 27, add this option, previously had not included 00039 * H2 as part of total */ 00040 ion = 2; 00041 /* hydrogen is special case since must include H2, 00042 * and must mult by 2 to conserve total H - will need to div 00043 * by two if column density ever used */ 00044 /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/ 00045 dc = 2.*hmi.H2_total*radius.drad_x_fillfac; 00046 mean.xIonMeans[0][nelem][ion] += dc; 00047 /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/ 00048 dc = 2.*hmi.H2_total*radius.drad_x_fillfac*dense.eden; 00049 mean.xIonEdenMeans[0][nelem][ion] += dc; 00050 /* xIonMeans[1][n][ion] is vol integration */ 00051 dcv = 2.*hmi.H2_total*radius.dVeff; 00052 mean.xIonMeans[1][nelem][ion] += dcv; 00053 /* xIonMeans[1][ion][n] is vol integration PLUS electron density */ 00054 dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden; 00055 mean.xIonEdenMeans[1][nelem][ion] += dcv; 00056 } 00057 for( ion=0; ion < (nelem + 2); ion++ ) 00058 { 00059 /* xIonMeans[0][ion][n] is radial integration */ 00060 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac; 00061 mean.xIonMeans[0][nelem][ion] += dc; 00062 00063 /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/ 00064 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden; 00065 mean.xIonEdenMeans[0][nelem][ion] += dc; 00066 00067 /* xIonMeans[1][n][ion] is vol integration */ 00068 /* >>chng 00 mar 28, r1r0sq had been dVeff, 00069 * bug discovered by Valentina Luridiana */ 00070 dcv = dense.xIonDense[nelem][ion]*radius.dVeff; 00071 mean.xIonMeans[1][nelem][ion] += dcv; 00072 00073 /* xIonMeans[1][ion][n] is vol integration PLUS electron density */ 00074 /* >>chng 00 mar 28, r1r0sq had been dVeff, 00075 * bug discovered by Valentina Luridiana */ 00076 dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden; 00077 mean.xIonEdenMeans[1][nelem][ion] += dcv; 00078 } 00079 /* these normalize the above, use gas_phase which includes molecular parts */ 00080 /* first two are means over radius */ 00081 dc = dense.gas_phase[nelem]*radius.drad_x_fillfac; 00082 mean.xIonMeansNorm[0][nelem] += dc; 00083 dc = dense.gas_phase[nelem]*radius.drad_x_fillfac*dense.eden; 00084 mean.xIonEdenMeansNorm[0][nelem] += dc; 00085 /* next two means over volume */ 00086 /* >>chng 04 jun 28, changed radius.dVeff to dVeff, 00087 * which is equivalent in thin zone limit, more accurate with thick zones, as per 00088 * Valentina Luridiana email */ 00089 /*dcv = dense.gas_phase[nelem]*radius.dVeff;*/ 00090 dcv = dense.gas_phase[nelem]*radius.dVeff; 00091 mean.xIonMeansNorm[1][nelem] += dcv; 00092 /*dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden;*/ 00093 dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden; 00094 mean.xIonEdenMeansNorm[1][nelem] += dcv; 00095 00096 /* this is mean temperature */ 00097 /* this is ionization on the IonFracs scale, offset 1 up since 00098 * zero is total abundance. This works well since the mean 00099 * ionization array xIonMeans is also offset up one, since 00100 * [0] is the normalization */ 00101 if( nelem==ipHYDROGEN ) 00102 { 00103 ion = 2; 00104 /* TempMeans[0][ion][n] is radial integration */ 00105 dc = 2.*hmi.H2_total*radius.drad_x_fillfac; 00106 mean.TempMeans[0][nelem][ion] += dc*phycon.te; 00107 mean.TempMeansNorm[0][nelem][ion] += dc; 00108 00109 /* TempMeans[0][ion][n] is radial integration PLUS electron density*/ 00110 dc = 2.*hmi.H2_total*radius.drad_x_fillfac*dense.eden; 00111 mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te; 00112 mean.TempEdenMeansNorm[0][nelem][ion] += dc; 00113 00114 /* TempMeans[1][ion+1][n] is vol integration */ 00115 /*>>chng 00 dec 18, following had dVeff rather than r1r0sq */ 00116 dcv = 2.*hmi.H2_total*radius.dVeff; 00117 mean.TempMeans[1][nelem][ion] += dcv*phycon.te; 00118 mean.TempMeansNorm[1][nelem][ion] += dcv; 00119 00120 /* TempMeans[1][ion+1][n] is vol integration PLUS electron density */ 00121 dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden; 00122 mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te; 00123 mean.TempEdenMeansNorm[1][nelem][ion] += dcv; 00124 } 00125 for( ion=0; ion < (nelem + 2); ion++ ) 00126 { 00127 /* TempMeans[0][ion][n] is radial integration */ 00128 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac; 00129 mean.TempMeans[0][nelem][ion] += dc*phycon.te; 00130 mean.TempMeansNorm[0][nelem][ion] += dc; 00131 00132 /* TempMeans[0][ion][n] is radial integration PLUS electron density*/ 00133 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden; 00134 mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te; 00135 mean.TempEdenMeansNorm[0][nelem][ion] += dc; 00136 00137 /* TempMeans[1][ion+1][n] is vol integration */ 00138 /*>>chng 00 dec 18, following had dVeff rather than r1r0sq */ 00139 dcv = dense.xIonDense[nelem][ion]*radius.dVeff; 00140 mean.TempMeans[1][nelem][ion] += dcv*phycon.te; 00141 mean.TempMeansNorm[1][nelem][ion] += dcv; 00142 00143 /* TempMeans[1][ion+1][n] is vol integration PLUS electron density */ 00144 dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden; 00145 mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te; 00146 mean.TempEdenMeansNorm[1][nelem][ion] += dcv; 00147 } 00148 } 00149 00150 /* ============================================================= 00151 * 00152 * these are some special quanties for the mean 00153 */ 00154 00155 /* >>chng 05 dec 28, add this 00156 * used to get magnetic field weighted wrt harmonic mean spin temperature, 00157 * * H0 over radius - as measured by 21cm temperature */ 00158 if( hyperfine.Tspin21cm > SMALLFLOAT ) 00159 { 00160 dc = dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac/phycon.te; 00161 } 00162 else 00163 { 00164 dc = 0.; 00165 } 00166 /* mean magnetic field weighted wrt 21 cm opacity, n(H0)/T */ 00167 mean.B_HarMeanTempRadius[0] += sqrt(fabs(magnetic.pressure)*PI8) * dc; 00168 /* this assumes that field is tangled */ 00169 mean.B_HarMeanTempRadius[1] += dc; 00170 00171 /* used to get harmonic mean temperature with respect to H over radius, 00172 * for comparison with 21cm temperature */ 00173 dc = dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac; 00174 00175 /* harmonic mean gas kinetic temp, over radius, for comparison with 21 cm obs */ 00176 /*HEADS UP - next four are the inverse of equation 3 of 00177 * >>refer Tspin 21 cm Abel, N.P., Brogan, C.L., Ferland, G.J., O'Dell, C.R., 00178 * >>refercon Shaw, G. & Troland, T.H. 2004, ApJ, 609, 247 */ 00179 mean.HarMeanTempRadius[0] += dc; 00180 mean.HarMeanTempRadius[1] += dc/phycon.te; 00181 00182 /* harmonic mean gas kinetic temp, over volume, for symmetry with above, for comp with 21 cm obs */ 00183 mean.HarMeanTempVolume[0] += dc*radius.r1r0sq; 00184 mean.HarMeanTempVolume[1] += dc/phycon.te*radius.r1r0sq; 00185 00186 /* harmonic mean of computed 21 cm spin temperature - this is what 21 cm actually measures */ 00187 mean.H_21cm_spin_mean_radius[0] += dc; 00188 mean.H_21cm_spin_mean_radius[1] += dc / SDIV( hyperfine.Tspin21cm ); 00189 00190 /* mean H2 temperature over radius */ 00191 dc = hmi.H2_total*radius.drad_x_fillfac; 00192 mean.H2MeanTempRadius[0] += dc*phycon.te; 00193 mean.H2MeanTempRadius[1] += dc; 00194 00195 /* mean H2 temperature over volume, for symmetry */ 00196 dcv = hmi.H2_total*radius.dVeff; 00197 mean.H2MeanTempVolume[0] += dc*phycon.te; 00198 mean.H2MeanTempVolume[1] += dc; 00199 00200 /* >>chng 05 dec 28, add mean temp over radius and vol */ 00201 dc = radius.drad_x_fillfac; 00202 mean.TempMeanRadius[0] += dc*phycon.te; 00203 mean.TempMeanRadius[1] += dc; 00204 00205 dcv = radius.dVeff; 00206 mean.TempMeanVolume[0] += dc*phycon.te; 00207 mean.TempMeanVolume[1] += dc; 00208 return; 00209 } 00210 00211 /*MeanZero zero mean of ionization fractions array */ 00212 void MeanZero(void) 00213 { 00214 long int ion, 00215 j, 00216 nelem, 00217 limit; 00218 static bool lgFirst=true; 00219 00220 DEBUG_ENTRY( "MeanZero()" ); 00221 00222 /* 00223 * MeanZero is called at start of calculation by zero, and by 00224 * startenditer to initialize the variables 00225 */ 00226 00227 if( lgFirst ) 00228 { 00229 lgFirst = false; 00230 /* both are [2], [nelem][ion] */ 00231 mean.xIonMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ); 00232 mean.xIonEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ); 00233 mean.xIonMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) ); 00234 mean.xIonEdenMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) ); 00235 mean.TempMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ); 00236 mean.TempEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ); 00237 mean.TempMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ); 00238 mean.TempEdenMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ); 00239 for( j=0; j<2; ++j ) 00240 { 00241 mean.xIonMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ); 00242 mean.xIonEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ); 00243 mean.xIonMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) ); 00244 mean.xIonEdenMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) ); 00245 mean.TempMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ); 00246 mean.TempEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ); 00247 mean.TempMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ); 00248 mean.TempEdenMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ); 00249 for( nelem=0; nelem<LIMELM; ++nelem ) 00250 { 00251 limit = MAX2(3,nelem+2); 00252 mean.xIonMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ); 00253 mean.xIonEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ); 00254 mean.TempMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ); 00255 mean.TempEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ); 00256 mean.TempMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ); 00257 mean.TempEdenMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ); 00258 } 00259 } 00260 } 00261 00262 for( j=0; j < 2; j++ ) 00263 { 00264 for( nelem=0; nelem < LIMELM; nelem++ ) 00265 { 00266 mean.xIonMeansNorm[j][nelem] = 0.; 00267 mean.xIonEdenMeansNorm[j][nelem] = 0.; 00268 /* >>chng 04 jun 27, H2 will be 3rd ion stage of H */ 00269 limit = MAX2(3,nelem+2); 00270 /*for( ion=0; ion <= (nelem + 2); ion++ )*/ 00271 for( ion=0; ion < limit; ion++ ) 00272 { 00273 mean.TempMeansNorm[j][nelem][ion] = 0.; 00274 mean.TempEdenMeansNorm[j][nelem][ion] = 0.; 00275 /* these are used to save info on temperature and ionization means */ 00276 mean.xIonMeans[j][nelem][ion] = 0.; 00277 mean.xIonEdenMeans[j][nelem][ion] = 0.; 00278 /* double here since [2] and [3] have norm for average */ 00279 mean.TempMeans[j][nelem][ion] = 0.; 00280 mean.TempEdenMeans[j][nelem][ion] = 0.; 00281 } 00282 } 00283 } 00284 00285 /* mean over radius, for comp with 21 cm obs */ 00286 mean.HarMeanTempRadius[0] = 0.; 00287 mean.HarMeanTempRadius[1] = 0.; 00288 00289 /* mean over volume, for symmetry */ 00290 mean.HarMeanTempVolume[0] = 0.; 00291 mean.HarMeanTempVolume[1] = 0.; 00292 00293 /* mena of calculated 21 cm spin temperature */ 00294 mean.H_21cm_spin_mean_radius[0] = 0.; 00295 mean.H_21cm_spin_mean_radius[1] = 0.; 00296 00297 /* H2 mean temp over radius */ 00298 mean.H2MeanTempRadius[0] = 0.; 00299 mean.H2MeanTempRadius[1] = 0.; 00300 /* H2 mean temp over volume */ 00301 mean.H2MeanTempVolume[0] = 0.; 00302 mean.H2MeanTempVolume[1] = 0.; 00303 00304 mean.TempMeanRadius[0] = 0.; 00305 mean.TempMeanRadius[1] = 0.; 00306 mean.TempMeanVolume[0] = 0.; 00307 mean.TempMeanVolume[1] = 0.; 00308 00309 mean.B_HarMeanTempRadius[0] = 0.; 00310 mean.B_HarMeanTempRadius[1] = 0.; 00311 return; 00312 } 00313 00314 /*MeanIonRadius derive mean ionization fractions over radius for some element */ 00315 void MeanIonRadius( 00316 /* either 'i' or 't' for ionization or temperature */ 00317 char chType, 00318 /* atomic number on c scale */ 00319 long int nelem, 00320 /* this will say how many ion stages in arlog have non-zero values */ 00321 long int *n, 00322 /* array of values, log both cases, 00323 * hydrogen is special case since [2] will be H2 */ 00324 realnum arlog[], 00325 /* true, include electron density, false do not*/ 00326 bool lgDensity ) 00327 { 00328 long int ion, 00329 limit; 00330 double abund_radmean, 00331 normalize; 00332 00333 DEBUG_ENTRY( "MeanIonRadius()" ); 00334 00335 ASSERT( chType=='i' || chType=='t' ); 00336 00337 /* >>chng 04 jun 27, add H2 to top of H ladder, 00338 * in this routine nelem is on physical scale, so 1 for H, 00339 * limit is on C scale, such that ion<limit */ 00340 limit = MAX2( 3, nelem+2 ); 00341 00342 /* fills in array arlog with log of ionization fractions 00343 * N is set to highest stage with non-zero abundance 00344 * N set to 0 if element turned off 00345 * 00346 * first call MeanZero to sero out save arrays 00347 * next call MeanInc in zone calc to enter ionziation fractions or temperature 00348 * finally this routine computes actual means 00349 * */ 00350 if( !dense.lgElmtOn[nelem] ) 00351 { 00352 /* no ionization for this element */ 00353 /* >>chng 04 jun 27, upper limit had been <nelem+1, had missed fully 00354 * stipped ion */ 00355 for( ion=0; ion < limit; ion++ ) 00356 { 00357 arlog[ion] = -30.f; 00358 } 00359 *n = 0; 00360 return; 00361 } 00362 00363 /* set high ion stages, with zero abundances, to -30, 00364 * limit is on c scale, such that ion<=limit */ 00365 *n = limit; 00366 while( *n > 0 && mean.xIonMeans[0][nelem][*n-1] <= 0. ) 00367 { 00368 arlog[*n-1] = -30.f; 00369 --*n; 00370 } 00371 00372 if( chType=='i' && lgDensity) 00373 { 00374 /* mean ionization with density*/ 00375 for( ion=0; ion < *n; ion++ ) 00376 { 00377 if( mean.xIonEdenMeans[0][nelem][ion] > 0. ) 00378 { 00379 abund_radmean = mean.xIonEdenMeans[0][nelem][ion]; 00380 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean / mean.xIonEdenMeansNorm[0][nelem])); 00381 } 00382 else 00383 { 00384 arlog[ion] = -30.f; 00385 } 00386 } 00387 } 00388 else if( chType=='i' ) 00389 { 00390 /* mean ionization no density*/ 00391 for( ion=0; ion < *n; ion++ ) 00392 { 00393 if( mean.xIonMeans[0][nelem][ion] > 0. ) 00394 { 00395 abund_radmean = mean.xIonMeans[0][nelem][ion]; 00396 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/mean.xIonMeansNorm[0][nelem])); 00397 } 00398 else 00399 { 00400 arlog[ion] = -30.f; 00401 } 00402 } 00403 } 00404 else if( chType=='t' && lgDensity ) 00405 { 00406 /* mean temperature wrt elec density */ 00407 for( ion=0; ion < *n; ion++ ) 00408 { 00409 normalize = mean.TempEdenMeansNorm[0][nelem][ion]; 00410 if( normalize > SMALLFLOAT ) 00411 { 00412 abund_radmean = mean.TempEdenMeans[0][nelem][ion]; 00413 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/normalize)); 00414 } 00415 else 00416 { 00417 arlog[ion] = -30.f; 00418 } 00419 } 00420 } 00421 else if( chType=='t' ) 00422 { 00423 /* mean temperature without elec density*/ 00424 for( ion=0; ion < *n; ion++ ) 00425 { 00426 normalize = mean.TempMeansNorm[0][nelem][ion]; 00427 if( normalize > SMALLFLOAT ) 00428 { 00429 abund_radmean = mean.TempMeans[0][nelem][ion]; 00430 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/normalize)); 00431 } 00432 else 00433 { 00434 arlog[ion] = -30.f; 00435 } 00436 } 00437 } 00438 else 00439 { 00440 fprintf(ioQQQ," MeanIonRadius called with insane job \n"); 00441 } 00442 return; 00443 } 00444 00445 /*MeanIonVolume do volume mean of ionization fractions over volume of any element */ 00446 void MeanIonVolume( 00447 /* either 'i' or 't' for ionization or temperature */ 00448 char chType, 00449 /* atomic number on c scale */ 00450 long int nelem, 00451 /* this will say how many of arlog have non-zero values */ 00452 long int *n, 00453 /* array of values, log both cases */ 00454 realnum arlog[], 00455 /* true, include electron density, false do not*/ 00456 bool lgDensity ) 00457 { 00458 long int ion; 00459 double abund_volmean, 00460 normalize; 00461 long int limit; 00462 00463 DEBUG_ENTRY( "MeanIonVolume()" ); 00464 00465 ASSERT( chType=='i' || chType=='t' ); 00466 00467 /* >>chng 04 jun 27, add H2 to top of H ladder, 00468 * in this routine nelem is on physical scale, so 1 for H*/ 00469 limit = MAX2( 3, nelem+2 ); 00470 00471 /* 00472 * fills in array arlog with log of ionization fractions 00473 * N is set to highest stage with non-zero abundance 00474 * N set to 0 if element turned off 00475 * 00476 * first call MeanZero to sero out save arrays 00477 * next call MeanInc in zone calc to enter ionziation fractions 00478 * finally this routine computes actual means 00479 */ 00480 if( !dense.lgElmtOn[nelem] ) 00481 { 00482 /* no ionization for this element */ 00483 for( ion=0; ion <= limit; ion++ ) 00484 { 00485 arlog[ion] = -30.f; 00486 } 00487 *n = 0; 00488 return; 00489 } 00490 00491 /* this is number of stages of ionization */ 00492 *n = limit; 00493 /* fill in higher stages with zero if non-existent, 00494 * also decrement n, the number with non-zero abundances */ 00495 while( *n > 0 && mean.xIonMeans[1][nelem][*n-1] <= 0. ) 00496 { 00497 arlog[*n-1] = -30.f; 00498 --*n; 00499 } 00500 /* n is now the limit to the number with positive abundances */ 00501 00502 if( chType=='i' && lgDensity) 00503 { 00504 /* mean ionization with electron density*/ 00505 /* this is denominator for forming a mean */ 00506 /* return log of means */ 00507 for( ion=0; ion < *n; ion++ ) 00508 { 00509 if( mean.xIonEdenMeans[1][nelem][ion] > 0. ) 00510 { 00511 abund_volmean = mean.xIonEdenMeans[1][nelem][ion]; 00512 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean) / mean.xIonEdenMeansNorm[1][nelem]); 00513 } 00514 else 00515 { 00516 arlog[ion] = -30.f; 00517 } 00518 } 00519 } 00520 else if( chType=='i' ) 00521 { 00522 /* mean ionization */ 00523 /* this is denominator for forming a mean */ 00524 /* return log of means */ 00525 for( ion=0; ion < *n; ion++ ) 00526 { 00527 if( mean.xIonMeans[1][nelem][ion] > 0. ) 00528 { 00529 abund_volmean = mean.xIonMeans[1][nelem][ion]; 00530 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean) / mean.xIonMeansNorm[1][nelem]); 00531 } 00532 else 00533 { 00534 arlog[ion] = -30.f; 00535 } 00536 } 00537 } 00538 else if( chType=='t' && lgDensity ) 00539 { 00540 /* mean temperature with density*/ 00541 /* this is denominator for forming a mean */ 00542 /* return log of means */ 00543 for( ion=0; ion < *n; ion++ ) 00544 { 00545 normalize = mean.TempEdenMeansNorm[1][nelem][ion]; 00546 if( normalize > SMALLFLOAT ) 00547 { 00548 abund_volmean = mean.TempEdenMeans[1][nelem][ion]; 00549 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean)/normalize); 00550 } 00551 else 00552 { 00553 arlog[ion] = -30.f; 00554 } 00555 } 00556 } 00557 else if( chType=='t' ) 00558 { 00559 /* mean temperature with NO density*/ 00560 /* this is denominator for forming a mean */ 00561 /* return log of means */ 00562 for( ion=0; ion < *n; ion++ ) 00563 { 00564 normalize = mean.TempMeansNorm[1][nelem][ion]; 00565 if( normalize > SMALLFLOAT ) 00566 { 00567 abund_volmean = mean.TempMeans[1][nelem][ion]; 00568 arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean)/normalize); 00569 } 00570 else 00571 { 00572 arlog[ion] = -30.f; 00573 } 00574 } 00575 } 00576 else 00577 { 00578 fprintf(ioQQQ," MeanIonVolume called with insane job\n"); 00579 } 00580 return; 00581 } 00582 00583 /*aver compute average of various quantities over the computed geometry 00584 * called by startenditer to initialize, radinc to increment, and prtfinal for final results */ 00585 void aver( 00586 /* 4 char + null string giving job, prin, zero, zone, doit*/ 00587 const char *chWhat, 00588 /* the quantity to average, like the temperature */ 00589 double quan, 00590 /* what we weight against, like the o++ abundance */ 00591 double weight, 00592 /* 10 char + null string giving title for this average */ 00593 const char *chLabl) 00594 { 00595 /* NAVER is the limit to the number than can be averaged */ 00596 # define NAVER 20 00597 static char chLabavr[NAVER][11]; 00598 long int i, 00599 ioff, 00600 j; 00601 static long int n; 00602 double raver[NAVER]={0.}, 00603 vaver[NAVER]={0.}; 00604 static double aversv[4][NAVER]; 00605 00606 DEBUG_ENTRY( "aver()" ); 00607 00608 if( strcmp(chWhat,"zero") == 0 ) 00609 { 00610 /* chWhat='zero' - zero out the save array */ 00611 for( i=0; i < NAVER; i++ ) 00612 { 00613 for( j=0; j < 4; j++ ) 00614 { 00615 aversv[j][i] = 0.; 00616 } 00617 } 00618 n = 0; 00619 } 00620 else if( strcmp(chWhat,"zone") == 0 ) 00621 { 00622 n = 0; 00623 } 00624 else if( strcmp(chWhat,"doit") == 0 ) 00625 { 00626 /* chWhat='doit' - enter values to average */ 00627 if( n >= NAVER ) 00628 { 00629 fprintf( ioQQQ, " Too many values entered into AVER, increase NAVER\n" ); 00630 cdEXIT(EXIT_FAILURE); 00631 } 00632 aversv[0][n] += weight*quan*radius.drad_x_fillfac; 00633 aversv[1][n] += weight*radius.drad_x_fillfac; 00634 aversv[2][n] += weight*quan*radius.dVeff; 00635 aversv[3][n] += weight*radius.dVeff; 00636 00637 /* this is the label for this particular average, like " TeNe "*/ 00638 strcpy( chLabavr[n], chLabl ); 00639 n += 1; 00640 } 00641 else if( strcmp(chWhat,"prin") == 0 ) 00642 { 00643 /* set things up to get answers */ 00644 ioff = n*10/2 + 1; 00645 00646 /* space out to center the label on the numbers */ 00647 fprintf( ioQQQ,"\n"); 00648 for( i=0; i < ioff; i++ ) 00649 { 00650 /*chTitAvr[i] = ' ';*/ 00651 fprintf( ioQQQ, " " ); 00652 } 00653 00654 /* now print header with cr */ 00655 fprintf( ioQQQ, "Averaged Quantities\n " ); 00656 00657 fprintf( ioQQQ, " " ); 00658 for( i=0; i < n; i++ ) 00659 { 00660 fprintf( ioQQQ, "%10.10s", chLabavr[i] ); 00661 } 00662 fprintf( ioQQQ, " \n" ); 00663 00664 for( i=0; i < n; i++ ) 00665 { 00666 if( aversv[1][i] > 0. ) 00667 { 00668 raver[i] = aversv[0][i]/aversv[1][i]; 00669 } 00670 else 00671 { 00672 raver[i] = 0.; 00673 } 00674 if( aversv[3][i] > 0. ) 00675 { 00676 vaver[i] = aversv[2][i]/aversv[3][i]; 00677 } 00678 else 00679 { 00680 vaver[i] = 0.; 00681 } 00682 } 00683 00684 fprintf( ioQQQ, " Radius:" ); 00685 for( i=0; i < n; i++ ) 00686 { 00687 /*fprintf( ioQQQ, "%11.3e", raver[i] );*/ 00688 fprintf( ioQQQ, " "); 00689 fprintf(ioQQQ,PrintEfmt("%9.2e", raver[i] ) ); 00690 } 00691 fprintf( ioQQQ, "\n" ); 00692 00693 /* only print vol aver is lgSphere is set */ 00694 if( geometry.lgSphere ) 00695 { 00696 fprintf( ioQQQ, " Volume:" ); 00697 for( i=0; i < n; i++ ) 00698 { 00699 /*fprintf( ioQQQ, "%11.3e", vaver[i] );*/ 00700 fprintf( ioQQQ, " "); 00701 fprintf(ioQQQ,PrintEfmt("%9.2e", vaver[i] ) ); 00702 } 00703 fprintf( ioQQQ, "\n" ); 00704 } 00705 } 00706 else 00707 { 00708 fprintf( ioQQQ, " AVER does not understand chWhat=%4.4s\n", 00709 chWhat ); 00710 ShowMe(); 00711 cdEXIT(EXIT_FAILURE); 00712 } 00713 return; 00714 }