#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>

#define MAX 2000000


void files(void);

FILE *pFile[4];
char filename[4][FILENAME_MAX];

long i,j,last;

double t,a[MAX],peak[50000];

double ave,ms,s,std;
double lower,upper;


void main()
{
	printf("\n\n\n\n peak_dist.cpp, ver 1.0  July 15, 2003");
	printf("\n\n by Tom Irvine ");
	printf("\n Email:  tomirvine@aol.com ");


	files();


	i=0;
	j=0;

	while(1)
	{
        if( fscanf(pFile[0],"%lf %lf \n",&t,&a[i]) <= 0 ){break;}

        ave+=a[i];

		ms+=pow( a[i],2.);


		i++;

		if(i==MAX){break;}
	}
	last = i;

	ave/=last;

	std = sqrt( (ms/last) - pow(ave,2.) );


	for(i=0; i<last; i++)
	{

		a[i]-=ave;
	}


    double lower,upper;

    lower=0.;
	upper=0.;

	for(i=0; i < last; i++)
	{
		if(a[i] > upper ){ upper = a[i]; }
		if(a[i] < lower ){ lower = a[i]; }
	}

	printf("\n Total samples = %ld \n",last);

	printf("\n minimum = %12.4g   maximum = %12.4g\n",lower,upper);


	printf("\n Standard Deviation = %12.4g \n",std);

    if(fabs(upper)< fabs(lower)){upper=fabs(lower);}

	printf("\n Highest Sigma = %12.4g sigma\n",upper/std);


    for(i=1; i< (last-1); i++)
	{

		if( a[i] >= a[i-1] && a[i] >= a[i+1] && a[i] > 0)
		{

			   peak[j]=fabs(a[i]);
               j++;

		       if(j==50000){break;}

		}
		if( a[i] <= a[i-1] && a[i] <= a[i+1] && a[i] < 0)
		{
			   peak[j]=fabs(a[i]);
               j++;

		       if(j==50000){break;}
		}

	}
    last =j;

    printf("\n Total peaks = %ld \n",last);
 
	upper=0.;

	for(i=0; i < last; i++)
	{
		peak[i]/=std;
	}


//	printf("\n maximum = %12.4g\n\n",upper);


	double m[60];
	long n[60];


	for(i=0; i < 60; i++)
	{
		n[i]= 0;
	}

	long nb=54;

	double bw = 0.25;

	for(i=0; i <= nb; i++)
	{
		m[i]= i*bw;
	}



	for(i=0; i< last; i++)
	{

		for( j = 0; j<=nb; j++)
		{
			if( peak[i] > m[j] && peak[i] <= m[j+1] ){ n[j]++; }
		}
	}
	 
	s=0.;

	fprintf(pFile[1],"%12.4g %ld\n",s,long(s));

	for( j = 0; j< nb; j++)
	{
		 s= (m[j]+m[j+1])/2.;

         if( s > 6.5 ){break;}

		 fprintf(pFile[1],"%12.4g %12.4g\n",s,(double(n[j])/double(last))/bw);
	}


    long n3=0;
    long n4=0;

	for( j=0; j<last; j++)
	{

		if( peak[j]>3){n3++;}
		if( peak[j]>4){n4++;}
	}

	printf("\n number of peaks > 3 sigma = %ld \n",n3);
	printf("\n number of peaks > 4 sigma = %ld \n",n4);

   printf("\n\n The output file is:  peak_dist.out \n");
   printf("\n The file has two columns:  nsigma  & probability density \n\n");

   printf("\n Press any key to exit.\n");
   getch();


}
void files(void)
{

    printf( "\n The base input file must contain two columns: \n");
    printf( " time(sec)  accel(G)    \n");

	printf("\n The program will remove the mean from the input data. \n");

    printf( " \n An error message will occur if the file does not exist.\n");
    printf( " Input filename \n");


	scanf("%s",filename[0]);
	pFile[0]=fopen(filename[0], "rb");

 
    i=0;


    if(pFile[i] == NULL )
    {
         printf("\n Failed to open file: %s \n", filename[i]);
         fclose(pFile[i]);
    }
    else
    {
         printf("\n File: %s opened. \n", filename[i]);
    }

	strcpy(filename[1],"peak_dist.out");
	pFile[1]=fopen(filename[1], "w");

}