/*MeanInc increment mean ionization fractions and temperatures over computed structure, in tauinc */
/*MeanZero zero mean of ionization fractions array */
/*RadMean derive mean ionization fractions or temperature over ravius for any element */
/*VolMean do volume mean ionization fractions or temperature over volume for any element */
#include "cddefines.h"
#include "ionfracs.h"
#include "filfac.h"
#include "radius.h"
#include "nomole.h"
#include "hmi.h"
#include "phycon.h"
#include "mean.h"

void MeanInc(void)
{
	long int i, 
	  nelem;
	double dc, 
	  dcv, 
	  sumrad, 
	  sumvol;

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

	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* this is mean ionization */
		sumrad = 0.;
		sumvol = 0.;
		for( i=0; i < (nelem + 2); i++ )
		{
			/* xIonMeans[0][i+1][n] is radial integration */
			dc = xIonFracs[i+1][nelem]*radius.drad*filfac.FillFac;
			IonMeans.xIonMeans[0][i+1][nelem] += dc;

			/* xIonMeans[1][i+1][n] is vol integration */
			dcv = dc*radius.dVeff;
			IonMeans.xIonMeans[1][i+1][nelem] += dcv;

			sumrad += dc;
			sumvol += dcv;
		}
		IonMeans.xIonMeans[0][0][nelem] += sumrad;
		IonMeans.xIonMeans[1][0][nelem] += sumvol;

		/* this is mean temperature */
		sumrad = 0.;
		sumvol = 0.;
		for( i=0; i < (nelem + 2); i++ )
		{
			/* TempMeans[0][i+1][n] is radial integration */
			dc = xIonFracs[i+1][nelem]*radius.drad*filfac.FillFac;
			IonMeans.TempMeans[0][i+1][nelem] += dc*phycon.te;
			IonMeans.TempMeans[0+2][i+1][nelem] += dc;

			/* TempMeans[1][i+1][n] is vol integration */
			dcv = dc*radius.dVeff;
			IonMeans.TempMeans[1][i+1][nelem] += dcv*phycon.te;
			IonMeans.TempMeans[1+2][i+1][nelem] += dcv;
		}
	}
	
	/* used to get harmonic mean temperature with respect to H over radius,
	 * for comparison with 21cm temperature */
	/* >>chng 99 jul 7, no h2 option, include h2 as it were two H's */
	if( nomole.lgNoH2Mole )
	{
		dc = (xIonFracs[1][0]+2.*hmi.htwo)*radius.drad*filfac.FillFac;
	}
	else
	{
		dc = xIonFracs[1][0]*radius.drad*filfac.FillFac;
	}
	IonMeans.HarMeanTemp[0] += dc/phycon.te;
	IonMeans.HarMeanTemp[1] += dc;

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

/*MeanZero zero mean of ionization fractions array */
void MeanZero(void)
{
	long int i, 
	  j, 
	  nelem;

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

	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		for( i=0; i <= (nelem + 2); i++ )
		{
			for( j=0; j < 2; j++ )
			{
				/* these are used to save info on temperature and ionization means */
				IonMeans.xIonMeans[j][i][nelem] = 0.;
				/* double here since [2] and [3] have norm for average */
				IonMeans.TempMeans[j][i][nelem] = 0.;
				IonMeans.TempMeans[j+2][i][nelem] = 0.;
			}
		}
	}

	/* used to get harmonic mean temperature with respect to H over radius,
	 * for comparison with 21cm temperature */
	IonMeans.HarMeanTemp[0] = 0.;
	IonMeans.HarMeanTemp[1] = 0.;

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

/*RadMean derive mean ionization fractions over ravius for some element */
void RadMean(
	/* either 'i' or 't' for ionization or temperature */
	char chType,
	/* atomic number on physical, no c, scale */
	long int nelem, 
	/* this will say how many of arlog have non-zero values */
	long int *n, 
	/* array of values, log bth cases */
	float arlog[])
{
	long int i;
	double abund, 
	  normalize;

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

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

	/* 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( IonMeans.xIonMeans[0][0][nelem-1] <= 0. )
	{
		/* no ionization for this element */
		for( i=0; i < (nelem + 1); i++ )
		{
			arlog[i] = -30.f;
		}
		*n = 0;
		
#		ifdef DEBUG_FUN
		fputs( " <->RadMean()\n", debug_fp );
#		endif
		return;
	}

	*n = nelem + 1;
	while( *n > 1 && IonMeans.xIonMeans[0][*n][nelem-1] <= 0. )
	{
		arlog[*n-1] = -30.f;
		*n -= 1;
	}

	if( chType=='i' )
	{
		/* mean ionization */
		normalize = IonMeans.xIonMeans[0][0][nelem-1];
		for( i=0; i < *n; i++ )
		{
			if( IonMeans.xIonMeans[0][i+1][nelem-1] > 0. )
			{
				abund = IonMeans.xIonMeans[0][i+1][nelem-1];
				arlog[i] = (float)log10(MAX2(1e-30,abund/normalize));
			}
			else
			{
				arlog[i] = -30.f;
			}
		}
	}
	else
	{
		/* mean temperature */
		for( i=0; i < *n; i++ )
		{
			normalize = IonMeans.TempMeans[0+2][i+1][nelem-1];
			if( normalize > SMALLFLOAT )
			{
				abund = IonMeans.TempMeans[0][i+1][nelem-1];
				arlog[i] = (float)log10(MAX2(1e-30,abund/normalize));
			}
			else
			{
				arlog[i] = -30.f;
			}
		}
	}


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

/*VolMean do volume mean of ionization fractions over volumn of any element */
void VolMean(
	/* either 'i' or 't' for ionization or temperature */
	char chType,
	/* atomic number on physical, no 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[])
{
	long int i;
	double abund, 
	  normalize;

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

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

	/* 
	 * 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( IonMeans.xIonMeans[1][0][nelem-1] <= 0. )
	{
		/* no ionization for this element */
		for( i=0; i < (nelem + 1); i++ )
		{
			arlog[i] = -30.f;
		}
		*n = 0;
		
#		ifdef DEBUG_FUN
		fputs( " <->VolMean()\n", debug_fp );
#		endif
		return;
	}

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

	if( chType=='i' )
	{
		/* mean ionization */
		/* this is denominator for forming a mean */
		normalize = IonMeans.xIonMeans[1][0][nelem-1];
		/* return log of means */
		for( i=0; i < *n; i++ )
		{
			if( IonMeans.xIonMeans[1][i+1][nelem-1] > 0. )
			{
				abund = IonMeans.xIonMeans[1][i+1][nelem-1];
				arlog[i] = (float)log10(MAX2(1e-30,abund)/normalize);
			}
			else
			{
				arlog[i] = -30.f;
			}
		}
	}
	else
	{
		/* mean temperature */
		/* this is denominator for forming a mean */
		/* return log of means */
		for( i=0; i < *n; i++ )
		{
			normalize = IonMeans.TempMeans[1+2][i+1][nelem-1];
			if( normalize > SMALLFLOAT )
			{
				abund = IonMeans.TempMeans[1][i+1][nelem-1];
				arlog[i] = (float)log10(MAX2(1e-30,abund)/normalize);
			}
			else
			{
				arlog[i] = -30.f;
			}
		}
	}

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