#include <stdio.h>
#include <math.h>
#include "stats.h"
#include "dynavars.h"

void stats(xn,sr) /* running min,max,mean and variance */
double xn;
struct statrec *sr;
/*
	To use: Initialize n to zero, then pass values thru xn. When done, min, max,
		mean and variance will all have correct values (nb. standard deviation
		equals sqrt(variance)). All values must be preserved between calls.
						   ----
		The formulas for running mean and variance were derived from formulas
		for mean and variance given on page 597 of "The VNR Concise Encyclopedia
		of Mathematics," Van Nostrand Reinhold: Toronto, 1975. Note that the
		quoted mean and variance for the data on page 596 is wrong. (awk got
		mean = 42.11 and variance = 8.60876 not 42.14 and 6.84.)
		
		Written by Allen Supynuk, University of Alberta, October 1990
*/
{
	double xnsq = xn * xn;

	(sr->n)++;
	if( sr->n < 1 ) {
		fprintf(stderr,"*** Error: stat called with n < 1");
		exit(sr->n);
	}
	if( sr->n == 1 ) {
		sr->min      = xn;
		sr->max      = xn;
		sr->variance = 0.0;
		sr->sum      = xn;
		sr->sumsq    = xnsq;
		sr->mean     = xn;
	}
	else if( sr->n == 2 ) {
		if( xn < sr->min ) sr->min = xn;
		if( xn > sr->max ) sr->max = xn;
		sr->mean    += (xn - sr->mean) / sr->n;
		sr->variance = (   xn   - sr->mean) * (   xn   - sr->mean)
		             + (sr->sum - sr->mean) * (sr->sum - sr->mean);
		sr->sum     += xn;
		sr->sumsq   += xnsq;
	}
	else {
		if( xn < sr->min ) sr->min = xn;
		if( xn > sr->max ) sr->max = xn;
		sr->variance += xnsq / sr->n
			      - 2.0 * xn * sr->mean / sr->n
			      - sr->sumsq / ((sr->n - 1.0) * (sr->n - 2.0))
			      + 2.0 * sr->sum * sr->mean / (sr->n * (sr->n - 2.0));
		sr->sum      += xn;
		sr->sumsq    += xnsq;
		sr->mean     += (xn - sr->mean) / sr->n;
	}
}

void prsr(sr)
struct statrec sr;
{
	printf("n = %7d, min = %le, max = %le, mean = %le, std dev = %le, sum = %le, sumsq = %le, variance = %le\n",
		   sr.n,    sr.min,    sr.max,    sr.mean,sqrt(sr.variance), sr.sum,    sr.sumsq,    sr.variance    );
}

void print_stats()
{
	printf("     vHR stats:\n"); prsr(vHR_stats);
	printf("      aH stats:\n"); prsr(aH_stats);
	printf("   omega stats:\n"); prsr(omega_stats);
	printf("omegadot stats:\n"); prsr(omegadot_stats);
	printf("     FER stats:\n"); prsr(FER_stats);
	printf("    GEHR stats:\n"); prsr(GEHR_stats);
	printf("      gH stats:\n"); prsr(gH_stats);
	printf("     ROT stats:\n"); prsr(ROT_stats);
}

void init_stats()
{
	vHR_stats.n = 0;
	aH_stats.n = 0;
	omega_stats.n = 0;
	omegadot_stats.n = 0;
	FER_stats.n = 0;
	GEHR_stats.n = 0;
	gH_stats.n = 0;
	ROT_stats.n = 0;
}

void stats_vec(array_vector, sr, fr) /* collect stats on an array of vectors */
double array_vector[MAX][3];
struct statrec *sr;
enum frame fr;
{
	int link,i;
	double tempv[3];

	for(link = 1; link <= 3; link++) {
		if( fr == InertialFrame ) {
			mxv(RIT[link],array_vector[link],tempv);
			for(i = 0; i < 3; i++) {
				stats(tempv[i],sr);
		}	}
		else {
			for(i = 0; i < 3; i++) {
				stats(array_vector[link][i],sr);
}	}	}	}

void stats_mat(array_matrix, sr) /* collect stats on an array of matrices */
double array_matrix[MAX][3][3];
struct statrec *sr;
{
	int link,i,j;

	for(link = 1; link <= 3; link++) {
		for(i = 0; i < 3; i++) {
			for(j = 0; j < 3; j++) {
				stats(array_matrix[link][i][j],sr);
}	}	}	}
