cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mean.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*MeanInc increment mean ionization fractions and temperatures over computed structure, in radius_increment */
4 /*MeanZero zero mean of ionization fractions array */
5 /*MeanIonRadius do radius mean ionization fractions or temperature over radius for any element */
6 /*MeanIonVolume do volume mean ionization fractions or temperature over volume for any element */
7 /*aver compute average of various quantities over the computed geometry
8  * called by startenditer to initialize, radinc to increment, and prtfinal for final results */
9 #include "cddefines.h"
10 #include "physconst.h"
11 #include "radius.h"
12 #include "dense.h"
13 #include "hyperfine.h"
14 #include "magnetic.h"
15 #include "hmi.h"
16 #include "phycon.h"
17 #include "geometry.h"
18 #include "mean.h"
19 
20 /*MeanInc increment mean ionization fractions and temperatures over computed structure, in radius_increment */
21 void MeanInc(void)
22 {
23  long int ion,
24  nelem;
25  double dc,
26  dcv;
27 
28  /* this routine is called by radius_increment during the calculation to update the
29  * sums that will become the vol and rad weighted mean ionizations */
30 
31  DEBUG_ENTRY( "MeanInc()" );
32 
33  for( nelem=0; nelem < LIMELM; nelem++ )
34  {
35  /* this is mean ionization */
36  if( nelem==ipHYDROGEN )
37  {
38  /* >>chng 04 jun 27, add this option, previously had not included
39  * H2 as part of total */
40  ion = 2;
41  /* hydrogen is special case since must include H2,
42  * and must mult by 2 to conserve total H - will need to div
43  * by two if column density ever used */
44  /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/
46  mean.xIonMeans[0][nelem][ion] += dc;
47  /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/
49  mean.xIonEdenMeans[0][nelem][ion] += dc;
50  /* xIonMeans[1][n][ion] is vol integration */
51  dcv = 2.*hmi.H2_total*radius.dVeff;
52  mean.xIonMeans[1][nelem][ion] += dcv;
53  /* xIonMeans[1][ion][n] is vol integration PLUS electron density */
54  dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden;
55  mean.xIonEdenMeans[1][nelem][ion] += dcv;
56  }
57  for( ion=0; ion < (nelem + 2); ion++ )
58  {
59  /* xIonMeans[0][ion][n] is radial integration */
60  dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac;
61  mean.xIonMeans[0][nelem][ion] += dc;
62 
63  /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/
64  dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden;
65  mean.xIonEdenMeans[0][nelem][ion] += dc;
66 
67  /* xIonMeans[1][n][ion] is vol integration */
68  /* >>chng 00 mar 28, r1r0sq had been dVeff,
69  * bug discovered by Valentina Luridiana */
70  dcv = dense.xIonDense[nelem][ion]*radius.dVeff;
71  mean.xIonMeans[1][nelem][ion] += dcv;
72 
73  /* xIonMeans[1][ion][n] is vol integration PLUS electron density */
74  /* >>chng 00 mar 28, r1r0sq had been dVeff,
75  * bug discovered by Valentina Luridiana */
76  dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden;
77  mean.xIonEdenMeans[1][nelem][ion] += dcv;
78  }
79  /* these normalize the above, use gas_phase which includes molecular parts */
80  /* first two are means over radius */
82  mean.xIonMeansNorm[0][nelem] += dc;
84  mean.xIonEdenMeansNorm[0][nelem] += dc;
85  /* next two means over volume */
86  /* >>chng 04 jun 28, changed radius.dVeff to dVeff,
87  * which is equivalent in thin zone limit, more accurate with thick zones, as per
88  * Valentina Luridiana email */
89  /*dcv = dense.gas_phase[nelem]*radius.dVeff;*/
90  dcv = dense.gas_phase[nelem]*radius.dVeff;
91  mean.xIonMeansNorm[1][nelem] += dcv;
92  /*dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden;*/
93  dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden;
94  mean.xIonEdenMeansNorm[1][nelem] += dcv;
95 
96  /* this is mean temperature */
97  /* this is ionization on the IonFracs scale, offset 1 up since
98  * zero is total abundance. This works well since the mean
99  * ionization array xIonMeans is also offset up one, since
100  * [0] is the normalization */
101  if( nelem==ipHYDROGEN )
102  {
103  ion = 2;
104  /* TempMeans[0][ion][n] is radial integration */
106  mean.TempMeans[0][nelem][ion] += dc*phycon.te;
107  mean.TempMeansNorm[0][nelem][ion] += dc;
108 
109  /* TempMeans[0][ion][n] is radial integration PLUS electron density*/
111  mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
112  mean.TempEdenMeansNorm[0][nelem][ion] += dc;
113 
114  /* TempMeans[1][ion+1][n] is vol integration */
115  /*>>chng 00 dec 18, following had dVeff rather than r1r0sq */
116  dcv = 2.*hmi.H2_total*radius.dVeff;
117  mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
118  mean.TempMeansNorm[1][nelem][ion] += dcv;
119 
120  /* TempMeans[1][ion+1][n] is vol integration PLUS electron density */
121  dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden;
122  mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te;
123  mean.TempEdenMeansNorm[1][nelem][ion] += dcv;
124  }
125  for( ion=0; ion < (nelem + 2); ion++ )
126  {
127  /* TempMeans[0][ion][n] is radial integration */
128  dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac;
129  mean.TempMeans[0][nelem][ion] += dc*phycon.te;
130  mean.TempMeansNorm[0][nelem][ion] += dc;
131 
132  /* TempMeans[0][ion][n] is radial integration PLUS electron density*/
133  dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden;
134  mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
135  mean.TempEdenMeansNorm[0][nelem][ion] += dc;
136 
137  /* TempMeans[1][ion+1][n] is vol integration */
138  /*>>chng 00 dec 18, following had dVeff rather than r1r0sq */
139  dcv = dense.xIonDense[nelem][ion]*radius.dVeff;
140  mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
141  mean.TempMeansNorm[1][nelem][ion] += dcv;
142 
143  /* TempMeans[1][ion+1][n] is vol integration PLUS electron density */
144  dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden;
145  mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te;
146  mean.TempEdenMeansNorm[1][nelem][ion] += dcv;
147  }
148  }
149 
150  /* =============================================================
151  *
152  * these are some special quanties for the mean
153  */
154 
155  /* >>chng 05 dec 28, add this
156  * used to get magnetic field weighted wrt harmonic mean spin temperature,
157  * * H0 over radius - as measured by 21cm temperature */
159  {
161  }
162  else
163  {
164  dc = 0.;
165  }
166  /* mean magnetic field weighted wrt 21 cm opacity, n(H0)/T */
167  mean.B_HarMeanTempRadius[0] += sqrt(fabs(magnetic.pressure)*PI8) * dc;
168  /* this assumes that field is tangled */
169  mean.B_HarMeanTempRadius[1] += dc;
170 
171  /* used to get harmonic mean temperature with respect to H over radius,
172  * for comparison with 21cm temperature */
174 
175  /* harmonic mean gas kinetic temp, over radius, for comparison with 21 cm obs */
176  /*HEADS UP - next four are the inverse of equation 3 of
177  * >>refer Tspin 21 cm Abel, N.P., Brogan, C.L., Ferland, G.J., O'Dell, C.R.,
178  * >>refercon Shaw, G. & Troland, T.H. 2004, ApJ, 609, 247 */
179  mean.HarMeanTempRadius[0] += dc;
180  mean.HarMeanTempRadius[1] += dc/phycon.te;
181 
182  /* harmonic mean gas kinetic temp, over volume, for symmetry with above, for comp with 21 cm obs */
185 
186  /* harmonic mean of computed 21 cm spin temperature - this is what 21 cm actually measures */
187  mean.H_21cm_spin_mean_radius[0] += dc;
189 
190  /* mean H2 temperature over radius */
192  mean.H2MeanTempRadius[0] += dc*phycon.te;
193  mean.H2MeanTempRadius[1] += dc;
194 
195  /* mean H2 temperature over volume, for symmetry */
196  dcv = hmi.H2_total*radius.dVeff;
197  mean.H2MeanTempVolume[0] += dc*phycon.te;
198  mean.H2MeanTempVolume[1] += dc;
199 
200  /* >>chng 05 dec 28, add mean temp over radius and vol */
201  dc = radius.drad_x_fillfac;
202  mean.TempMeanRadius[0] += dc*phycon.te;
203  mean.TempMeanRadius[1] += dc;
204 
205  dcv = radius.dVeff;
206  mean.TempMeanVolume[0] += dc*phycon.te;
207  mean.TempMeanVolume[1] += dc;
208  return;
209 }
210 
211 /*MeanZero zero mean of ionization fractions array */
212 void MeanZero(void)
213 {
214  long int ion,
215  j,
216  nelem,
217  limit;
218  static bool lgFirst=true;
219 
220  DEBUG_ENTRY( "MeanZero()" );
221 
222  /*
223  * MeanZero is called at start of calculation by zero, and by
224  * startenditer to initialize the variables
225  */
226 
227  if( lgFirst )
228  {
229  lgFirst = false;
230  /* both are [2], [nelem][ion] */
231  mean.xIonMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
232  mean.xIonEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
233  mean.xIonMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) );
234  mean.xIonEdenMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) );
235  mean.TempMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
236  mean.TempEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
237  mean.TempMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
238  mean.TempEdenMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
239  for( j=0; j<2; ++j )
240  {
241  mean.xIonMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
242  mean.xIonEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
243  mean.xIonMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) );
244  mean.xIonEdenMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) );
245  mean.TempMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
246  mean.TempEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
247  mean.TempMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
248  mean.TempEdenMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
249  for( nelem=0; nelem<LIMELM; ++nelem )
250  {
251  limit = MAX2(3,nelem+2);
252  mean.xIonMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
253  mean.xIonEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
254  mean.TempMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
255  mean.TempEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
256  mean.TempMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
257  mean.TempEdenMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
258  }
259  }
260  }
261 
262  for( j=0; j < 2; j++ )
263  {
264  for( nelem=0; nelem < LIMELM; nelem++ )
265  {
266  mean.xIonMeansNorm[j][nelem] = 0.;
267  mean.xIonEdenMeansNorm[j][nelem] = 0.;
268  /* >>chng 04 jun 27, H2 will be 3rd ion stage of H */
269  limit = MAX2(3,nelem+2);
270  /*for( ion=0; ion <= (nelem + 2); ion++ )*/
271  for( ion=0; ion < limit; ion++ )
272  {
273  mean.TempMeansNorm[j][nelem][ion] = 0.;
274  mean.TempEdenMeansNorm[j][nelem][ion] = 0.;
275  /* these are used to save info on temperature and ionization means */
276  mean.xIonMeans[j][nelem][ion] = 0.;
277  mean.xIonEdenMeans[j][nelem][ion] = 0.;
278  /* double here since [2] and [3] have norm for average */
279  mean.TempMeans[j][nelem][ion] = 0.;
280  mean.TempEdenMeans[j][nelem][ion] = 0.;
281  }
282  }
283  }
284 
285  /* mean over radius, for comp with 21 cm obs */
286  mean.HarMeanTempRadius[0] = 0.;
287  mean.HarMeanTempRadius[1] = 0.;
288 
289  /* mean over volume, for symmetry */
290  mean.HarMeanTempVolume[0] = 0.;
291  mean.HarMeanTempVolume[1] = 0.;
292 
293  /* mena of calculated 21 cm spin temperature */
296 
297  /* H2 mean temp over radius */
298  mean.H2MeanTempRadius[0] = 0.;
299  mean.H2MeanTempRadius[1] = 0.;
300  /* H2 mean temp over volume */
301  mean.H2MeanTempVolume[0] = 0.;
302  mean.H2MeanTempVolume[1] = 0.;
303 
304  mean.TempMeanRadius[0] = 0.;
305  mean.TempMeanRadius[1] = 0.;
306  mean.TempMeanVolume[0] = 0.;
307  mean.TempMeanVolume[1] = 0.;
308 
309  mean.B_HarMeanTempRadius[0] = 0.;
310  mean.B_HarMeanTempRadius[1] = 0.;
311  return;
312 }
313 
314 /*MeanIonRadius derive mean ionization fractions over radius for some element */
316  /* either 'i' or 't' for ionization or temperature */
317  char chType,
318  /* atomic number on c scale */
319  long int nelem,
320  /* this will say how many ion stages in arlog have non-zero values */
321  long int *n,
322  /* array of values, log both cases,
323  * hydrogen is special case since [2] will be H2 */
324  realnum arlog[],
325  /* true, include electron density, false do not*/
326  bool lgDensity )
327 {
328  long int ion,
329  limit;
330  double abund_radmean,
331  normalize;
332 
333  DEBUG_ENTRY( "MeanIonRadius()" );
334 
335  ASSERT( chType=='i' || chType=='t' );
336 
337  /* >>chng 04 jun 27, add H2 to top of H ladder,
338  * in this routine nelem is on physical scale, so 1 for H,
339  * limit is on C scale, such that ion<limit */
340  limit = MAX2( 3, nelem+2 );
341 
342  /* fills in array arlog with log of ionization fractions
343  * N is set to highest stage with non-zero abundance
344  * N set to 0 if element turned off
345  *
346  * first call MeanZero to sero out save arrays
347  * next call MeanInc in zone calc to enter ionziation fractions or temperature
348  * finally this routine computes actual means
349  * */
350  if( !dense.lgElmtOn[nelem] )
351  {
352  /* no ionization for this element */
353  /* >>chng 04 jun 27, upper limit had been <nelem+1, had missed fully
354  * stipped ion */
355  for( ion=0; ion < limit; ion++ )
356  {
357  arlog[ion] = -30.f;
358  }
359  *n = 0;
360  return;
361  }
362 
363  /* set high ion stages, with zero abundances, to -30,
364  * limit is on c scale, such that ion<=limit */
365  *n = limit;
366  while( *n > 0 && mean.xIonMeans[0][nelem][*n-1] <= 0. )
367  {
368  arlog[*n-1] = -30.f;
369  --*n;
370  }
371 
372  if( chType=='i' && lgDensity)
373  {
374  /* mean ionization with density*/
375  for( ion=0; ion < *n; ion++ )
376  {
377  if( mean.xIonEdenMeans[0][nelem][ion] > 0. )
378  {
379  abund_radmean = mean.xIonEdenMeans[0][nelem][ion];
380  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean / mean.xIonEdenMeansNorm[0][nelem]));
381  }
382  else
383  {
384  arlog[ion] = -30.f;
385  }
386  }
387  }
388  else if( chType=='i' )
389  {
390  /* mean ionization no density*/
391  for( ion=0; ion < *n; ion++ )
392  {
393  if( mean.xIonMeans[0][nelem][ion] > 0. )
394  {
395  abund_radmean = mean.xIonMeans[0][nelem][ion];
396  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/mean.xIonMeansNorm[0][nelem]));
397  }
398  else
399  {
400  arlog[ion] = -30.f;
401  }
402  }
403  }
404  else if( chType=='t' && lgDensity )
405  {
406  /* mean temperature wrt elec density */
407  for( ion=0; ion < *n; ion++ )
408  {
409  normalize = mean.TempEdenMeansNorm[0][nelem][ion];
410  if( normalize > SMALLFLOAT )
411  {
412  abund_radmean = mean.TempEdenMeans[0][nelem][ion];
413  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/normalize));
414  }
415  else
416  {
417  arlog[ion] = -30.f;
418  }
419  }
420  }
421  else if( chType=='t' )
422  {
423  /* mean temperature without elec density*/
424  for( ion=0; ion < *n; ion++ )
425  {
426  normalize = mean.TempMeansNorm[0][nelem][ion];
427  if( normalize > SMALLFLOAT )
428  {
429  abund_radmean = mean.TempMeans[0][nelem][ion];
430  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/normalize));
431  }
432  else
433  {
434  arlog[ion] = -30.f;
435  }
436  }
437  }
438  else
439  {
440  fprintf(ioQQQ," MeanIonRadius called with insane job \n");
441  }
442  return;
443 }
444 
445 /*MeanIonVolume do volume mean of ionization fractions over volume of any element */
447  /* either 'i' or 't' for ionization or temperature */
448  char chType,
449  /* atomic number on c scale */
450  long int nelem,
451  /* this will say how many of arlog have non-zero values */
452  long int *n,
453  /* array of values, log both cases */
454  realnum arlog[],
455  /* true, include electron density, false do not*/
456  bool lgDensity )
457 {
458  long int ion;
459  double abund_volmean,
460  normalize;
461  long int limit;
462 
463  DEBUG_ENTRY( "MeanIonVolume()" );
464 
465  ASSERT( chType=='i' || chType=='t' );
466 
467  /* >>chng 04 jun 27, add H2 to top of H ladder,
468  * in this routine nelem is on physical scale, so 1 for H*/
469  limit = MAX2( 3, nelem+2 );
470 
471  /*
472  * fills in array arlog with log of ionization fractions
473  * N is set to highest stage with non-zero abundance
474  * N set to 0 if element turned off
475  *
476  * first call MeanZero to sero out save arrays
477  * next call MeanInc in zone calc to enter ionziation fractions
478  * finally this routine computes actual means
479  */
480  if( !dense.lgElmtOn[nelem] )
481  {
482  /* no ionization for this element */
483  for( ion=0; ion <= limit; ion++ )
484  {
485  arlog[ion] = -30.f;
486  }
487  *n = 0;
488  return;
489  }
490 
491  /* this is number of stages of ionization */
492  *n = limit;
493  /* fill in higher stages with zero if non-existent,
494  * also decrement n, the number with non-zero abundances */
495  while( *n > 0 && mean.xIonMeans[1][nelem][*n-1] <= 0. )
496  {
497  arlog[*n-1] = -30.f;
498  --*n;
499  }
500  /* n is now the limit to the number with positive abundances */
501 
502  if( chType=='i' && lgDensity)
503  {
504  /* mean ionization with electron density*/
505  /* this is denominator for forming a mean */
506  /* return log of means */
507  for( ion=0; ion < *n; ion++ )
508  {
509  if( mean.xIonEdenMeans[1][nelem][ion] > 0. )
510  {
511  abund_volmean = mean.xIonEdenMeans[1][nelem][ion];
512  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean) / mean.xIonEdenMeansNorm[1][nelem]);
513  }
514  else
515  {
516  arlog[ion] = -30.f;
517  }
518  }
519  }
520  else if( chType=='i' )
521  {
522  /* mean ionization */
523  /* this is denominator for forming a mean */
524  /* return log of means */
525  for( ion=0; ion < *n; ion++ )
526  {
527  if( mean.xIonMeans[1][nelem][ion] > 0. )
528  {
529  abund_volmean = mean.xIonMeans[1][nelem][ion];
530  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean) / mean.xIonMeansNorm[1][nelem]);
531  }
532  else
533  {
534  arlog[ion] = -30.f;
535  }
536  }
537  }
538  else if( chType=='t' && lgDensity )
539  {
540  /* mean temperature with density*/
541  /* this is denominator for forming a mean */
542  /* return log of means */
543  for( ion=0; ion < *n; ion++ )
544  {
545  normalize = mean.TempEdenMeansNorm[1][nelem][ion];
546  if( normalize > SMALLFLOAT )
547  {
548  abund_volmean = mean.TempEdenMeans[1][nelem][ion];
549  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean)/normalize);
550  }
551  else
552  {
553  arlog[ion] = -30.f;
554  }
555  }
556  }
557  else if( chType=='t' )
558  {
559  /* mean temperature with NO density*/
560  /* this is denominator for forming a mean */
561  /* return log of means */
562  for( ion=0; ion < *n; ion++ )
563  {
564  normalize = mean.TempMeansNorm[1][nelem][ion];
565  if( normalize > SMALLFLOAT )
566  {
567  abund_volmean = mean.TempMeans[1][nelem][ion];
568  arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean)/normalize);
569  }
570  else
571  {
572  arlog[ion] = -30.f;
573  }
574  }
575  }
576  else
577  {
578  fprintf(ioQQQ," MeanIonVolume called with insane job\n");
579  }
580  return;
581 }
582 
583 /*aver compute average of various quantities over the computed geometry
584  * called by startenditer to initialize, radinc to increment, and prtfinal for final results */
585 void aver(
586  /* 4 char + null string giving job, prin, zero, zone, doit*/
587  const char *chWhat,
588  /* the quantity to average, like the temperature */
589  double quan,
590  /* what we weight against, like the o++ abundance */
591  double weight,
592  /* 10 char + null string giving title for this average */
593  const char *chLabl)
594 {
595  /* NAVER is the limit to the number than can be averaged */
596 # define NAVER 20
597  static char chLabavr[NAVER][11];
598  long int i,
599  ioff,
600  j;
601  static long int n;
602  double raver[NAVER]={0.},
603  vaver[NAVER]={0.};
604  static double aversv[4][NAVER];
605 
606  DEBUG_ENTRY( "aver()" );
607 
608  if( strcmp(chWhat,"zero") == 0 )
609  {
610  /* chWhat='zero' - zero out the save array */
611  for( i=0; i < NAVER; i++ )
612  {
613  for( j=0; j < 4; j++ )
614  {
615  aversv[j][i] = 0.;
616  }
617  }
618  n = 0;
619  }
620  else if( strcmp(chWhat,"zone") == 0 )
621  {
622  n = 0;
623  }
624  else if( strcmp(chWhat,"doit") == 0 )
625  {
626  /* chWhat='doit' - enter values to average */
627  if( n >= NAVER )
628  {
629  fprintf( ioQQQ, " Too many values entered into AVER, increase NAVER\n" );
630  cdEXIT(EXIT_FAILURE);
631  }
632  aversv[0][n] += weight*quan*radius.drad_x_fillfac;
633  aversv[1][n] += weight*radius.drad_x_fillfac;
634  aversv[2][n] += weight*quan*radius.dVeff;
635  aversv[3][n] += weight*radius.dVeff;
636 
637  /* this is the label for this particular average, like " TeNe "*/
638  strcpy( chLabavr[n], chLabl );
639  n += 1;
640  }
641  else if( strcmp(chWhat,"prin") == 0 )
642  {
643  /* set things up to get answers */
644  ioff = n*10/2 + 1;
645 
646  /* space out to center the label on the numbers */
647  fprintf( ioQQQ,"\n");
648  for( i=0; i < ioff; i++ )
649  {
650  /*chTitAvr[i] = ' ';*/
651  fprintf( ioQQQ, " " );
652  }
653 
654  /* now print header with cr */
655  fprintf( ioQQQ, "Averaged Quantities\n " );
656 
657  fprintf( ioQQQ, " " );
658  for( i=0; i < n; i++ )
659  {
660  fprintf( ioQQQ, "%10.10s", chLabavr[i] );
661  }
662  fprintf( ioQQQ, " \n" );
663 
664  for( i=0; i < n; i++ )
665  {
666  if( aversv[1][i] > 0. )
667  {
668  raver[i] = aversv[0][i]/aversv[1][i];
669  }
670  else
671  {
672  raver[i] = 0.;
673  }
674  if( aversv[3][i] > 0. )
675  {
676  vaver[i] = aversv[2][i]/aversv[3][i];
677  }
678  else
679  {
680  vaver[i] = 0.;
681  }
682  }
683 
684  fprintf( ioQQQ, " Radius:" );
685  for( i=0; i < n; i++ )
686  {
687  /*fprintf( ioQQQ, "%11.3e", raver[i] );*/
688  fprintf( ioQQQ, " ");
689  fprintf(ioQQQ,PrintEfmt("%9.2e", raver[i] ) );
690  }
691  fprintf( ioQQQ, "\n" );
692 
693  /* only print vol aver is lgSphere is set */
694  if( geometry.lgSphere )
695  {
696  fprintf( ioQQQ, " Volume:" );
697  for( i=0; i < n; i++ )
698  {
699  /*fprintf( ioQQQ, "%11.3e", vaver[i] );*/
700  fprintf( ioQQQ, " ");
701  fprintf(ioQQQ,PrintEfmt("%9.2e", vaver[i] ) );
702  }
703  fprintf( ioQQQ, "\n" );
704  }
705  }
706  else
707  {
708  fprintf( ioQQQ, " AVER does not understand chWhat=%4.4s\n",
709  chWhat );
710  ShowMe();
711  cdEXIT(EXIT_FAILURE);
712  }
713  return;
714 }

Generated for cloudy by doxygen 1.8.4