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

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

/*******************************************************************************/
void export_sample_dates(
FILE *fp,
ExportData export,
double xp[],
double yp[],
double smooth2[],
double trend2[],
double deriv2[],
int np,
FilterData filter
)

{

/*-------------------------------------
 * export data at same dates as 
 * original samples
 *-------------------------------------*/
	static char *format;
	int 	i, year, month, day, hour, minute, second;
	double  p, f, h, t;

	format = export.format;
	for (i=0; i<np; i++) {
		if (export.cal_format) {
			GetCalendarDate(xp[i], &year, &month, &day, &hour, &minute, &second);
			if (export.use_hour) {
				fprintf (fp, "%4d %2d %2d %2d",year, month, day, hour);
			} else {
				fprintf (fp, "%4d %2d %2d",year, month, day);
			}
		} else {
			fprintf (fp, "%13.8lf",xp[i]);
		}

		t = xp[i] - filter.timezero;
		p = userfunc (filter.pm, t, 0, filter.numpolyterms);
		f = userfunc (filter.pm, t, filter.numharmonics, filter.numpolyterms);
		h = f-p;
		if (export.orig)  	fprintf (fp, format, yp[i]);			/* data */
		if (export.func)  	fprintf (fp, format, f);			/* function */
		if (export.poly)  	fprintf (fp, format, p);			/* polynomial */
		if (export.sm) 		fprintf (fp, format, smooth2[i]+f);		/* smoothed curve */
		if (export.tr) 		fprintf (fp, format, trend2[i]+p);		/* trend */
		if (export.detrend) 	fprintf (fp, format, yp[i]-trend2[i]-p);	/* detrended data */
		if (export.smcycle) 	fprintf (fp, format, smooth2[i]+h-trend2[i]);/* detrended cycle */
		if (export.harm) 	fprintf (fp, format, h);			/* harmonics */
		if (export.res) 	fprintf (fp, format, yp[i] - f);		/* residuals */
		if (export.smres) 	fprintf (fp, format, smooth2[i]);		/* smoothed residuals */
		if (export.trres) 	fprintf (fp, format, trend2[i]);		/* trend of residuals */
		if (export.ressm) 	fprintf (fp, format, yp[i]-(f+smooth2[i]));	/* residuals from smoothed curve */
		if (export.gr) 		fprintf (fp, format, deriv2[i]);		/* growth rate */
		fprintf (fp, "\n");
	}
}

/**********************************************************************************/
/*----------------------------------------------
 * export data at equally spaced intervals
 *----------------------------------------------*/
