/*--------------------------------------------------------------------
 *	$Id: sample1d.c,v 1.3.4.2 2002/02/27 17:58:55 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * sample1d reads a 1-D dataset, and resamples the values on (1) user
 * supplied time values or (2) equidistant time-values based on <timestart> and
 * <dt>, both supplied at the command line. Choose among linear, cubic
 * spline, and Akima's spline.  sample1d will handle multiple column files,
 * user must choose which column contains the independent, monotonically
 * increasing variable.
 *
 * Author:	Paul Wessel
 * Date:	05-JUL--2000
 * Version:	3.4.1
 *
 */
 
#include "gmt.h"

main (int argc, char **argv)
{
	int i, j, k, t = 0, n, m, n_alloc, m_alloc, result, n_col, rows = 1, n_fields, n_expected_fields;
	int n_files = 0, fno, n_args, n_req, m_supplied;
	
	BOOLEAN error = FALSE, equidist = FALSE, t_supplied = FALSE, *nan_flag;
	BOOLEAN got_t0 = FALSE, got_dt = FALSE, nofile = TRUE, done = FALSE;
	
	double *t_supplied_out, *t_out, *ttime, *data, **col, **out;
	double dt = 0.0, t0 = 0.0, tt, low_t, high_t, *in, *dout;
	
	char line[BUFSIZ], type[3];
	
	FILE *fp = NULL, *fpt = NULL;
	
	type[0] = 'l';	type[1] = 'a';	type[2] = 'c';
	
	argc = GMT_begin (argc, argv);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'H':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;

				/* Supplemental parameters */
				
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'F':
					switch (argv[i][2]) {
						case 'l':
						case 'L':
							gmtdefs.interpolant = 0;
							break;
						case 'a':
						case 'A':
							gmtdefs.interpolant = 1;
							break;
						case 'c':
						case 'C':
							gmtdefs.interpolant = 2;
							break;
						default:	/* Use GMT defaults */
							break;
					}
					break;
				case 'I':
					dt = atof (&argv[i][2]);
					equidist = TRUE;
					got_dt = TRUE;
					break;
				case 'T':
					t = atoi (&argv[i][2]);
					break;
				case 'M':               /* Multiple line segments */
					GMT_multisegment (&argv[i][2]);
					break;
				case 'N':
					if ((fpt = fopen (&argv[i][2], "r")) == NULL) {
						fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, &argv[i][2]);
						exit (EXIT_FAILURE);
					}
					t_supplied = TRUE;
					break;
				case 'S':
					t0 = atof (&argv[i][2]);
					equidist = TRUE;
					got_t0 = TRUE;
					break;
					
				/* For backward compatibility for now */

				case 'L':
					gmtdefs.interpolant = 0;
					break;
				case 'A':
					gmtdefs.interpolant = 1;
					break;
				case 'C':
					gmtdefs.interpolant = 2;
					break;

				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
		
	}

	if (argc == 1 || GMT_quick) {	/* Display usage */
		fprintf (stderr,"sample1d %s - Resampling of 1-D data sets\n\n", GMT_VERSION);
		fprintf (stderr, "usage: sample1d <infile(s)> [-Fl|a|c] [-H[<nrec>]] [-I<t_inc>] [-M[<flag>]] [-N<knotfile>]\n");
		fprintf (stderr, "\t[-S<xstart>] [-T<time_col>] [-V] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<infile> is one or more multicolumn ASCII (or binary, see -b) tables. [Default is standard input]\n");
		fprintf (stderr, "\tThe independent variable (see -T) must be monotonically in/de-creasing\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-F sets the interpolation mode.  Choose from:\n");
		fprintf (stderr, "\t   l Linear interpolation\n");
		fprintf (stderr, "\t   a Akima spline interpolation\n");
		fprintf (stderr, "\t   c Cubic spline interpolation\n");
		fprintf (stderr, "\t   [Default is -F%c]\n", type[gmtdefs.interpolant]);
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I <x_inc> sets equidistant grid interval\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N <knotfile> is an ASCII table with the desired time positions in column 0\n");
		fprintf (stderr, "\t   Overrides the -I and -S settings.  If none of -I, -S, and -N is set\n");
		fprintf (stderr, "\t   then <tstart> = first input point, <t_inc> = (t[1] - t[0])\n");
		fprintf (stderr, "\t-S <xstart> sets the first output point\n");
		fprintf (stderr, "\t-T gives column number of the independent variable (time) [Default is 0 (first)]\n");
		GMT_explain_option ('V');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		GMT_explain_option ('o');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (t < 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -T option: Column number cannot be negative\n", GMT_program);
		error++;
	}
	if (t_supplied && equidist) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR: Specify only one of -N and -S\n", GMT_program);
		error++;
	}
	if (!t_supplied && !equidist) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify one of -N and -S\n", GMT_program);
		error++;
	}
	if (got_dt && dt <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option: Must specify positive increment\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[0] && gmtdefs.io_header) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
        if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 2;
	n_req = (t >= 2) ? t + 1 : 2;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < n_req) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least %d columns\n", GMT_program, n_req);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */
	
	if (GMT_io.binary[0] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	t_out = (double *)NULL;

	if (t_supplied) {	/* read file with abscissae */
		m_alloc = GMT_CHUNK;
		t_supplied_out = (double *) GMT_memory (VNULL, (size_t)m_alloc, sizeof (double), GMT_program);
		m = 0;
		if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, BUFSIZ, fpt);
		while (fgets (line, BUFSIZ, fpt)) {
			if (line[0] == '#' || line[0] == GMT_io.EOF_flag) continue;
			sscanf (line, "%lf", &t_supplied_out[m]);
			m++;
			if (m == m_alloc) {	/* Get more memory */
				m_alloc += GMT_CHUNK;
				t_supplied_out = (double *) GMT_memory ((void *)t_supplied_out, (size_t)m_alloc, sizeof (double), GMT_program);
			}
		}
		fclose (fpt);
		m_supplied = m;
		t_supplied_out = (double *) GMT_memory ((void *)t_supplied_out, (size_t)m_supplied, sizeof (double), GMT_program);
		t_out = (double *) GMT_memory (VNULL, (size_t)m_supplied, sizeof (double), GMT_program);
                if (gmtdefs.verbose) fprintf (stderr, "%s: Read %d knots from file\n", GMT_program, m_supplied);
	}

	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;

	n_args = (argc > 1) ? argc : 2;

	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over input files, if any */
		if (!nofile && argv[fno][0] == '-') continue;
		
		if (nofile) {	/* Just read standard input */
			fp = GMT_stdin;
			done = TRUE;
#ifdef SET_IO_MODE
			GMT_setmode (0);
#endif
		}
		else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) {
			fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]);
			continue;
		}

		if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]);

		if (gmtdefs.io_header) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (line, BUFSIZ, fp);
				fprintf (GMT_stdout, "%s", line);
			}
		}
	
		n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : BUFSIZ;
		i = n_col = 0;
		n_fields = GMT_input (fp, &n_expected_fields, &in);
	
		while (! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

			while (GMT_io.status & GMT_IO_SEGMENT_HEADER) {
				GMT_write_segmentheader (GMT_stdout, n_expected_fields);
				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
			if (GMT_io.status & GMT_IO_EOF) continue;	/* At EOF */

			i++;
			n = 0;
	
			while (! (GMT_io.status & (GMT_IO_SEGMENT_HEADER | GMT_IO_EOF))) {	/* Keep going until FALSE or = 2 segment header */
				if (GMT_io.status & GMT_IO_MISMATCH) {
					fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, i);
					exit (EXIT_FAILURE);
				}

				if (n_col == 0) {	/* Allocate memory first time around */
					n_col = n_expected_fields;
					if (t >= n_col) {
						fprintf (stderr, "%s: time_col = %d exceeds range of columns (%d)!\n", GMT_program, t, n_col);
						exit (EXIT_FAILURE);
					}
					dout = (double *) GMT_memory (VNULL, (size_t)n_col, sizeof (double), GMT_program);
					col = (double **) GMT_memory (VNULL, (size_t)n_col, sizeof (double *), GMT_program);
					out = (double **) GMT_memory (VNULL, (size_t)n_col, sizeof (double *), GMT_program);
					nan_flag = (BOOLEAN *) GMT_memory (VNULL, (size_t)n_col, sizeof (BOOLEAN), GMT_program);
					n_alloc = GMT_CHUNK;
					for (j = 0; j < n_col; j++) col[j] = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
				}
		
				for (j = 0; j < n_col; j++) {
					if (GMT_is_dnan (in[j])) {
						nan_flag[j] = TRUE;
						col[j][n] = GMT_d_NaN;
					}
					else
						col[j][n] = in[j];
				}
				n++;
		
				if (n == n_alloc) {	/* Get more memory */
					n_alloc += GMT_CHUNK;
					for (j = 0; j < n_col; j++) col[j] = (double *) GMT_memory ((void *)col[j], (size_t)n_alloc, sizeof (double), GMT_program);
				}

				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
			
			/* Here we have one segment to work on */

			/* If we didn't get input abscissa, now's the time to generate them */

			if (t_supplied) {	/* Get relevant t_out segment */
				low_t  = MIN (col[t][0], col[t][n-1]);
				high_t = MAX (col[t][0], col[t][n-1]);
				for (i = m = 0; i < m_supplied; i++) {
					if (t_supplied_out[i] < low_t || t_supplied_out[i] > high_t) continue;
					t_out[m++] = t_supplied_out[i];
				}
				if (m == 0) fprintf (stderr, "%s: Warning: No output points for range %lg to %lg\n", GMT_program, col[t][0], col[t][n-1]);
			}
			else {	/* Generate evenly spaced grid */
				if (!got_dt) dt = col[t][1] - col[t][0];
				if (got_dt && (col[t][1] - col[t][0]) < 0.0 && dt > 0.0) dt = -dt;	/* For monotonically decreasing data */
				if (!got_t0) {
					if (dt > 0.0) {
						t0 = floor (col[t][0] / dt) * dt;
						if (t0 < col[t][0]) t0 += dt;
					}
					else {
						t0 = ceil (col[t][0] / (-dt)) * (-dt);
						if (t0 > col[t][0]) t0 += dt;
					}
				}
				m = m_alloc = irint (fabs((col[t][n-1] - t0) / dt)) + 1;
				t_out = (double *) GMT_memory ((void *)t_out, (size_t)m_alloc, sizeof (double), GMT_program);
				t_out[0] = t0;
				i = 1;
				if (dt > 0.0) {
					while (i < m && (tt = t0 + i * dt) <= col[t][n-1]) {
						t_out[i] = tt;
						i++;
					}
				}
				else {
					while (i < m && (tt = t0 + i * dt) >= col[t][n-1]) {
						t_out[i] = tt;
						i++;
					}
				}
				m = i;
				if (fabs (t_out[m-1]-col[t][n-1]) < SMALL) {	/* Fix roundoff */
					t_out[m-1] = col[t][n-1];
				}
			}

			if (nan_flag[t]) {
				fprintf (stderr, "%s: Independent column has NaN's!\n", GMT_program);
				exit (EXIT_FAILURE);
			}

			for (j = 0; m && j < n_col; j++) {

				if (j == t) continue;	/* Skip the time column */

				out[j] = (double *) GMT_memory (VNULL, (size_t)m, sizeof (double), GMT_program);
		
				if (nan_flag[j]) {	/* NaN's present, need "clean" time and data columns */

					ttime = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
					data = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
					for (i = k = 0; i < n; i++) {
						if ( GMT_is_dnan (col[j][i]) ) continue;
						ttime[k] = col[t][i];
						data[k] = col[j][i];
						k++;
					}
					result = GMT_intpol (ttime, data, k, m, t_out, out[j], gmtdefs.interpolant);
					GMT_free ((void *)ttime);
					GMT_free ((void *)data);
				}
				else
					result = GMT_intpol (col[t], col[j], n, m, t_out, out[j], gmtdefs.interpolant);
				
				if (result != 0) {
					fprintf (stderr, "%s: Error from GMT_intpol near row %d!\n", GMT_program, rows+result+1);
					exit (EXIT_FAILURE);
				}
			}
	
			out[t] = t_out;

			for (i = 0; i < m; i++) {
				for (j = 0; j < n_col; j++) dout[j] = out[j][i];
				GMT_output (GMT_stdout, n_col, dout);
			}
			for (j = 0; m && j < n_col; j++) if (j != t) GMT_free ((void *)out[j]);
		
			rows += n + 1;
		}
	
		if (fp != GMT_stdin) GMT_fclose(fp);
	}

	for (j = 0; j < n_col; j++) GMT_free ((void *)col[j]);

	GMT_free ((void *)t_out);
	GMT_free ((void *)out);
	GMT_free ((void *)col);
	
	GMT_end (argc, argv);
}