/*******************************************************************

Program: autocor.c   Eric R. Weeks  3/10/97

link to this program: http://www.physics.emory.edu/~weeks/software/autocor.html
email: weeks@physics.emory.edu

This program is public domain, but you are not allowed to remove
the web or email information.  No guarantees are made, but I think
the program is bug-free.

-------autocor.c-----------
03-10-97: started from cor.c
		add -t option

 *******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX 500100

char *prgname;

main(argc,argv)
int argc;
char **argv;
{
	int r,t;
	int c,numpts;
	FILE *fp, *fopen();
	extern int optind;
	extern char *optarg;
	double sum[2],avg[2],sumsq[2],std[2],max[2],min[2];
	double dat[MAX][2];
	double sqrt();
	double temp,cor;
	int tau;
	float deltat=-1;

	tau = 100;
    while ((c = getopt(argc, argv, "ht:d:")) != EOF)
        switch (c)
        {
        case 'h':  
fprintf(stderr,"input should be one column of data\n\n");
fprintf(stderr,"Usage: cor [<] datafile\n");
fprintf(stderr,"\nOptions:\n");
fprintf(stderr," -h   : this help option\n");
fprintf(stderr," -d # : set maximum delay for autocorrelation [%d]\n",tau);
fprintf(stderr," -t # : set delta-t between each point [1.0]\n");
               exit(1);
               break;
	case 'd': tau = atoi(optarg);
			break;
	case 't': deltat = atof(optarg);
			break;
        }

    argc -= (optind-1) ; argv += (optind-1) ;
    fp = (argc > 1) ? fopen(*++argv, "r") : stdin;

	numpts = 0;
	for (r=0;r<2;r++)  {
		sum[r] = 0.0;
		sumsq[r] = 0.0;
		max[r] = -9.999e99;
		min[r] = +9.999e99;
	}

	while (fscanf(fp,"%lf",&dat[numpts][0]) == 1)  {
		for (r=0;r<1;r++)  {
			sum[r] = sum[r] + dat[numpts][r];
			sumsq[r] = sumsq[r] + dat[numpts][r]*dat[numpts][r];
			max[r] = (max[r] > dat[numpts][r]) ? max[r] : dat[numpts][r];
			min[r] = (min[r] < dat[numpts][r]) ? min[r] : dat[numpts][r];
		}
		++numpts;
		if (numpts >= MAX)  {
			fprintf(stderr,"too many data points, limit is %d\n",MAX);
			exit(1);
		}
	}

	fprintf(stderr,"   avg        std        max         min        number\n");
	for (r=0;r<1;r++)  {
		avg[r] = sum[r]/numpts;
		std[r] = sqrt((sumsq[r]/numpts)-avg[r]*avg[r]);
		fprintf(stderr,"%10g %10g %10g",avg[r],std[r],max[r]);
		fprintf(stderr," %10g %10d\n",min[r],numpts);
		std[r] = 1.0/std[r];
	}
	cor = 0.0;
	for (t=0;t<numpts;t++)
		dat[t][0] = (dat[t][0]-avg[0])*std[0];

	for (r=0;r<tau;r++)  {
		cor = 0.0;
		for (t=0;t<(numpts-r);t++)
			cor += dat[t][0]*dat[t+r][0];
		cor = cor / ((double)(numpts-r));
		if (deltat<0.0)
			printf("%d %g\n",r,cor);
		else
			printf("%f %g\n",((float)r)*deltat,cor);
	}
	exit(0);
}