/* This file is part of Cloudy and is copyright (C) 1978-2004 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 */
#include "cddefines.h"
#include "radius.h"
#include "dense.h"
#include "hyperfine.h"
#include "hmi.h"
#include "phycon.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 = 3;
			/* 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.TempMeans[0+2][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.TempEdenMeans[0+2][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.TempMeans[1+2][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.TempEdenMeans[1+2][nelem][ion] += dcv;
		}
		for( ion=1; ion <= (nelem + 2); ion++ )
		{
			/* TempMeans[0][ion][n] is radial integration */
			dc = dense.xIonDense[nelem][ion-1]*radius.drad_x_fillfac;
			mean.TempMeans[0][nelem][ion] += dc*phycon.te;
			mean.TempMeans[0+2][nelem][ion] += dc;

			/* TempMeans[0][ion][n] is radial integration PLUS electron density*/
			dc = dense.xIonDense[nelem][ion-1]*radius.drad_x_fillfac*dense.eden;
			mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
			mean.TempEdenMeans[0+2][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-1]*radius.dVeff;
			mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
			mean.TempMeans[1+2][nelem][ion] += dcv;

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

	/* =============================================================
	 *
	 * these are some special quanties for the mean 
	 */
	
	/* 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 */
	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.Hmolec[ipMH2g]*radius.drad_x_fillfac;*/
	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.Hmolec[ipMH2g]*radius.dVeff;*/
	dcv = hmi.H2_total*radius.dVeff;
	mean.H2MeanTempVolume[0] += dc*phycon.te;
	mean.H2MeanTempVolume[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();
		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();
			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();
			}
		}
	}

	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++ )
			{
				/* 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.TempMeans[j+2][nelem][ion] = 0.;
				mean.TempEdenMeans[j][nelem][ion] = 0.;
				mean.TempEdenMeans[j+2][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.;

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

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

