/*-----------------------------------------------------------------
 *	$Id: x2sys_binlist.c,v 1.9 2004/08/27 21:30:08 pwessel Exp $
 *
 *      Copyright (c) 1999-2004 by P. Wessel
 *      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: www.soest.hawaii.edu/pwessel
 *--------------------------------------------------------------------*/
/* x2sys_binlist will read one or several data files and dump their
 * contents to stdout in ascii or binary (double precision) mode.
 * Input data file formats are determined by the definition file
 * given by the -D option.
 *
 * Author:	Paul Wessel
 * Date:	15-JUN-2004
 * Version:	1.1, based on the spirit of the old xsystem code
 *
 */

#include "x2sys.h"

int outside (double x, double y, struct X2SYS_BIX *B, BOOLEAN geo);
unsigned int get_data_flag (double *data[], int j, struct X2SYS_INFO *s);

main (int argc, char **argv)
{
	char *TAG = CNULL;

	int i, j, this_j, last_j, ij, bi, bj, ii, jj, start_i, end_i;
	int last_bin_i;	/* Previous i node for bin */
	int last_bin_j;	/* Previous j node for bin */
	int last_bin_ij;/* Previous bin */
	int this_bin_i;	/* This i node for bin */
	int this_bin_j;	/* This j node for bin */
	int this_bin_ij;/* This bin */
	int nx;
	unsigned int nav_flag;

	BOOLEAN error = FALSE;

	double **data, time_gap = DBL_MAX, x, y, dx, dx2, xc[2], yc[2];

	struct X2SYS_INFO *s;
	struct X2SYS_FILE_INFO p;		/* File information */
	struct X2SYS_BIX B;
	
	GMT_begin (argc, argv);
	for (i = strlen(argv[0]); i >= 0 && argv[0][i] != '/'; i--);
	X2SYS_program = &argv[0][i+1];	/* Name without full path */

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* Common parameters */
			
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], NULL, NULL, NULL, NULL);
					break;
				
				/* Supplemental parameters */
				
				case 'T':
					TAG = &argv[i][2];
					break;
				default:
					error = TRUE;
					break;
			}
		}
	}

	if (argc == 1 || error || GMT_quick) {
		fprintf (stderr, "x2sys_binlist %s - Create bin index listing from data files\n\n", X2SYS_VERSION);
		fprintf (stderr, "usage: x2sys_binlist <files> -T<TAG> [-V]\n\n");

		if (GMT_quick) exit (EXIT_FAILURE);

		fprintf (stderr, "	<files> is one or more datafiles\n");
		fprintf (stderr, "	-T <TAG> is the system tag for this compilation\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}
	
	x2sys_set_system (TAG, &s, &B, &GMT_io);
		
	if (s->geographic) {
		GMT_io.out_col_type[0] = GMT_IS_LON;
		GMT_io.out_col_type[1] = GMT_IS_LAT;
	}
	else
		GMT_io.out_col_type[0] = GMT_io.out_col_type[1] = GMT_IS_FLOAT;

	x2sys_bix_init (&B, TRUE);
	nav_flag = (1 << s->x_col) + (1 << s->y_col);	/* For bins just cut by track but no points inside the bin */

	fprintf (GMT_stdout, "# %s\n", TAG);
	
	for (i = 1; i < argc; i++) {

		if (argv[i][0] == '-') continue;

		if (gmtdefs.verbose) fprintf (stderr, "x2sys_binlist: Reading file %s\n", argv[i]);

		if (((s->read_file) (argv[i], &data, s, &p, &GMT_io)) < 0) {
			fprintf (stderr, "x2sys_binlist: Error reading file %s\n", argv[i]);
			exit (EXIT_FAILURE);
		}

		/* Reset bin flags */
		
		memset ((void *)B.binflag, 0, (size_t)(B.nm_bin * sizeof (unsigned int)));
		
		j = 0;
		while (j < p.n_rows && outside (data[s->x_col][j], data[s->y_col][j], &B, s->geographic)) j++;	/* Goto first point inside region */
		
		this_bin_ij = x2sys_bix_get_ij (data[s->x_col][j], data[s->y_col][j], &this_bin_i, &this_bin_j, &B);
		B.binflag[this_bin_ij] |= get_data_flag (data, j, s);
		this_j = j;

		for (j = this_j + 1; j < p.n_rows; j++) {
			if (outside (data[s->x_col][j], data[s->y_col][j], &B, s->geographic)) continue;
			
			last_bin_i = this_bin_i;
			last_bin_j = this_bin_j;
			last_bin_ij = this_bin_ij;
			last_j = this_j;
			this_j = j;
			this_bin_ij = x2sys_bix_get_ij (data[s->x_col][j], data[s->y_col][j], &this_bin_i, &this_bin_j, &B);
			
			/* While this may be the same bin as the last bin, the data available may have changed so we keep
			 * turning the data flags on again and again. */
			 
			B.binflag[this_bin_ij] |= get_data_flag (data, j, s);
			
			/* If there is time we will check for time gaps */
			if (s->t_col >= 0 && (data[s->t_col][j] - data[s->t_col][j-1]) > time_gap) continue;

			/* If any of the code inside these 2 for-loops are executed it is because the track between
			 * this and previous point cuts a neighboring bin so we include those bins as well.  However
			 * There are no data in those bins so we only assign the navigation flag */
			
			if (s->geographic && (this_bin_i - last_bin_i) > 180) {		/* Jumped from east to west across Greenwich */
				start_i = this_bin_i + 1;
				end_i = last_bin_i + 360;
				dx = (data[s->x_col][this_j] - 360.0) - data[s->x_col][last_j];
			}
			else if (s->geographic && (this_bin_i - last_bin_i) < -180) {	/* Jumped from west to east across Greenwich */
				start_i = last_bin_i + 1;
				end_i = this_bin_i + 360;
				dx = data[s->x_col][this_j] - (data[s->x_col][last_j] - 360.0);
			}
			else {								/* Did no such thing */
				start_i = MIN (last_bin_i, this_bin_i) + 1;
				end_i = MAX (last_bin_i, this_bin_i);
				dx = data[s->x_col][this_j] - data[s->x_col][last_j];
			}
			nx = 0;
			for (bj = MIN (last_bin_j, this_bin_j) + 1; bj <= MAX (last_bin_j, this_bin_j); bj++) {	/* If we go in here we know dy is non-zero */
				y = B.y_min + bj * B.bin_y;
				x = data[s->x_col][last_j] + (y - data[s->y_col][last_j]) * dx / (data[s->y_col][this_j] - data[s->y_col][last_j]);
				ij = x2sys_bix_get_ij (x, y, &ii, &jj, &B);
				B.binflag[ij] |= nav_flag;		/* Only update nav flags we havent been here already */
				if (nx < 2) {
					xc[nx] = x;	yc[nx] = y;
					nx++;
				}
			}
			for (bi = start_i; bi <= end_i; bi++) {	/* If we go in here we know dx is non-zero */
				x = B.x_min + bi * B.bin_x;
				if (s->geographic && x >= 360.0) x -= 360.0;
				dx2 = x - data[s->x_col][last_j];
				if (fabs (dx2) > 180.0) dx2 = copysign (360.0 - fabs (dx2), -dx2);
				y = data[s->y_col][last_j] + dx2 * (data[s->y_col][this_j] - data[s->y_col][last_j]) / dx;
				ij = x2sys_bix_get_ij (x, y, &ii, &jj, &B);
				B.binflag[ij] |= nav_flag;		/* Only update nav flags we havent been here already */
				if (nx < 2) {
					xc[nx] = x;	yc[nx] = y;
					nx++;
				}
			}
			if (nx == 2) {	/* Cut a corner, try half-point */
				dx = xc[0] - xc[1];
				if (s->geographic && dx < -180.0)
					x = 0.5 * (xc[0] + (xc[1] - 360.0));
				else if (s->geographic && dx > +180.0)
					x = 0.5 * (xc[0] - 360.0 + xc[1]);
				else
					x = 0.5 * (xc[0] + xc[1]);
				y = 0.5 * (yc[0] + yc[1]);
				ij = x2sys_bix_get_ij (x, y, &ii, &jj, &B);
				B.binflag[ij] |= nav_flag;		/* Only update nav flags we havent been here already */
			}			
		}

		x2sys_free_data (data, s->n_fields);
	
		/* Time for bin index output */
	
		fprintf (GMT_stdout, "> %s\n", argv[i]);
		for (ij = 0; ij < B.nm_bin; ij++) {
			if (B.binflag[ij] == 0) continue;
			x = B.x_min + ((ij % B.nx_bin) + 0.5) * B.bin_x;
			y = B.y_min + ((ij / B.nx_bin) + 0.5) * B.bin_y;
			GMT_ascii_output_one (GMT_stdout, x, 0);
			fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
			GMT_ascii_output_one (GMT_stdout, y, 1);
			fprintf (GMT_stdout, "%s%d%s%u\n", gmtdefs.field_delimiter, ij, gmtdefs.field_delimiter, B.binflag[ij]);
		}
	}

	x2sys_free_info (s);
	GMT_free ((void *)B.binflag);

	GMT_end (argc, argv);
}

int outside (double x, double y, struct X2SYS_BIX *B, BOOLEAN geo)
{
	if (y < B->y_min || y > B->y_max) return (1);
	if (geo) {	/* Geographic data with periodic longitudes */
		while (x < B->x_min) x += 360.0;
		while (x > B->x_max) x -= 360.0;
		if (x < B->x_min) return (1);
	}
	else {	/* Plain Cartesian test */
		if (x < B->x_min || x > B->x_max) return (1);
	}
	return (0);	/* Inside */
}

unsigned int get_data_flag (double *data[], int j, struct X2SYS_INFO *s)
{
	int i;
	unsigned int bit, flag;
	for (i = flag = 0, bit = 1; i < s->n_fields; i++, bit <<= 1) {
		if (GMT_is_dnan (data[i][j])) continue;	/* NaN, so no data here */
		flag |= bit;
	}
	return (flag);
}

