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

#include "filter.h"
#include "dates.h"


#define NINT(t)         ((int) ((t<0) ? t-.5 : t+.5))

int parsearg();
int read_data();
int check_export();
void set_export();
void compute_curves();
void print_coef();
void print_stats();
void print_ampstats();
void print_monthstats();
void get_amplitudes();
static void print_averages ();

double	xp[MAX_OBS]; 			/* array holding dates of input data */
double  yp[MAX_OBS];			/* array holding mixing ratio of input data */
double	sig[MAX_OBS];			/* array holding uncertainties of points in yp array */
int	np;				/* number of points in xp, yp arrays */
double	xfilt[MAX_OBS];			/* dates of filter results, equally spaced */
double  smooth[MAX_OBS];		/* short term filter results, equally spaced  */
double  *smooth2;			/* short term filter results, specified at data points  */
double  trend[MAX_OBS];			/* long term filter results, equally spaced */
double  *trend2;			/* long term filter results, specified at data points */
double  deriv[MAX_OBS];			/* growth rate results, equally spaced */
double  *deriv2;			/* growth rate results, specified at data points */
double	work[MAX_OBS];			/* work array */
int	nfilt;				/* number of points in xfilt, smooth, trend, deriv */
double	**covar; 			/* covariance matrix from function fit */


/***************************************************************/
static void usage() 
{
	printf ("Usage: ccgcrv [OPTION] [FILE]\n\n");
	printf (" Filtering Options: \n");
	printf ("   -npoly num		Number of polynomial terms in function.\n");
	printf ("   -nharm num		Number of harmonic terms in function.\n");
	printf ("   -interv num		Sampling interval in days.\n");
	printf ("   -short num		Short-term filter cutoff in days.\n");
	printf ("   -long num		Long-term filter cutoff in days.\n");
	
	printf ("\n");
	printf (" Output Options: \n");
	printf ("   -f file             Write output data to file instead of stdout.\n");
	printf ("   -equal              Output data at equal intervals.\n");
	printf ("   -cal                Output dates in calendar format.\n");
	printf ("   -hr                 Include hour in calendar format.\n");
	printf ("   -date yy:mm:dd      Output data starting at date.\n");
	printf ("   -user file          Output data based on user supplied dates in file.\n");
	printf ("   -orig               Print original data points.\n");
	printf ("   -func               Print function values.\n");
	printf ("   -poly               Print polynomial values.\n");
	printf ("   -smooth             Print smoothed data.\n");
	printf ("   -trend              Print long term trend values.\n");
	printf ("   -detrend            Print detrended values.\n");
	printf ("   -smcycle            Print smoothed, detrended annual cycle.\n");
	printf ("   -harm               Print values of annual harmonic functions.\n");
	printf ("   -res                Print residuals from the functions.\n");
	printf ("   -smres              Print smoothed residuals from the functions.\n");
	printf ("   -trres              Print long-term smoothed residuals from the function.\n");
	printf ("   -ressm              Print residuals from the smoothed curve.\n");
	printf ("   -gr                 Print growth rate values.\n");
	printf ("   -coef num1,num2     Print coefficients from index num1 to index num2.\n");
	printf ("   -stats              Print table of summary statistics for curve fit.\n");
	printf ("   -amp                Print table of statistics for annual amplitudes.\n");
		
}

