/home66/gary/public_html/cloudy/c08_branch/source/mean.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*MeanInc increment mean ionization fractions and temperatures over computed structure, in radius_increment */
/*MeanZero zero mean of ionization fractions array */
/*MeanIonRadius do radius mean ionization fractions or temperature over radius for any element */
/*MeanIonVolume do volume mean ionization fractions or temperature over volume for any element */
/*aver compute average of various quantities over the computed geometry
 * called by startenditer to initialize, radinc to increment, and prtfinal for final results */
#include "cddefines.h"
#include "physconst.h"
#include "radius.h"
#include "dense.h"
#include "hyperfine.h"
#include "magnetic.h"
#include "hmi.h"
#include "phycon.h"
#include "geometry.h"
#include "mean.h"

/*MeanInc increment mean ionization fractions and temperatures over computed structure, in radius_increment */
void MeanInc(void)
{
        long int ion, 
          nelem;
        double dc, 
          dcv;

        /* this routine is called by radius_increment during the calculation to update the 
         * sums that will become the vol and rad weighted mean ionizations */

        DEBUG_ENTRY( "MeanInc()" );

        for( nelem=0; nelem < LIMELM; nelem++ )
        {
                /* this is mean ionization */
                if( nelem==ipHYDROGEN )
                {
                        /* >>chng 04 jun 27, add this option, previously had not included
                         * H2 as part of total */
                        ion = 2;
                        /* hydrogen is special case since must include H2,
                         * and must mult by 2 to conserve total H - will need to div
                         * by two if column density ever used */
                        /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/
                        dc = 2.*hmi.H2_total*radius.drad_x_fillfac;
                        mean.xIonMeans[0][nelem][ion] += dc;
                        /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/
                        dc = 2.*hmi.H2_total*radius.drad_x_fillfac*dense.eden;
                        mean.xIonEdenMeans[0][nelem][ion] += dc;
                        /* xIonMeans[1][n][ion] is vol integration */
                        dcv = 2.*hmi.H2_total*radius.dVeff;
                        mean.xIonMeans[1][nelem][ion] += dcv;
                        /* xIonMeans[1][ion][n] is vol integration PLUS electron density */
                        dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden;
                        mean.xIonEdenMeans[1][nelem][ion] += dcv;
                }
                for( ion=0; ion < (nelem + 2); ion++ )
                {
                        /* xIonMeans[0][ion][n] is radial integration */
                        dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac;
                        mean.xIonMeans[0][nelem][ion] += dc;

                        /* xIonEdenMeans[0][ion][n] is radial integration PLUS electron density*/
                        dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden;
                        mean.xIonEdenMeans[0][nelem][ion] += dc;

                        /* xIonMeans[1][n][ion] is vol integration */
                        /* >>chng 00 mar 28, r1r0sq had been dVeff,
                         * bug discovered by Valentina Luridiana */
                        dcv = dense.xIonDense[nelem][ion]*radius.dVeff;
                        mean.xIonMeans[1][nelem][ion] += dcv;

                        /* xIonMeans[1][ion][n] is vol integration PLUS electron density */
                        /* >>chng 00 mar 28, r1r0sq had been dVeff,
                         * bug discovered by Valentina Luridiana */
                        dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden;
                        mean.xIonEdenMeans[1][nelem][ion] += dcv;
                }
                /* these normalize the above, use gas_phase which includes molecular parts */
                /* first two are means over radius */
                dc = dense.gas_phase[nelem]*radius.drad_x_fillfac;
                mean.xIonMeansNorm[0][nelem] += dc;
                dc = dense.gas_phase[nelem]*radius.drad_x_fillfac*dense.eden;
                mean.xIonEdenMeansNorm[0][nelem] += dc;
                /* next two means over volume */
                /* >>chng 04 jun 28, changed radius.dVeff to dVeff,
                 * which is equivalent in thin zone limit, more accurate with thick zones, as per
                 * Valentina Luridiana email */
                /*dcv = dense.gas_phase[nelem]*radius.dVeff;*/
                dcv = dense.gas_phase[nelem]*radius.dVeff;
                mean.xIonMeansNorm[1][nelem] += dcv;
                /*dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden;*/
                dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden;
                mean.xIonEdenMeansNorm[1][nelem] += dcv;

                /* this is mean temperature */
                /* this is ionization on the IonFracs scale, offset 1 up since
                 * zero is total abundance.  This works well since the mean
                 * ionization array xIonMeans is also offset up one, since 
                 * [0] is the normalization */
                if( nelem==ipHYDROGEN )
                {
                        ion = 2;
                        /* TempMeans[0][ion][n] is radial integration */
                        dc = 2.*hmi.H2_total*radius.drad_x_fillfac;
                        mean.TempMeans[0][nelem][ion] += dc*phycon.te;
                        mean.TempMeansNorm[0][nelem][ion] += dc;

                        /* TempMeans[0][ion][n] is radial integration PLUS electron density*/
                        dc = 2.*hmi.H2_total*radius.drad_x_fillfac*dense.eden;
                        mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
                        mean.TempEdenMeansNorm[0][nelem][ion] += dc;

                        /* TempMeans[1][ion+1][n] is vol integration */
                        /*>>chng 00 dec 18, following had dVeff rather than r1r0sq */
                        dcv = 2.*hmi.H2_total*radius.dVeff;
                        mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
                        mean.TempMeansNorm[1][nelem][ion] += dcv;

                        /* TempMeans[1][ion+1][n] is vol integration PLUS electron density */
                        dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden;
                        mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te;
                        mean.TempEdenMeansNorm[1][nelem][ion] += dcv;
                }
                for( ion=0; ion < (nelem + 2); ion++ )
                {
                        /* TempMeans[0][ion][n] is radial integration */
                        dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac;
                        mean.TempMeans[0][nelem][ion] += dc*phycon.te;
                        mean.TempMeansNorm[0][nelem][ion] += dc;

                        /* TempMeans[0][ion][n] is radial integration PLUS electron density*/
                        dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden;
                        mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
                        mean.TempEdenMeansNorm[0][nelem][ion] += dc;

                        /* TempMeans[1][ion+1][n] is vol integration */
                        /*>>chng 00 dec 18, following had dVeff rather than r1r0sq */
                        dcv = dense.xIonDense[nelem][ion]*radius.dVeff;
                        mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
                        mean.TempMeansNorm[1][nelem][ion] += dcv;

                        /* TempMeans[1][ion+1][n] is vol integration PLUS electron density */
                        dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden;
                        mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te;
                        mean.TempEdenMeansNorm[1][nelem][ion] += dcv;
                }
        }

        /* =============================================================
         *
         * these are some special quanties for the mean 
         */

        /* >>chng 05 dec 28, add this
         * used to get magnetic field weighted wrt harmonic mean spin temperature, 
         * * H0 over radius - as measured by 21cm temperature */
        if( hyperfine.Tspin21cm > SMALLFLOAT )
        {
                dc = dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac/phycon.te;
        }
        else
        {
                dc = 0.;
        }
        /* mean magnetic field weighted wrt 21 cm opacity, n(H0)/T */
        mean.B_HarMeanTempRadius[0] += sqrt(fabs(magnetic.pressure)*PI8) * dc;
        /* this assumes that field is tangled */
        mean.B_HarMeanTempRadius[1] += dc;

        /* used to get harmonic mean temperature with respect to H over radius,
         * for comparison with 21cm temperature */
        dc = dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac;

        /* harmonic mean gas kinetic temp, over radius, for comparison with 21 cm obs */
        /*HEADS UP - next four are the inverse of equation 3 of
         * >>refer      Tspin   21 cm   Abel, N.P., Brogan, C.L., Ferland, G.J., O'Dell, C.R., 
         * >>refercon   Shaw, G. & Troland, T.H. 2004, ApJ, 609, 247 */
        mean.HarMeanTempRadius[0] += dc;
        mean.HarMeanTempRadius[1] += dc/phycon.te;

        /* harmonic mean gas kinetic temp, over volume, for symmetry with above, for comp with 21 cm obs  */
        mean.HarMeanTempVolume[0] += dc*radius.r1r0sq;
        mean.HarMeanTempVolume[1] += dc/phycon.te*radius.r1r0sq;

        /* harmonic mean of computed 21 cm spin temperature - this is what 21 cm actually measures */
        mean.H_21cm_spin_mean_radius[0] += dc;
        mean.H_21cm_spin_mean_radius[1] += dc / SDIV( hyperfine.Tspin21cm );

        /* mean H2 temperature over radius */
        dc = hmi.H2_total*radius.drad_x_fillfac;
        mean.H2MeanTempRadius[0] += dc*phycon.te;
        mean.H2MeanTempRadius[1] += dc;

        /* mean H2 temperature over volume, for symmetry */
        dcv = hmi.H2_total*radius.dVeff;
        mean.H2MeanTempVolume[0] += dc*phycon.te;
        mean.H2MeanTempVolume[1] += dc;

        /* >>chng 05 dec 28, add mean temp over radius and vol */
        dc = radius.drad_x_fillfac;
        mean.TempMeanRadius[0] += dc*phycon.te;
        mean.TempMeanRadius[1] += dc;

        dcv = radius.dVeff;
        mean.TempMeanVolume[0] += dc*phycon.te;
        mean.TempMeanVolume[1] += dc;
        return;
}