/**********************************************************************************/
void export_equally_spaced(
FILE *fp,
ExportData export, 
double xfilt[], 
double smooth[], 
double trend[], 
double deriv[], 
int nfilt, 
FilterData filter
)
{
	char	*format;
	int 	i, year, month, day, hour, minute, second, j;
	double 	rinterval, dp;
	double  p, f, h, t;
	double	t0, t1, s1, d1, *s2, *t2, *d2;


	format = export.format;
	rinterval = filter.sampleinterval/365.0;
	i=0;

	s2 = (double *) calloc (nfilt, sizeof(double));
	if (s2 == NULL) {
		fprintf (stderr, "Can't allocate memory to interpolate values.");
		return;
	}
	t2 = (double *) calloc (nfilt, sizeof(double));
	if (t2 == NULL) {
		fprintf (stderr,"Can't allocate memory to interpolate values.");
		return;
	}
	d2 = (double *) calloc (nfilt, sizeof(double));
	if (d2 == NULL) {
		fprintf (stderr, "Can't allocate memory to interpolate values.");
		return;
	}

	spline (xfilt, smooth, nfilt, s2);
	spline (xfilt, trend, nfilt, t2);
	spline (xfilt, deriv, nfilt, d2);

	while ((t=export.startdate+i*rinterval) <= xfilt[nfilt-1]) {
/*printf ("t = %f, rinterval = %f, i = %d\n", t, rinterval, i); */
		if (export.cal_format) {
			GetCalendarDate (t, &year, &month, &day, &hour, &minute, &second);
			if (export.use_hour) {
				fprintf (fp, "%4d %2d %2d %2d",year, month, day, hour);
			} else {
				fprintf (fp, "%4d %2d %2d",year, month, day);
			}
		} else {
			fprintf (fp, "%13.8f",t);
		}

		t0 = t - filter.timezero;
		p = userfunc (filter.pm, t0, 0, filter.numpolyterms);
		f = userfunc (filter.pm, t0, filter.numharmonics, filter.numpolyterms);
		h = f-p;

		s1 = interpolate(xfilt, smooth, nfilt, s2, t);
		t1 = interpolate(xfilt, trend, nfilt, t2, t);

		if (export.func)    fprintf (fp, format, f);
		if (export.poly)    fprintf (fp, format, p);
		if (export.sm)      fprintf (fp, format, f + s1);
		if (export.tr)      fprintf (fp, format, p+t1);
		if (export.smcycle) fprintf (fp, format, h+s1-t1);
		if (export.harm)    fprintf (fp, format, h);
		if (export.smres)   fprintf (fp, format, s1);
		if (export.trres)   fprintf (fp, format, t1);
		if (export.gr) {
			if (t < xfilt[0] || t > xfilt[nfilt-1]) {
			/* calculate derivative of polynomial */
				if (filter.numpolyterms > 0 ) {
					p=filter.pm[filter.numpolyterms-1];
					dp=0.0;
					for (j=filter.numpolyterms-2;j>=0; j--) {
						dp = dp*t0+p; 
						p = p*t0+filter.pm[j];
					}
				} else {
					dp=0;
				}
				fprintf (fp, format, dp);
			} else {
				d1 = interpolate(xfilt, deriv, nfilt, d2, t);
				fprintf (fp, format, d1);
			}
		}

		fprintf (fp, "\n");
		i++;
	}
	free (s2);
	free (t2);
	free (d2);
}
/**********************************************************************************/
void export_user(
char *datefile,
FILE *fp,
ExportData export, 
double xfilt[], 
double smooth[], 
double trend[], 
double deriv[], 
int nfilt, 
FilterData filter
)
{
	char	*format, s[256];
	int 	i, year, month, day, hour, minute, second, j;
	double 	rinterval, dp;
	double  p, f, h, t;
	double	t0, t1, s1, d1, *s2, *t2, *d2;
	FILE 	*infile;


	infile = fopen (datefile, "r");
	if ( infile == NULL) {
		fprintf (stderr, "Can't open %s.\n", datefile);
		return;
	}

	format = export.format;
	rinterval = filter.sampleinterval/365.0;
	i=0;

	s2 = (double *) calloc (nfilt, sizeof(double));
	if (s2 == NULL) {
		fprintf (stderr, "Can't allocate memory to interpolate values.");
		return;
	}
	t2 = (double *) calloc (nfilt, sizeof(double));
	if (t2 == NULL) {
		fprintf (stderr,"Can't allocate memory to interpolate values.");
		return;
	}
	d2 = (double *) calloc (nfilt, sizeof(double));
	if (d2 == NULL) {
		fprintf (stderr, "Can't allocate memory to interpolate values.");
		return;
	}

	spline (xfilt, smooth, nfilt, s2);
	spline (xfilt, trend, nfilt, t2);
	spline (xfilt, deriv, nfilt, d2);


	while (fgets (s, 256, infile) != NULL) {
		sscanf (s, "%lf", &t);
/*printf ("t = %f, rinterval = %f, i = %d\n", t, rinterval, i); */
		if (export.cal_format) {
			GetCalendarDate (t, &year, &month, &day, &hour, &minute, &second);
			if (export.use_hour) {
				fprintf (fp, "%4d %2d %2d %2d",year, month, day, hour);
			} else {
				fprintf (fp, "%4d %2d %2d",year, month, day);
			}
		} else {
			fprintf (fp, "%13.8f",t);
		}

		t0 = t - filter.timezero;
		p = userfunc (filter.pm, t0, 0, filter.numpolyterms);
		f = userfunc (filter.pm, t0, filter.numharmonics, filter.numpolyterms);
		h = f-p;

		s1 = interpolate(xfilt, smooth, nfilt, s2, t);
		t1 = interpolate(xfilt, trend, nfilt, t2, t);

		if (export.func)    fprintf (fp, format, f);
		if (export.poly)    fprintf (fp, format, p);
		if (export.sm)      fprintf (fp, format, f + s1);
		if (export.tr)      fprintf (fp, format, p+t1);
		if (export.smcycle) fprintf (fp, format, h+s1-t1);
		if (export.harm)    fprintf (fp, format, h);
		if (export.smres)   fprintf (fp, format, s1);
		if (export.trres)   fprintf (fp, format, t1);
		if (export.gr) {
			if (t < xfilt[0] || t > xfilt[nfilt-1]) {
			/* calculate derivative of polynomial */
				if (filter.numpolyterms > 0 ) {
					p=filter.pm[filter.numpolyterms-1];
					dp=0.0;
					for (j=filter.numpolyterms-2;j>=0; j--) {
						dp = dp*t0+p; 
						p = p*t0+filter.pm[j];
					}
				} else {
					dp=0;
				}
				fprintf (fp, format, dp);
			} else {
				d1 = interpolate(xfilt, deriv, nfilt, d2, t);
				fprintf (fp, format, d1);
			}
		}

		fprintf (fp, "\n");
		i++;
	}
	free (s2);
	free (t2);
	free (d2);
}