.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include "sacdata.h"
 
long int read_sac(char *fname, SAC_DATA *sdata, float *beg, float *dt);
 
SAC_DATA ew, ns;
 
float XCOR(long int npts, float *ew.fdat[], float *ns.fdat[]);
 
int main(int argc, char *argv[])
{
	int i, j, k;
	float bege, begn, dt;
	float coeff;
 
	coeff = 0;
 
	if( argc < 3 )
	{
		fprintf(stderr,"Usage: %s e.sac n.sac\n",argv[0]);
		return -1;
	}
 
	read_sac(argv[1],&ew,&bege,&dt);
	read_sac(argv[2],&ns,&begn,&dt);
 
	coeff = XCOR(ew.npts, &ew.fdat[], &ns.fdat[]);
 
	fprintf(stdout, "  coe : %.2f\n", coeff);
 
	return 0;
}
 
float XCOR(long int npts, float *ew.fdat[], float *ns.fdat[])
{
	int i, j, k;
	float bege, begn, dt;
	float ewsum, nssum, ewavg, nsavg;
	float ewsq, nssq, ewd, nsd;
	float ewns, r;
	float npts;
 
	ewsum = 0.,   nssum = 0.,   ewavg = 0.,   nsavg = 0.;
	ewsq = 0.,   nssq = 0.,   ewd = 0.,   nsd = 0.;
	ewns = 0.,   r = 0.,   npts = 0.;
 
	/////////////// Phase info. ///////////////
 
	npts = ew.npts;
 
	for (i=0; i<npts; i++)
	{
		if (ew.fdat[i] > 0)
		{
		    	ew.fdat[i] = 1.;
		}
			else if (ew.fdat[i] < 0)
		{
			ew.fdat[i] = -1.;
		}
			else
		{
			ew.fdat[i] = 0;
		}
	}
 
	for (i=0; i<npts; i++)
	{
		if (ns.fdat[i] > 0)
		{
			ns.fdat[i] = 1.;
		}
		else if (ns.fdat[i] < 0)
		{
			ns.fdat[i] = -1.;
		}
		else
		{
			ns.fdat[i] = 0;
		}
	}
 
 
	////////// Average //////////
 
	for (i=0; i<npts; i++)
	{
		ewsum += ew.fdat[i];
	}
 
	for (i=0; i<npts; i++)
	{
		nssum += ns.fdat[i];
	}
 
	ewavg = ewsum / npts;
	nsavg = nssum / npts;
 
 
	////////// Standard Deviation ///////////
 
	for (i=0; i<npts; i++)
	{
		ewsq += (ew.fdat[i] - ewavg) * (ew.fdat[i] - ewavg);
	}
 
	for (i=0; i<npts; i++)
	{
		nssq += (ns.fdat[i] - nsavg) * (ns.fdat[i] - nsavg);
	}
 
	ewd = sqrt (ewsq / npts);
	nsd = sqrt (nssq / npts);
 
 
	///////// crosscorrelation coefficient /////////
 
	for (i=0; i<npts; i++)
	{
		ewns += ((ew.fdat[i] - ewavg) / ewd) * ((ns.fdat[i] - nsavg) / nsd);
	}
 
	r = ewns / npts;
 
	return r;
}