/*MeanZero zero mean of ionization fractions array */
void MeanZero(void)
{
        long int ion, 
          j, 
          nelem,
          limit;
        static bool lgFirst=true;

        DEBUG_ENTRY( "MeanZero()" );

        /* 
         * MeanZero is called at start of calculation by zero, and by
         * startenditer to initialize the variables 
         */

        if( lgFirst )
        {
                lgFirst = false;
                /* both are [2], [nelem][ion] */
                mean.xIonMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
                mean.xIonEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
                mean.xIonMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) );
                mean.xIonEdenMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) );
                mean.TempMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
                mean.TempEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
                mean.TempMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
                mean.TempEdenMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
                for( j=0; j<2; ++j )
                {
                        mean.xIonMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
                        mean.xIonEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
                        mean.xIonMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) );
                        mean.xIonEdenMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) );
                        mean.TempMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
                        mean.TempEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
                        mean.TempMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
                        mean.TempEdenMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
                        for( nelem=0; nelem<LIMELM; ++nelem )
                        {
                                limit = MAX2(3,nelem+2);
                                mean.xIonMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
                                mean.xIonEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
                                mean.TempMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
                                mean.TempEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
                                mean.TempMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
                                mean.TempEdenMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
                        }
                }
        }

        for( j=0; j < 2; j++ )
        {
                for( nelem=0; nelem < LIMELM; nelem++ )
                {
                        mean.xIonMeansNorm[j][nelem] = 0.;
                        mean.xIonEdenMeansNorm[j][nelem] = 0.;
                        /* >>chng 04 jun 27, H2 will be 3rd ion stage of H */
                        limit = MAX2(3,nelem+2);
                        /*for( ion=0; ion <= (nelem + 2); ion++ )*/
                        for( ion=0; ion < limit; ion++ )
                        {
                                mean.TempMeansNorm[j][nelem][ion] = 0.;
                                mean.TempEdenMeansNorm[j][nelem][ion] = 0.;
                                /* these are used to save info on temperature and ionization means */
                                mean.xIonMeans[j][nelem][ion] = 0.;
                                mean.xIonEdenMeans[j][nelem][ion] = 0.;
                                /* double here since [2] and [3] have norm for average */
                                mean.TempMeans[j][nelem][ion] = 0.;
                                mean.TempEdenMeans[j][nelem][ion] = 0.;
                        }
                }
        }

        /* mean over radius, for comp with 21 cm obs */
        mean.HarMeanTempRadius[0] = 0.;
        mean.HarMeanTempRadius[1] = 0.;

        /* mean over volume, for symmetry */
        mean.HarMeanTempVolume[0] = 0.;
        mean.HarMeanTempVolume[1] = 0.;

        /* mena of calculated 21 cm spin temperature */
        mean.H_21cm_spin_mean_radius[0] = 0.;
        mean.H_21cm_spin_mean_radius[1] = 0.;

        /* H2 mean temp over radius */
        mean.H2MeanTempRadius[0] = 0.;
        mean.H2MeanTempRadius[1] = 0.;
        /* H2 mean temp over volume */
        mean.H2MeanTempVolume[0] = 0.;
        mean.H2MeanTempVolume[1] = 0.;

        mean.TempMeanRadius[0] = 0.;
        mean.TempMeanRadius[1] = 0.;
        mean.TempMeanVolume[0] = 0.;
        mean.TempMeanVolume[1] = 0.;

        mean.B_HarMeanTempRadius[0] = 0.;
        mean.B_HarMeanTempRadius[1] = 0.;
        return;
}

