/*gett2o3 analyze computed [OIII] spectrum to get t^2 */
#include "cddefines.h"
#include "struc.h"
#include "zonecnt.h"
#include "gett2o3.h"

void gett2o3(float *tsqr)
{
	long int i, 
	  limit;
	double tmean;
	double a, 
	  as, 
	  b;

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

	/* get T, t^2 */
	a = 0.;
	b = 0.;
	limit = MIN2(NZLIM,ZoneCnt.nzone);
	for( i=0; i < limit; i++ )
	{
		as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
		  (double)(struc.ednstr[i]);
		a += (double)(struc.testr[i])*as;

		/* B is used twice below */
		b += as;
	}

	if( b <= 0. )
	{
		*tsqr = 0.;
	}

	else
	{
		/* following is H+ weighted mean temp over vol */
		tmean = a/b;
		a = 0.;
		limit = MIN2(NZLIM,ZoneCnt.nzone);
		for( i=0; i < limit; i++ )
		{
			as = (double)(struc.volstr[i])*(double)(struc.o3str[i])*
			  struc.ednstr[i];
			a += (POW2((double)(struc.testr[i]-tmean)))*as;
		}
		*tsqr = (float)(a/(b*(POW2(tmean))));
	}


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