/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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 */

#	ifdef DEBUG_FUN
	fputs( "<+>MeanInc()\n", debug_fp );
#	endif

	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(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;

#	ifdef DEBUG_FUN
	fputs( " <->MeanInc()\n", debug_fp );
#	endif
	return;
}

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

#	ifdef DEBUG_FUN
	fputs( "<+>MeanZero()\n", debug_fp );
#	endif

	/* 
	 * 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] */
		if( (mean.xIonMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.xIonEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.xIonMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.xIonEdenMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.TempMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.TempEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.TempMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		if( (mean.TempEdenMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) ) )==NULL )
			BadMalloc();
		for( j=0; j<2; ++j )
		{
			if( (mean.xIonMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.xIonEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.xIonMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.xIonEdenMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.TempMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.TempEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.TempMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			if( (mean.TempEdenMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) ) )==NULL )
				BadMalloc();
			for( nelem=0; nelem<LIMELM; ++nelem )
			{
				limit = MAX2(3,nelem+2);
				if( (mean.xIonMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ) )==NULL )
					BadMalloc();
				if( (mean.xIonEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ) )==NULL )
					BadMalloc();
				if( (mean.TempMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ) )==NULL )
					BadMalloc();
				if( (mean.TempEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ) )==NULL )
					BadMalloc();
				if( (mean.TempMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ) )==NULL )
					BadMalloc();
				if( (mean.TempEdenMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) ) )==NULL )
					BadMalloc();
			}
		}
	}

	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.;

#	ifdef DEBUG_FUN
	fputs( " <->MeanZero()\n", debug_fp );
#	endif
	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  */
	float arlog[],
	/* TRUE, include electron density, FALSE do not*/
	int lgDensity )
{
	long int ion,
		limit;
	double abund_radmean, 
	  normalize;

#	ifdef DEBUG_FUN
	fputs( "<+>MeanIonRadius()\n", debug_fp );
#	endif

	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;
		
#		ifdef DEBUG_FUN
		fputs( " <->MeanIonRadius()\n", debug_fp );
#		endif
		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] = (float)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] = (float)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] = (float)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] = (float)log10(MAX2(1e-30,abund_radmean/normalize));
			}
			else
			{
				arlog[ion] = -30.f;
			}
		}
	}
	else
	{
		fprintf(ioQQQ," MeanIonRadius called with insane job \n");
	}

#	ifdef DEBUG_FUN
	fputs( " <->MeanIonRadius()\n", debug_fp );
#	endif
	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 */
	float arlog[],
	/* TRUE, include electron density, FALSE do not*/
	int lgDensity )
{
	long int ion;
	double abund_volmean, 
	  normalize;
	long int limit;

#	ifdef DEBUG_FUN
	fputs( "<+>MeanIonVolume()\n", debug_fp );
#	endif

	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;
		
#		ifdef DEBUG_FUN
		fputs( " <->MeanIonVolume()\n", debug_fp );
#		endif
		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] = (float)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] = (float)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] = (float)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] = (float)log10(MAX2(1e-30,abund_volmean)/normalize);
			}
			else
			{
				arlog[ion] = -30.f;
			}
		}
	}
	else
	{
		fprintf(ioQQQ," MeanIonVolume called with insane job\n");
	}

#	ifdef DEBUG_FUN
	fputs( " <->MeanIonVolume()\n", debug_fp );
#	endif
	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];

#	ifdef DEBUG_FUN
	fputs( "<+>aver()\n", debug_fp );
#	endif

	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" );
			puts( "[Stop in aver]" );
			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();
		puts( "[Stop in aver]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->aver()\n", debug_fp );
#	endif
	return;
}