/*MeanIonRadius derive mean ionization fractions over radius for some element */
void MeanIonRadius(
        /* either 'i' or 't' for ionization or temperature */
        char chType,
        /* atomic number on c scale */
        long int nelem, 
        /* this will say how many ion stages in arlog have non-zero values */
        long int *n, 
        /* array of values, log both cases, 
         * hydrogen is special case since [2] will be H2  */
        realnum arlog[],
        /* true, include electron density, false do not*/
        bool lgDensity )
{
        long int ion,
                limit;
        double abund_radmean, 
          normalize;

        DEBUG_ENTRY( "MeanIonRadius()" );

        ASSERT( chType=='i' || chType=='t' );

        /* >>chng 04 jun 27, add H2 to top of H ladder,
         * in this routine nelem is on physical scale, so 1 for H,
         * limit is on C scale, such that ion<limit */
        limit = MAX2( 3, nelem+2 );

        /* fills in array arlog with log of ionization fractions
         * N is set to highest stage with non-zero abundance
         * N set to 0 if element turned off
         *
         * first call MeanZero to sero out save arrays
         * next call MeanInc in zone calc to enter ionziation fractions or temperature
         * finally this routine computes actual means
         * */
        if( !dense.lgElmtOn[nelem] )
        {
                /* no ionization for this element */
                /* >>chng 04 jun 27, upper limit had been <nelem+1, had missed fully
                 * stipped ion */
                for( ion=0; ion < limit; ion++ )
                {
                        arlog[ion] = -30.f;
                }
                *n = 0;
                return;
        }

        /* set high ion stages, with zero abundances, to -30,
         * limit is on c scale, such that ion<=limit */
        *n = limit;
        while( *n > 0 && mean.xIonMeans[0][nelem][*n-1] <= 0. )
        {
                arlog[*n-1] = -30.f;
                --*n;
        }

        if( chType=='i' && lgDensity)
        {
                /* mean ionization with density*/
                for( ion=0; ion < *n; ion++ )
                {
                        if( mean.xIonEdenMeans[0][nelem][ion] > 0. )
                        {
                                abund_radmean = mean.xIonEdenMeans[0][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean / mean.xIonEdenMeansNorm[0][nelem]));
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else if( chType=='i' )
        {
                /* mean ionization no density*/
                for( ion=0; ion < *n; ion++ )
                {
                        if( mean.xIonMeans[0][nelem][ion] > 0. )
                        {
                                abund_radmean = mean.xIonMeans[0][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/mean.xIonMeansNorm[0][nelem]));
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else if( chType=='t' && lgDensity )
        {
                /* mean temperature wrt elec density */
                for( ion=0; ion < *n; ion++ )
                {
                        normalize = mean.TempEdenMeansNorm[0][nelem][ion];
                        if( normalize > SMALLFLOAT )
                        {
                                abund_radmean = mean.TempEdenMeans[0][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/normalize));
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else if( chType=='t'  )
        {
                /* mean temperature without elec density*/
                for( ion=0; ion < *n; ion++ )
                {
                        normalize = mean.TempMeansNorm[0][nelem][ion];
                        if( normalize > SMALLFLOAT )
                        {
                                abund_radmean = mean.TempMeans[0][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_radmean/normalize));
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else
        {
                fprintf(ioQQQ," MeanIonRadius called with insane job \n");
        }
        return;
}

/*MeanIonVolume do volume mean of ionization fractions over volume of any element */
void MeanIonVolume(
        /* either 'i' or 't' for ionization or temperature */
        char chType,
        /* atomic number on c scale */
        long int nelem, 
        /* this will say how many of arlog have non-zero values */
        long int *n, 
        /* array of values, log both cases */
        realnum arlog[],
        /* true, include electron density, false do not*/
        bool lgDensity )
{
        long int ion;
        double abund_volmean, 
          normalize;
        long int limit;

        DEBUG_ENTRY( "MeanIonVolume()" );

        ASSERT( chType=='i' || chType=='t' );

        /* >>chng 04 jun 27, add H2 to top of H ladder,
         * in this routine nelem is on physical scale, so 1 for H*/
        limit = MAX2( 3, nelem+2 );

        /* 
         * fills in array arlog with log of ionization fractions
         * N is set to highest stage with non-zero abundance
         * N set to 0 if element turned off
         *
         * first call MeanZero to sero out save arrays
         * next call MeanInc in zone calc to enter ionziation fractions
         * finally this routine computes actual means
         */
        if( !dense.lgElmtOn[nelem] )
        {
                /* no ionization for this element */
                for( ion=0; ion <= limit; ion++ )
                {
                        arlog[ion] = -30.f;
                }
                *n = 0;
                return;
        }

        /* this is number of stages of ionization */
        *n = limit;
        /* fill in higher stages with zero if non-existent, 
         * also decrement n, the number with non-zero abundances */
        while( *n > 0 && mean.xIonMeans[1][nelem][*n-1] <= 0. )
        {
                arlog[*n-1] = -30.f;
                --*n;
        }
        /* n is now the limit to the number with positive abundances */

        if( chType=='i' && lgDensity)
        {
                /* mean ionization with electron density*/
                /* this is denominator for forming a mean */
                /* return log of means */
                for( ion=0; ion < *n; ion++ )
                {
                        if( mean.xIonEdenMeans[1][nelem][ion] > 0. )
                        {
                                abund_volmean = mean.xIonEdenMeans[1][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean) / mean.xIonEdenMeansNorm[1][nelem]);
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else if( chType=='i' )
        {
                /* mean ionization */
                /* this is denominator for forming a mean */
                /* return log of means */
                for( ion=0; ion < *n; ion++ )
                {
                        if( mean.xIonMeans[1][nelem][ion] > 0. )
                        {
                                abund_volmean = mean.xIonMeans[1][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean) / mean.xIonMeansNorm[1][nelem]);
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else if( chType=='t' && lgDensity )
        {
                /* mean temperature with density*/
                /* this is denominator for forming a mean */
                /* return log of means */
                for( ion=0; ion < *n; ion++ )
                {
                        normalize = mean.TempEdenMeansNorm[1][nelem][ion];
                        if( normalize > SMALLFLOAT )
                        {
                                abund_volmean = mean.TempEdenMeans[1][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean)/normalize);
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else if( chType=='t' )
        {
                /* mean temperature with NO density*/
                /* this is denominator for forming a mean */
                /* return log of means */
                for( ion=0; ion < *n; ion++ )
                {
                        normalize = mean.TempMeansNorm[1][nelem][ion];
                        if( normalize > SMALLFLOAT )
                        {
                                abund_volmean = mean.TempMeans[1][nelem][ion];
                                arlog[ion] = (realnum)log10(MAX2(1e-30,abund_volmean)/normalize);
                        }
                        else
                        {
                                arlog[ion] = -30.f;
                        }
                }
        }
        else
        {
                fprintf(ioQQQ," MeanIonVolume called with insane job\n");
        }
        return;
}

/*aver compute average of various quantities over the computed geometry
 * called by startenditer to initialize, radinc to increment, and prtfinal for final results */
void aver(
        /* 4 char + null string giving job, prin, zero, zone, doit*/
        const char *chWhat, 
        /* the quantity to average, like the temperature */
        double quan, 
        /* what we weight against, like the o++ abundance */
        double weight, 
        /* 10 char + null string giving title for this average */
        const char *chLabl)
{
        /* NAVER is the limit to the number than can be averaged */
#       define  NAVER   20
        static char chLabavr[NAVER][11];
        long int i, 
          ioff, 
          j;
        static long int n;
        double raver[NAVER]={0.}, 
          vaver[NAVER]={0.};
        static double aversv[4][NAVER];

        DEBUG_ENTRY( "aver()" );

        if( strcmp(chWhat,"zero") == 0 )
        {
                /* chWhat='zero' - zero out the save array */
                for( i=0; i < NAVER; i++ )
                {
                        for( j=0; j < 4; j++ )
                        {
                                aversv[j][i] = 0.;
                        }
                }
                n = 0;
        }
        else if( strcmp(chWhat,"zone") == 0 )
        {
                n = 0;
        }
        else if( strcmp(chWhat,"doit") == 0 )
        {
                /* chWhat='doit' - enter values to average */
                if( n >= NAVER )
                {
                        fprintf( ioQQQ, " Too many values entered into AVER, increase NAVER\n" );
                        cdEXIT(EXIT_FAILURE);
                }
                aversv[0][n] += weight*quan*radius.drad_x_fillfac;
                aversv[1][n] += weight*radius.drad_x_fillfac;
                aversv[2][n] += weight*quan*radius.dVeff;
                aversv[3][n] += weight*radius.dVeff;

                /* this is the label for this particular average, like " TeNe "*/
                strcpy( chLabavr[n], chLabl );
                n += 1;
        }
        else if( strcmp(chWhat,"prin") == 0 )
        {
                /* set things up to get answers */
                ioff = n*10/2 + 1;

                /* space out to center the label on the numbers */
                fprintf( ioQQQ,"\n");
                for( i=0; i < ioff; i++ )
                {
                        /*chTitAvr[i] = ' ';*/
                        fprintf( ioQQQ, " " );
                }

                /* now print header with cr */
                fprintf( ioQQQ, "Averaged Quantities\n " );

                fprintf( ioQQQ, "        " );
                for( i=0; i < n; i++ )
                {
                        fprintf( ioQQQ, "%10.10s", chLabavr[i] );
                }
                fprintf( ioQQQ, " \n" );

                for( i=0; i < n; i++ )
                {
                        if( aversv[1][i] > 0. )
                        {
                                raver[i] = aversv[0][i]/aversv[1][i];
                        }
                        else
                        {
                                raver[i] = 0.;
                        }
                        if( aversv[3][i] > 0. )
                        {
                                vaver[i] = aversv[2][i]/aversv[3][i];
                        }
                        else
                        {
                                vaver[i] = 0.;
                        }
                }

                fprintf( ioQQQ, " Radius:" );
                for( i=0; i < n; i++ )
                {
                        /*fprintf( ioQQQ, "%11.3e", raver[i] );*/
                        fprintf( ioQQQ, " ");
                        fprintf(ioQQQ,PrintEfmt("%9.2e", raver[i] ) );
                }
                fprintf( ioQQQ, "\n" );

                /* only print vol aver is lgSphere is set */
                if( geometry.lgSphere )
                {
                        fprintf( ioQQQ, " Volume:" );
                        for( i=0; i < n; i++ )
                        {
                                /*fprintf( ioQQQ, "%11.3e", vaver[i] );*/
                                fprintf( ioQQQ, " ");
                                fprintf(ioQQQ,PrintEfmt("%9.2e", vaver[i] ) );
                        }
                        fprintf( ioQQQ, "\n" );
                }
        }
        else
        {
                fprintf( ioQQQ, " AVER does not understand chWhat=%4.4s\n", 
                  chWhat );
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

---