/***************************************************************/
int main(argc, argv)
int argc;
char **argv;
{
	int n, r, hour, minute, second, do_coef = 0;
	int begcoef, endcoef;
	char *file = NULL;
	ExportData export = {0.0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	FilterData filter;
	int	equal=0, syear= -1, smonth= -1, sday= -1;
	int do_amp = 0, do_mm = 0, do_column = 0;
	int sample = 0, all=0, unique=0;
	char userfile[256], outfile1[256], outfile2[256];
	FILE *fp;

/* set some defaults */

	filter.shortterm = 80;
	filter.longterm = 667;
	filter.numpolyterms = 3;
	filter.numharmonics = 4;
	filter.timezero = -1.0;
	filter.sampleinterval = 0.0;
	covar = dmatrix(0, MAX_PARAM, 0, MAX_PARAM);

	n=1;
	while (n < argc) {
		if (strcmp (argv[n], "-npoly") == 0) {
			r = sscanf (argv[n+1], "%d", &filter.numpolyterms);
			if (r <= 0) {
				fprintf (stderr, "Cannot get -npoly argument:\n");
				exit (1);
			}
			if (filter.numpolyterms < 0 || filter.numpolyterms > 10) {
				fprintf (stderr, "Error in -npoly argument: value out of range (0-10) %s\n", argv[n+1]);
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-nharm") == 0) {
			r = sscanf (argv[n+1], "%d", &filter.numharmonics);
			if (r <= 0) {
				fprintf(stderr, "Cannot get -nharm argument:\n");
				exit (1);
			}
			if (filter.numharmonics < 0 || filter.numharmonics > 10) {
				fprintf (stderr, "Error in -nharm argument: value out of range (0-10) %s\n", argv[n+1]);
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-interv") == 0) {
			r = sscanf (argv[n+1], "%f", &filter.sampleinterval);
			if (r <= 0) {
				fprintf(stderr, "Cannot get -interv argument:\n");
				exit (1);
			}
			if (filter.sampleinterval < 0) {
				fprintf (stderr, "Error in -interv argument: %s\n", argv[n+1]);
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-short") == 0) {
			r = sscanf (argv[n+1], "%d", &filter.shortterm);
			if (r <= 0) {
				fprintf(stderr, "Cannot get -short argument:\n");
				exit (1);
			}
			if (filter.shortterm < 0) {
				fprintf (stderr, "Error in -short argument: %s\n", argv[n+1]);
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-long") == 0) {
			r = sscanf (argv[n+1], "%d", &filter.longterm);
			if (r <= 0) {
				fprintf(stderr, "Cannot get -long argument:\n");
				exit (1);
			}
			if (filter.longterm < 0) {
				fprintf (stderr, "Error in -long argument: %s\n", argv[n+1]);
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-timez") == 0) {
			r = sscanf (argv[n+1], "%f", &filter.timezero);
			if (r <= 0) {
				fprintf(stderr, "Cannot get -timez argument:\n");
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-equal") == 0) {
			equal = 1;
		} else if (strcmp (argv[n], "-sample") == 0) {
			sample = 1;
		} else if (strcmp (argv[n], "-cal") == 0) {
			export.cal_format = 1;
		} else if (strcmp (argv[n], "-hr") == 0) {
			export.use_hour = 1;
		} else if (strcmp (argv[n], "-date") == 0) {
			r = sscanf (argv[n+1], "%d:%d:%d", &syear, &smonth, &sday);
			if (r <= 0) {
				fprintf(stderr, "Can not get -date argument:\n");
				exit (1);
			}
			if (syear < 0 || smonth < 1 || smonth > 12 || sday < 1 || sday > 31) {
				fprintf (stderr, "Error in -date argument: %s\n", argv[n+1]);
				exit (1);
			}
			n++;
		} else if (strcmp (argv[n], "-user") == 0) {
			r = sscanf (argv[n+1], "%s", userfile);
			if (r <= 0) {
				fprintf (stderr, "Can't read user file name.\n");
				exit (1);
			}
			n++;
			export.user = 1;
		} else if (strcmp (argv[n], "-f") == 0) {
			r = sscanf (argv[n+1], "%s", outfile2);
			if (r <= 0) {
				fprintf (stderr, "Can't read -f output file name.\n");
				exit (1);
			}
			n++;
			export.outfile2 = 1;
		} else if (strcmp (argv[n], "-s") == 0) {
			r = sscanf (argv[n+1], "%s", outfile1);
			if (r <= 0) {
				fprintf (stderr, "Can't read -s output file name.\n");
				exit (1);
			}
			n++;
			export.outfile1 = 1;
		} else if (strcmp (argv[n], "-all") == 0) {
			all = 1;
		} else if (strcmp (argv[n], "-unique") == 0) {
			unique=1;
		} else if (strcmp (argv[n], "-orig") == 0) {
			export.orig = 1;
		} else if (strcmp (argv[n], "-func") == 0) {
			export.func = 1;
		} else if (strcmp (argv[n], "-poly") == 0) {
			export.poly = 1;
		} else if (strcmp (argv[n], "-smooth") == 0) {
			export.sm = 1;
		} else if (strcmp (argv[n], "-trend") == 0) {
			export.tr = 1;
		} else if (strcmp (argv[n], "-detrend") == 0) {
			export.detrend = 1;
		} else if (strcmp (argv[n], "-smcycle") == 0) {
			export.smcycle = 1;
		} else if (strcmp (argv[n], "-harm") == 0) {
			export.harm = 1;
		} else if (strcmp (argv[n], "-res") == 0) {
			export.res = 1;
		} else if (strcmp (argv[n], "-smres") == 0) {
			export.smres = 1;
		} else if (strcmp (argv[n], "-trres") == 0) {
			export.trres = 1;
		} else if (strcmp (argv[n], "-ressm") == 0) {
			export.ressm = 1;
		} else if (strcmp (argv[n], "-gr") == 0) {
			export.gr = 1;
		} else if (strcmp (argv[n], "-stats") == 0) {
			export.stats = 1;
		} else if (strcmp (argv[n], "-amp") == 0) {
			do_amp = 1;
		} else if (strcmp (argv[n], "-mm") == 0) {
			do_mm = 1;
		} else if (strcmp (argv[n], "-column") == 0) {
			do_column = 1;
		} else if (strcmp (argv[n], "-coef") == 0) {
			if (parsearg (argv[n+1], &begcoef, &endcoef)!=1) {
				fprintf(stderr, "Cannot get coefficient range:\n");
				exit (1);
			}
			do_coef = 1;
			n++;
		} else if (strcmp (argv[n], "-h") == 0) {
			usage();
			exit(0);
		} else {
			if (argv[n][0] == '-') {
				fprintf (stderr, "Unknown option %s.\n", argv[n]);
				usage();
				exit (1);
			} else 
				file = argv[n];
		}
		n++;
	}


	read_data(file, &filter);

#ifdef DEBUG
printf ("%d data points read.\n", np);
printf ("sample interval = %f\n", filter.sampleinterval);
#endif

	if (syear < 0 ) 
		GetCalendarDate (xp[0], &syear, &smonth, &sday, &hour, &minute, &second);

	export.startdate = DateValue (syear, smonth, sday, 12, 0, 0);

	compute_curves(&filter);

	export.format = " %.6e";



	if (export.user == 0 && equal==0) sample = 1;

	r = check_export(export, all, unique);
	if (r != 0) {
		if (export.outfile1) {
			fp = fopen (outfile1, "w");
			if (fp == NULL) {
				fprintf (stderr, "Can't open file %s for writing.\n", outfile1);
				exit (1);
			}
		} else {
			fp = stdout;
		}

		set_export (&export, all, unique, 0);
		if (sample) {
			export_sample_dates(fp, export, xp, yp, smooth2, trend2, deriv2, np, filter);
		}
		

		if (export.outfile2) {
			fp = fopen (outfile2, "w");
			if (fp == NULL) {
				fprintf (stderr, "Can't open file %s for writing.\n", outfile2);
				exit (1);
			}
		} else {
			fp = stdout;
		}

		set_export (&export, all, unique, 1);
		if (export.user) {
			export_user(userfile, fp, export, xfilt, smooth, trend, deriv, nfilt, filter);
		} else if (equal) {
			export_equally_spaced(fp, export, xfilt, smooth, trend, deriv, nfilt, filter);
		}
	}

	if (do_coef) 
		print_coef(begcoef, endcoef, filter);

	if (export.stats) 
		print_stats(file, filter);

	if (do_amp)
		print_ampstats(filter, do_column);

	if (do_mm)
		print_monthstats(filter, do_column);

	free_dmatrix(covar, 0, MAX_PARAM, 0, MAX_PARAM);
	exit(0);
}


/***************************************************************/
int read_data(file, filter)
char *file;
FilterData *filter;
{
	FILE *fp;
        char string[BUFSIZ];
        int k=0, n, sd, nh;
        double x,y, diff, sdiff, avginterval;


	if (file == NULL) 
		fp = stdin;
	else {
		fp = fopen (file, "r");
		if (fp == NULL) {
			fprintf (stderr, "Can not open file %s.\n", file);
				exit (1);
		}
	}

#ifdef DEBUG
if (file==NULL)
printf ("reading stdin\n");
else
printf ("reading %s\n", file);
#endif

        k=0; sd=0; sdiff=0.0;
        while ((fgets (string, BUFSIZ, fp)) != NULL) {
                n=sscanf (string, "%lf %lf", &x, &y);
                if (n!=2)  return (-1);
                xp[k] = x ; /*+ (1900*(a[0]<100)); */
                yp[k] = y;
                sig[k] = 1.0;
                if (k>0) {
                        diff = xp[k] - xp[k-1];
                        sdiff += diff;
                        sd++;
                }
                k++;
                if (k == MAX_OBS) {
                        fprintf (stderr, "Maximum number of data (%d) points reached.", MAX_OBS);
                        break;
                }
        }
        fclose (fp);

/* 
 * Calculate the average time interval between data points.
 * Set the sampleinterval variable if not set on the command line.
 */

	if (filter->sampleinterval == 0) {
		avginterval = sdiff/sd * 365.0;
		if (avginterval > 1.0)
			filter->sampleinterval = (float) NINT(avginterval);
		else
			filter->sampleinterval = avginterval;
	}

        np=k;

/* 
 * If the data is actually an average over a relatively large time period,
 * such as annual averages, change the number of harmonics to 
 * an appropriate value.
 */

	if (filter->timezero < 0) {
		filter->timezero = (int) xp[0]; /*-1900; */
	}
        nh = (int) 365.0/(filter->sampleinterval*2);
        if (nh < filter->numharmonics) {
                filter->numharmonics = nh;
        }


        return (0);
}

/************************************************************************/
void compute_curves (filter)
FilterData *filter;
{
	int i;
        double  t, u, p;
	float ave;


/*
	if ((xp[np-1]-xp[0]) < 2.0) filter->numpolyterms=2;
*/

/* compute filter
 * subtract date specified for time = 0 first, so that
 * coefficients of function fit are relative to this time.
 * Then add back this date so that arrays are relative to actual date */


	for (i=0; i<np; i++) {
		work[i] = xp[i] - filter->timezero;
/*		if (i < 10) printf ("%d %f\n",i, work[i]); */
	}
	filter_data (work, yp, sig, np, 
		filter->shortterm, filter->longterm, filter->sampleinterval, 
                filter->numpolyterms, filter->numharmonics, 
		filter->pm, covar, &filter->numpm,
                xfilt, smooth, trend, deriv, &nfilt, 
		&filter->varf1, &filter->varf2, &filter->chisq, &filter->rsd1);

	filter->funcvar = (float) varnce(covar, filter->numpm, xfilt[0]);
	for (i=0; i<nfilt; i++)
		xfilt[i] += filter->timezero;



/* compute smoothed results at every data point */
	if (smooth2) free(smooth2);
	smooth2 = calloc (np, sizeof(double));
	spline_f (xfilt, smooth, nfilt, xp, smooth2, work, np);

/* compute trend results at every data point */
	if (trend2) free(trend2);
	trend2 = calloc (np, sizeof(double));
	spline_f (xfilt, trend, nfilt, xp, trend2, work, np);

/* compute deriv results at every data point */
	if (deriv2) free(deriv2);
	deriv2 = calloc (np, sizeof(double));
	spline_f (xfilt, deriv, nfilt, xp, deriv2, work, np);

        for (i=0; i<np; i++) {
                t = xp[i] - filter->timezero;
                u = userfunc (filter->pm, t, filter->numharmonics, filter->numpolyterms); /* polynomial + harmonics */
                p = polyv(filter->pm, t, filter->numpolyterms); /* polynomial only */

                work[i] = yp[i] - (u + smooth2[i]);   /* residuals from smooth curve */
        }                                                                                                      

/* calculate the standard deviation of the residuals about the smooth curve */
        means (work, &ave, &filter->rsd2, np);


}

/**********************************************************************************/
void print_coef(beg, end, filter)
int beg, end;
FilterData filter;
{
	int i, e;

	e = (end > filter.numpm) ? filter.numpm : end; 
	for (i=beg; i<=e; i++) 
		printf (" %.6f", filter.pm[i]);

	printf ("\n");
}


/***********************************************************************/
int parsearg (arg, beg, end)
char *arg;
int *beg, *end;

/* get beginning and ending records from option argument.
   First argument must be a number, followed by a single character as a
   separator.  Non valid characters are 0-9, *, ;.  The last argument can
   be either a number or an *, in which case the end record is set to a large
   value, and the end of file is caught in the readdata routine   */
{
	char b2[4];
	int n, b1;

	n=sscanf (arg, "%d%*[^0123456789*]%[0123456789*]", &b1, b2);
	if (n == 1) {
		*beg = b1;
		*end = b1;
		return (1);
	}

	*beg = b1;
	if (strcmp(b2,"*") != 0) {
		n=sscanf (b2,"%d",&b1);
		if (n != 1) return (-1);
		*end = b1;
	}
	return (1);
}

/****************************************************************/
void print_stats(file, filter)
char *file;
FilterData filter;
{
#define rad *180.0/M_PI		/* radians to degrees conversion */
	int i, j, ix;
	double a, b, c, siga, sigtheta;
	FILE *fp;
	int year, month, day, hour, minute, second;

	if (np == 0) {
		printf ("No data points.  No Statistics available.");
		return;
	}

fp=stdout;

	GetCalendarDate(filter.timezero, &year, &month, &day, &hour, &minute, &second);

	fprintf (fp, "*****  Filter Statistics.  *****\n\n");

	if (file != NULL) 
	fprintf (fp, "File:   %s\n\n", file);

        fprintf (fp, "Beginning date:                 %.6f\n",xp[0]);
        fprintf (fp, "Ending date:                    %.6f\n",xp[np-1]);

	fprintf (fp, "Number of data points:          %d\n\n",np);
	fprintf (fp, "FUNCTION PARAMETERS\n");
	fprintf (fp, "Time = 0 on %d/%d/%d\n"        ,year, month, day);
	fprintf (fp, "Number of polynomial terms:     %d\n",filter.numpolyterms);
	fprintf (fp, "Number of harmonic terms:       %d\n",filter.numharmonics);
	fprintf (fp, "Total Number of parameters:     %d\n",filter.numpm);
	fprintf (fp, "------------------------------------------------------\n");
	fprintf (fp, "Parameter          Value          Standard Deviation\n");
	for (i=0; i<filter.numpm; i++) {
		fprintf (fp, "%5d %20.6f %20.6f\n",i, filter.pm[i], sqrt(covar[i][i]));
	}
	fprintf (fp, "------------------------------------------------------\n");
	fprintf (fp, "Harmonic   Amplitude  Std. Dev.    Phase (degrees)  Std. Dev.\n");
        for (i=1; i<=filter.numharmonics; i++) {
		ix = 2*(i-1)+filter.numpolyterms;
		a = filter.pm[ix];
		b = filter.pm[ix+1];
		c = a*a+b*b;
		siga = (a*a*covar[ix][ix]+b*b*covar[ix+1][ix+1])/c;
		sigtheta = (b*b*covar[ix][ix]+a*a*covar[ix+1][ix+1])/(c*c);
		fprintf (fp, "%5d %11.2f %10.2f %16.2f %12.2f\n", i, sqrt(c), sqrt(siga), 
							      atan2(b,a) rad, sqrt(sigtheta) rad);
	}
	fprintf (fp, "------------------------------------------------------\n");
	fprintf (fp, "Full covariance matrix:\n");
        for (i=0; i<filter.numpm; i++) {
		for (j=0; j<filter.numpm; j++)
			fprintf (fp, "%13.4e", covar[i][j]);
		fprintf (fp, "\n");
	}
	fprintf (fp, "------------------------------------------------------\n");
/*	fprintf (fp, "Chi squared value of function fit: %f\n",chisq); */
	fprintf (fp, "Residual standard deviation about function: %f\n",filter.rsd1);
	fprintf (fp, "------------------------------------------------------\n");
	fprintf (fp, "\n\nFILTER PARAMETERS\n");
	fprintf (fp, "Short term filter cutoff:       %d days\n",filter.shortterm);
	fprintf (fp, "Long term filter cutoff:        %d days\n",filter.longterm);
	fprintf (fp, "Sampling interval:              %g days\n\n",filter.sampleinterval);
	fprintf (fp, "Function Standard Deviation:          %8.4f\n",sqrt(filter.funcvar));
	fprintf (fp, "Short Term Filter Standard Deviation: %8.4f\n",sqrt(filter.varf1));
	fprintf (fp, "Long  Term Filter Standard Deviation: %8.4f\n",sqrt(filter.varf2));
	fprintf (fp, "Smoothed curve Standard Deviation:    %8.4f\n",sqrt(filter.varf1 + filter.funcvar));
	fprintf (fp, "Trend curve Standard Deviation:       %8.4f\n",sqrt(filter.varf2 + filter.funcvar));
	fprintf (fp, "Detrended Cycle Standard Deviation:   %8.4f\n",sqrt(filter.varf2 + filter.varf1 + 2*filter.funcvar));
	fprintf (fp, "Growth Rate Standard Deviation:       %8.4f\n\n",sqrt(2*(filter.varf2 + filter.funcvar)));
	fprintf (fp, "Residual standard deviation about smooth curve: %f\n",filter.rsd2);
}

/***************************************************************/
int check_export (export, all, unique)
ExportData export;
int all, unique;
{
	
	if (all || unique
		 || export.orig
		 || export.func
		 || export.poly
		 || export.sm 	
		 || export.tr 
		 || export.detrend
		 || export.smcycle
		 || export.harm
		 || export.res 
		 || export.smres
		 || export.trres
		 || export.ressm
		 || export.gr ) {
		return (1);
	} else {
		return (0);
	}
}

/***************************************************************/
void set_export (export, all, unique, which)
ExportData *export;
int all, unique, which;
{
	if (all == 0 && unique == 0) return;

	if (all) {
		export->orig = 1;
		export->func = 1;
		export->poly = 1;
		export->sm = 1;
		export->tr = 1;
		export->detrend = 1;
		export->smcycle = 1;
		export->harm = 1;
		export->res = 1;
		export->smres = 1;
		export->trres = 1;
		export->ressm = 1;
		export->gr = 1;
	} 
	if (unique) {
		if (which == 0) {  /* sample dates */
			export->orig = 1;
			export->func = 0;
			export->poly = 0;
			export->sm = 0;
			export->tr = 0;
			export->detrend = 1;
			export->smcycle = 0;
			export->harm = 0;
			export->res = 1;
			export->smres = 0;
			export->trres = 0;
			export->ressm = 1;
			export->gr = 0;
		}
		if (which == 1) {
			export->orig = 0;
			export->func = 1;
			export->poly = 1;
			export->sm = 1;
			export->tr = 1;
			export->detrend = 0;
			export->smcycle = 1;
			export->harm = 1;
			export->res = 0;
			export->smres = 1;
			export->trres = 1;
			export->ressm = 0;
			export->gr = 1;
		}
	} 
}



/****************************************************************/
void print_ampstats(filter, do_column)
FilterData filter;
int do_column;
{
	int 	i, na, yr, mm1, mm2, dd1, dd2, hr, nc1, nc2;
	int	m, s;
	float 	a;
	float	xwork[200], datemin[200], datemax[200], c1[200], c2[200];
	float	ampmin[200], ampmax[200];
	FILE 	*fp;
	char  	s1[10], s2[10];
	static char *mon[] = { "Jan", "Feb", "Mar", "Apr", "May", "Jun",
			       "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" };

	if (np == 0) {
		printf ("No data points.  No Statistics available.");
		return;
	}

	get_amplitudes (filter, xwork, ampmin, ampmax, datemin, datemax, &na, c1, c2, &nc1, &nc2);

fp=stdout;

	if (!do_column) {
		fprintf (fp, "*****  Season Cycle Statistics.  *****\n\n");
		fprintf (fp, " Year      Amplitude     Maximum   Date     Minimum   Date\n");
		fprintf (fp, "-----------------------------------------------------------\n");
		for (i=0; i<na; i++) {
			GetCalendarDate((double) datemin[i], &yr, &mm1, &dd1, &hr, &m, &s);
			GetCalendarDate((double) datemax[i], &yr, &mm2, &dd2, &hr, &m, &s);
			fprintf (fp, "%5.0f %12.2f %12.2f   %3s %2d %9.2f   %3s %2d\n",
					xwork[i], 
					ampmax[i]-ampmin[i], 
					ampmax[i], mon[mm2-1], dd2, 
					ampmin[i], mon[mm1-1], dd1);
		}
		fprintf (fp, "-----------------------------------------------------------\n");
		a = 2*(filter.varf2+filter.varf1+2*filter.funcvar);
		fprintf (fp, "The uncertainty of each amplitude (1 Standard Deviation) =  %8.4f\n\n",sqrt(a));
		fprintf (fp, "Trend Crossing Dates\n");
		fprintf (fp, "          - to +         + to -\n");
		fprintf (fp, "-----------------------------------------------------------\n");
		na = (nc1 > nc2) ? nc1 : nc2;
		for (i=0; i<na; i++) {
			if (i < nc1) {
				GetCalendarDate((double) c1[i], &yr, &mm1, &dd1, &hr, &m, &s);
				sprintf (s1, "%02d %3s %2d", dd1, mon[mm1-1], yr);
			} else {
				sprintf (s1, "         ");
			}
			if (i<nc2) {
				GetCalendarDate((double) c2[i], &yr, &mm2, &dd2, &hr, &m, &s);
				sprintf (s2, "%02d %3s %2d", dd2, mon[mm2-1], yr);
			} else {
				sprintf (s2, "         ");
			}
			fprintf (fp, "%5d %12s %12s\n", i+1, s1, s2);
		}
	} else {
		for (i=0; i<na; i++) {
/*
			GetCalendarDate((double) datemin[i], &yr, &mm1, &dd1, &hr, &m, &s);
			GetCalendarDate((double) datemax[i], &yr, &mm2, &dd2, &hr, &m, &s);
*/
			fprintf (fp, "%5.0f %7.2f %7.2f %7.2f\n",
					xwork[i], 
					ampmax[i]-ampmin[i], 
					ampmax[i], 
					ampmin[i]);
		}
	}
}

/***********************************************************************/
void get_amplitudes(filter, x, ampmin, ampmax, d_min, d_max, n, cdateup, cdatedown, nc1, nc2)
FilterData filter;
float x[], ampmin[], ampmax[];
float d_min[], d_max[];
float cdateup[], cdatedown[];
int *n, *nc1, *nc2;
{
/* get peak to peak amplitudes for each year */

	int 	i, na, nd1, nd2;
	float	year, tyear, a, amax=0.0, amin=1000.0;
	float	dmax = 0.0, dmin=0.0, ta = 0.0;
	double  t;

	na=0; nd1=0; nd2=0;
	*n = 0;
	tyear = (int) xfilt[0];
	for (i=0; i<nfilt; i++) {
		t = xfilt[i] - filter.timezero;
		year = (int) xfilt[i];
		a = userfunc (filter.pm, t, filter.numharmonics, filter.numpolyterms) - polyv (filter.pm, t, filter.numpolyterms) + smooth[i] - trend[i];
		if (i == 0) ta = a;
		if (year != tyear) {
			ampmin[na] = amin;
			ampmax[na] = amax;
			d_max[na] = dmax;
			d_min[na] = dmin;
			x[na] = tyear;
			na++;
			amax = 0.0;
			amin = 1000.0;
			tyear=year;
		}
		if (a > amax) {amax = a; dmax = xfilt[i];}
		if (a < amin) {amin = a; dmin = xfilt[i];}
		if (ta < 0.0 && a >= 0.0) { cdateup[nd1] = xfilt[i]; nd1++; }
		if (ta > 0.0 && a <= 0.0) { cdatedown[nd2] = xfilt[i]; nd2++; }
		ta = a;
			
	}
	ampmin[na] = amin;
	ampmax[na] = amax;
	d_max[na] = dmax;
	d_min[na] = dmin;
	x[na] = tyear;
	na++;
	*n = na;
	*nc1 = nd1;
	*nc2 = nd2;
}


/********************************************************************/
void print_monthstats(filter, do_column)
FilterData filter;
int do_column;
{
	int	year, hour, month, day, i, tyear, j;
	float	annsum, sum[12], a, t;
	int	count[12], anncount, min, sec;
	FILE 	*fp;
	double	u;

	for (i=0; i<12; i++) {
		sum[i] = 0.0;
		count[i] = 0;
	}
	annsum = 0.0;
	anncount = 0;

fp=stdout;

	if (!do_column) {
		fprintf (fp, "YEAR |    JAN     FEB     MAR     APR     MAY     JUN     JUL     AUG     SEP     OCT     NOV     DEC     ANNUAL\n");
		fprintf (fp,"================================================================================================================\n");
	}

/* -----------------------------------------------------------------
 * find appropriate month for each data point, 
 * and sum up.  After each year, print monthly averages
 * ----------------------------------------------------------------*/

	tyear = (int) xfilt[0];
	for (i=0; i<nfilt; i++) {
		GetCalendarDate (xfilt[i], &year, &month, &day, &hour, &min, &sec);
		if (year != tyear) {
			print_averages (fp, do_column, tyear, sum, count, annsum, anncount);

			for (j=0; j<12; j++) {
				sum[j] = 0.0;
				count[j] = 0;
			}
			annsum = 0.0;
			anncount = 0;
		}

		t = xfilt[i] - filter.timezero;
                u = userfunc (filter.pm, t, filter.numharmonics, filter.numpolyterms);
		a = u + smooth[i];
		sum[month-1] += a;
		count[month-1] += 1;

		annsum += a;
		anncount += 1;

		tyear=year;
	}

	print_averages (fp, do_column, tyear, sum, count, annsum, anncount);
	
}

/*********************************************************************/
static void print_averages (fp, do_column, tyear, sum, count, annsum, anncount)
FILE *fp;
int do_column;
float sum[], annsum;
int count[], tyear, anncount;
{
	int i;
	float mean, annual;

	if (do_column) {
		for (i=0; i<12; i++) {
			if (count[i] == 0) 
				mean = -9.99;
			else
				mean = sum[i]/count[i];
			fprintf (fp,"%4d %02d %7.2f\n",tyear, i+1, mean);
		}

	} else {
		fprintf (fp,"%4d |",tyear);
		for (i=0; i<12; i++) {
			if (count[i] == 0) 
				mean = -9.99;
			else
				mean = sum[i]/count[i];
			fprintf (fp," %7.2f",mean);
		}

		annual = annsum/anncount;
		fprintf (fp,"   %6.2f\n",annual);
	}
}