/*-----------------------------------------------------------------
 *	$Id: x2sys_cross.c,v 1.14 2004/08/19 18:56:00 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_cross will calculate crossovers generated by the
 * intersections of two tracks.  Optionally, it will also evaluate
 * the interpolated datafields at the crossover locations.
 *
 * Author:	Paul Wessel
 * Date:	15-JUN-2004
 * Version:	1.0, based on the spirit of the old xsystem code,
 *		but with a smarter algorithm based on the book
 *		"Algorithms in C" by R. Sedgewick.
 *
 */

#include "x2sys.h"

struct PAIR {				/* Used with -Kkombinations.lis option */
	char *id1, *id2;
} *pair;

BOOLEAN combo_ok (char *name_1, char *name_2, struct PAIR *pair, int n_pairs);

main (int argc, char **argv)
{
	char *sfile, *def = "x2sys";		/* Name (and default value) of the definition file */
	char **file;				/* Filename for leg */
	char line[BUFSIZ];			/* Input buffer */
	char *header_format;			/* Pointer to header format */
	char t_or_i;				/* t = time, i = dummy node time */
	char *fp_list = CNULL;			/* File pointer to data list */
	char *fp_combo = CNULL;			/* File pointer to combinations list */
	char *fflags = CNULL;			/* Pointer to list of column names (in -F) */
	char name1[80], name2[80];		/* Name of two files to be examined */
	char *TAG = CNULL;			/* System TAG name for these kinds of data */

	int n_rec[2];				/* Number of data records for both files */
	int half_window_width = 3;		/* Number of points on either side in the interpolation */
	int window_width;			/* Max number of points to use in the interpolation */
	int n_legs = 0;				/* Total number of data sets to compare */
	int nx;					/* Number of crossovers found for this pair */
	int *col_number;			/* Array with the column numbers of the data fields */
	int n_output;				/* Number of columns on output */
	int dist_flag = 0;			/* 0 = Cartesian, 1 = Flat Earth, 2 = Spherical */
	int n_pairs = 0;			/* Number of acceptable combinations */
	int geodetic = 2;			/* 0 for [0,360>, 1 for <-360,0], else [-180,+180]  */
	int n_bin_header;
	int A, B, i, j, col, k, start;		/* Misc. counters and local variables */
	int end, first, n_ok, *ok, n_alloc = 1;
	int n_data_col, left[2], t_left;
	int n_left, right[2], t_right, n_right;

	size_t col_alloc;

	BOOLEAN xover_locations_only = FALSE;	/* TRUE if only x,y (and possible indices) to be output */
	BOOLEAN internal = TRUE;		/* FALSE if only external xovers are needed */
	BOOLEAN external = TRUE;		/* FALSE if only internal xovers are needed */
	BOOLEAN error = FALSE;			/* TRUE for invalid arguments */
	BOOLEAN geographic = FALSE;		/* Data coordinates are geographic locations */
	BOOLEAN do_project = FALSE;		/* TRUE if we must mapproject first */
	BOOLEAN got_time = FALSE;		/* TRUE if there is a time column */
	BOOLEAN xover_mode = FALSE;		/* Backwards compatibility with xover */
	BOOLEAN speed_check = FALSE;		/* TRUE if speed must be checked */
	BOOLEAN first_header = TRUE;		/* TRUE for very first crossover */
	BOOLEAN first_crossover;		/* TRUE for first crossover between two data sets */
	BOOLEAN combinations = FALSE;		/* TRUE if a combination list has been given with -K */
	BOOLEAN two_values = FALSE;		/* TRUE if we want the two track values rather than crossover and mean returned */
	BOOLEAN t_given = FALSE;		/* TRUE when -T aoption is set */
	BOOLEAN d_given = FALSE;		/* TRUE when -D aoption is set */

	double max_time_separation = DBL_MAX;	/* Largest time separation between two points in timeunits */
	double dt;				/* Time between crossover and previous node */
	double dist_x[2];			/* Distance(s) along track at the crossover point */
	double time_x[2];			/* Time(s) along track at the crossover point */
	double delx, dely, deld, delt;		/* Differences in x, y, dist, and time across xpoint */
	double speed;				/* speed across the xpoint ( = deld/delt) */
	double **data[2];			/* Data matrices for the two data sets to be checked */
	double *xdata[2];			/* Data vectors with estimated values at crossover points */
	double *dist[2];			/* Data vectors with along-track distances */
	double *time[2];			/* Data vectors with along-track times (or dummy node indices) */
	double *t, *y;				/* Interpolation y(t) arrays */
	double *out;				/* Output record array */
	double west, east, south, north;	/* Used when projecting data before xover calculations */
	double X2SYS_NaN;			/* Value to write out when result is NaN */
	double lower_speed = 0.0;		/* Ignore crossovers on segments that implies speed lower than this */
	double upper_speed = DBL_MAX;		/* Ignore crossovers on segments that implies speed higher than this */
	double max_gap = DBL_MAX;		/* Ignore crossovers on segments that have longer length than this */
	double xx, yy;				/* Temporary projection variables */


	struct X2SYS_INFO *s;				/* Data format information  */
	struct GMT_XSEGMENT *ylist_A, *ylist_B;		/* y-indices sorted in increasing order */
	struct GMT_XOVER XC;				/* Structure with resulting crossovers */
	struct X2SYS_FILE_INFO data_set[2];		/* File information */

	PFI out_record;			/* Pointer to function that writes the crossover record */
	FILE *fp;

/*----------------------------------END OF VARIBLE DECLARATIONS-----------------------------------------------*/

	argc = 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 */

	sfile = def;
	gmtdefs.interpolant = 0;	/* Linear interpolation default here */

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {

				/* Common parameters */
			
				case 'R':
				case 'J':
				case 'V':
				case 'b':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
				
				case 'A':	/* Get new distance gap */
					max_gap = atof (&argv[i][2]);
					break;
				case 'C':	/* Distance calculation flag */
					if (argv[i][2] == 'c') dist_flag = 0;
					if (argv[i][2] == 'f') dist_flag = 1;
					if (argv[i][2] == 'g') dist_flag = 2;
					break;
				case 'D':
					sfile = &argv[i][2];
					d_given = TRUE;
					break;
				case 'F':
					fflags = &argv[i][2];
					break;
				case 'G':	/* Geographical coordinates, set discontinuity */
					geographic = TRUE;
					geodetic = 0;
					if (argv[i][2] == 'd') geodetic = 2;
					break;
				case 'I':
					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 'L':	/* Get list of files to check */
					fp_list = &argv[i][2];
					break;
				case 'K':	/* Get list of approved filepair combinations to check */
					fp_combo = &argv[i][2];
					combinations = TRUE;
					break;
				case 'M':               /* Multiple line segments */
					GMT_multisegment (&argv[i][2]);
					break;
				case 'N':	/* Get new timegap */
					half_window_width = atoi (&argv[i][2]);
					break;
				case 'O':	/* Compatibility with old xover mode */
					xover_mode = TRUE;
					break;
				case 'S':	/* Speed checks */
					switch (argv[i][2]) {
						case 'L':
						case 'l':	/* Lower cutoff speed */
							lower_speed = atof (&argv[i][3]);
							speed_check = TRUE;
							break;
						case 'U':
						case 'u':	/* Upper cutoff speed */
							upper_speed = atof (&argv[i][3]);
							speed_check = TRUE;
							break;
						default:
							fprintf (stderr, "%s: Syntax Error: -S<l|u><speed>\n", GMT_program);
							error++;
							break;
					}
					break;
				case 'T':
					TAG = &argv[i][2];
					t_given = TRUE;
					break;
				case 'W':	/* Get new timegap */
					max_time_separation = atof (&argv[i][2]);
					if (max_time_separation == 0.0) max_time_separation = DBL_MAX;	/* default is no gap */
					break;
				case 'Q':	/* Specify internal or external only */
					if (argv[i][2] == 'e') internal = FALSE;
					if (argv[i][2] == 'i') external = FALSE;
					break;
				case '2':	/* Return z1, z1 rather than (z1-z1) and 0.5 * (z1 + z2) */
					two_values = TRUE;
					break;
				default:
					error = TRUE;
					break;
			}
		}
	}

	if (argc == 1 || error || GMT_quick) {
		fprintf (stderr, "%s %s - calculate crossovers\n\n", GMT_program, X2SYS_VERSION);
		fprintf (stderr, "usage: x2sys_cross <files> or -L<filelist> [-A<gap>] [-Cc|f|g] [-D<deffile>] [-F<fields>]\n");
		fprintf (stderr, "	[-G[d/g]]  [-Il|a|c] [-J<params>] [-K<combi.lis>] [M[<flag>]] [-N<window>] [-O]\n");
		fprintf (stderr, "	[-Qe|i] [-Sl|u<speed>] [-T<TAG>] [-W<size>] [-V] [-bo[s]] [-2]\n\n");

		if (GMT_quick) exit (EXIT_FAILURE);

		fprintf (stderr, "	<files> is one or more datafiles, OR\n");
		fprintf (stderr,"	-L Supply file with list of filename\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr,"	-A Ignore crossovers on segments with point separation > <gap> [Default is no gap checking]\n");
		fprintf (stderr,"	-C Select procedure for along-track distance calculation:\n");
		fprintf (stderr,"	   c Cartesian distances [Default]\n");
		fprintf (stderr,"	   f Flat Earth distances\n");
		fprintf (stderr,"	   g Great circle distances\n");
		fprintf (stderr, "	-D definition file for this data set [Default is x2sys]\n");
		fprintf (stderr, "	-F is comma-separated list of column names to output [Default are all fields]\n");
		fprintf (stderr,"	-G for geographical coordinates.  Append g for discontinuity at Greenwich (output 0/360 [Default])\n");
		fprintf (stderr,"	   and append d for discontinuity at Dateline (output -180/+180)\n");
		fprintf (stderr, "	-I sets the interpolation mode.  Choose among:\n");
		fprintf (stderr, "	   l Linear interpolation [Default]\n");
		fprintf (stderr, "	   a Akima spline interpolation\n");
		fprintf (stderr, "	   c Acubic spline interpolation\n");
		GMT_explain_option ('J');
		fprintf (stderr,"	-K Gives list of file pairs that are ok to compare [Default is all combinations]\n");
		GMT_explain_option ('M');
		fprintf (stderr,"	-N Maximum points on either side of xover to use in interpolation [Default is 6]\n");
		fprintf (stderr,"	-O Output data using old XOVER format [Default is X2SYS format]\n");
		fprintf (stderr,"	-Q Append e for external crossovers\n");
		fprintf (stderr,"	   Append i for internal crossovers [Default is all crossovers]\n");
		GMT_explain_option ('R');
		fprintf (stderr, "	-S Sets limits on lower and upper speeds:\n");
		fprintf (stderr, "	   -Sl sets lower speed [Default is 0]\n");
		fprintf (stderr, "	   -Sl sets upper speed [Default is Infinity]\n");
		fprintf (stderr, "	-T <TAG> is the system tag for the data set.  If not given, please set -D\n");
		GMT_explain_option ('V');
		fprintf (stderr,"	-W maximum time gap (in user units) allowed at crossover [Default is infinite]\n");
		fprintf (stderr,"	-2 Return z-values for each track [Default is crossover and mean value]\n");
		GMT_explain_option ('o');
		exit (EXIT_FAILURE);
	}

	if (t_given && d_given) {
		fprintf (stderr, "%s: ERROR: Only give one of -D and -T\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (!(internal || external)) {
		fprintf (stderr, "%s: Error: Only one of -Qe -Qi can be specified!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (half_window_width < 1) {
		fprintf (stderr, "%s: Error -N: window must be at least 1\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (max_time_separation < 0.0) {
		fprintf (stderr, "%s: Error -W: maximum gap must be > 0!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (max_gap < 0.0) {
		fprintf (stderr, "%s: Error -A: minimum gap must be >= 0!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (lower_speed > upper_speed) {
		fprintf (stderr, "%s: Error -S: lower speed cutoff higher than upper cutoff!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (t_given)
		x2sys_set_system (TAG, &s, NULL, &GMT_io);
	else
		s = x2sys_initialize (sfile, &GMT_io);	/* Initialize X2SYS and info structure */
	if (GMT_io.multi_segments) {	/* Files have multiple segments; make sure this is also set in s */
		s->multi_segment = TRUE;
		s->ms_flag = GMT_io.EOF_flag;
	}
	if (!geographic) GMT_io.in_col_type[0] = GMT_io.in_col_type[1] = GMT_io.out_col_type[0] = GMT_io.out_col_type[1] = GMT_IS_UNKNOWN;
	if (geographic) s->geodetic = geodetic;	/* Override setting */
	if (fflags) x2sys_pick_fields (fflags, s);

	if (fp_list)	/* Read file names */
	{
		n_legs = x2sys_read_list (fp_list, &file);
	}
	else {	/* Files given on command line */

		n_alloc = GMT_CHUNK;
		file = (char **)GMT_memory (VNULL, n_alloc, sizeof (char *), GMT_program);

		for (i = 1, A = 0; i < argc; i++) {
	
			if (argv[i][0] == '-') continue;	/* Skip options */
			
			file[A] = GMT_memory (VNULL, (size_t)(strlen (argv[i])+1), sizeof (char), GMT_program);
			strcpy (file[A], argv[i]);
			A++;
			if (A == n_alloc) {
				n_alloc += GMT_CHUNK;
				file = (char **)GMT_memory ((void *)file, n_alloc, sizeof (char *), GMT_program);
			}
		}
		n_legs = A;
		file = (char **)GMT_memory ((void *)file, n_legs, sizeof (char *), GMT_program);
	}

	if (n_legs == 0) {
		fprintf (stderr, "%s: Error: Must give at least one data set!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (combinations) {	/* Read list of acceptable file combinations */

		if ((fp = fopen (fp_combo, "r")) == NULL) {
			fprintf (stderr, "%s: Error: Could not open combinations file %s!\n", GMT_program, fp_combo);
			exit (EXIT_FAILURE);
		}

		n_alloc = GMT_CHUNK;
		pair = (struct PAIR *) GMT_memory (VNULL, n_alloc, sizeof (struct PAIR), GMT_program);

		while (fgets (line, BUFSIZ, fp)) {

			if (line[0] == '#' || line[0] == '\n') continue;	/* Skip comments and blanks */

			if (sscanf (line, "%s %s", name1, name2) != 2) {
				fprintf (stderr, "%s: Error: Error decoding combinations file for pair %d!\n", GMT_program, n_pairs);
				exit (EXIT_FAILURE);
			}
			pair[n_pairs].id1 = GMT_memory (VNULL, (size_t)(strlen (name1)+1), sizeof (char), GMT_program);
			strcpy (pair[n_pairs].id1, name1);
			pair[n_pairs].id2 = GMT_memory (VNULL, (size_t)(strlen (name2)+1), sizeof (char), GMT_program);
			strcpy (pair[n_pairs].id2, name2);
			n_pairs++;
			if (n_pairs == n_alloc) {
				n_alloc += GMT_CHUNK;
				pair = (struct PAIR *) GMT_memory ((void *)pair, n_alloc, sizeof (struct PAIR), GMT_program);
			}
		}
		fclose (fp);

		if (!n_pairs) {
			fprintf (stderr, "%s: Error: No combinations found in file %s!\n", GMT_program, fp_combo);
			exit (EXIT_FAILURE);
		}
		pair = (struct PAIR *) GMT_memory ((void *)pair, n_pairs, sizeof (struct PAIR), GMT_program);
	}

	if (xover_mode && !GMT_io.binary[1]) {	/* Backwards compatibility */
		out_record = (PFI) x2sys_xover_output;
		header_format = x2sys_xover_header;
		X2SYS_NaN = (double) GMTMGG_NODATA;
		first_header = FALSE;	/* No special header */
	}
	else {
		out_record = (PFI) GMT_output;
		header_format = x2sys_header;
		X2SYS_NaN = GMT_d_NaN;
	}

	if (gmtdefs.interpolant == 0) half_window_width = 1;
	window_width = 2 * half_window_width;
	n_data_col = x2sys_n_data_cols (s);
	got_time = (s->t_col >= 0);

	n_output = 8 + 2 * n_data_col;
	n_bin_header = (n_output - 1) * ((GMT_io.single_precision[1]) ? sizeof (float) : sizeof (double));
	if (n_data_col == 0){
		xover_locations_only = TRUE;
	}
	else {	/* Set the actual column numbers with data fields */
		t = (double *) GMT_memory (VNULL, (size_t)window_width, sizeof (double), GMT_program);
		y = (double *) GMT_memory (VNULL, (size_t)window_width, sizeof (double), GMT_program);
		col_number = (int *) GMT_memory (VNULL, (size_t)n_data_col, sizeof (int), GMT_program);
		ok = (BOOLEAN *) GMT_memory (VNULL, (size_t)n_data_col, sizeof (BOOLEAN), GMT_program);
		for (col = k = 0; col < s->n_fields; col++) {
			if (col == s->x_col || col == s->y_col || col == s->t_col) continue;
			if (!s->use_column[col]) continue;	/* Turned off by -F */
			col_number[k++] = col;
		}
		col_alloc = n_data_col * sizeof (int);
		if (s->t_col < 0 && gmtdefs.verbose) fprintf (stderr, "%s: No time column, use dummy times\n", GMT_program);
	}

	out = (double *) GMT_memory (VNULL, (size_t)n_output, sizeof (double), GMT_program);
	xdata[0] = (double *) GMT_memory (VNULL, (size_t)s->n_fields, sizeof (double), GMT_program);
	xdata[1] = (double *) GMT_memory (VNULL, (size_t)s->n_fields, sizeof (double), GMT_program);

	if (project_info.region_supplied && project_info.projection >= 0) {
		do_project = TRUE;
		geographic = FALSE;	/* Since we then have x,y projected coordinates, not lon,lat */
		dist_flag = 0;
		GMT_map_setup (west, east, south, north);
	}

	for (A = 0; A < n_legs; A++) {	/* Loop over all files */

		if (s->x_col < 0 || s->x_col < 0) {
			fprintf (stderr, "%s: Error: x and/or y column not found for file %s!\n", GMT_program, file[A]);
			exit (EXIT_FAILURE);
		}

		if ((n_rec[0] = (s->read_file) (file[A], &data[0], s, &data_set[0], &GMT_io)) < 0) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file[A]);
			exit (EXIT_FAILURE);
		}

		if (do_project) {	/* Convert all the coordinates */
			for (i = 0; i < n_rec[0]; i++) {
				GMT_geo_to_xy (data[0][s->x_col][i], data[0][s->y_col][i], &xx, &yy);
				data[0][s->x_col][i] = xx;
				data[0][s->y_col][i] = yy;
			}
		}

		dist[0] = x2sys_distances (data[0][s->x_col], data[0][s->y_col], n_rec[0], dist_flag);

		time[0] = (s->t_col < 0) ? x2sys_dummytimes (n_rec[0]) : data[0][s->t_col];
			
		ylist_A = GMT_init_track (data[0][s->x_col], data[0][s->y_col], n_rec[0]);

		for (B = A; B < n_legs; B++) {

			if (!internal && (A == B)) continue;	/* Only do internal errors */
			if (!external && (A != B)) continue;	/* Only do external errors */

			if (combinations && !combo_ok (file[A], file[B], pair, n_pairs)) continue;	/* Do not want this combo */
			if (A == B) {	/* Just set pointers */
				data[1] = data[0];
				dist[1] = dist[0];
				time[1] = time[0];
				n_rec[1] = n_rec[0];
				ylist_B = ylist_A;
				data_set[1] = data_set[0];
			}
			else {	/* Must read a second file */

				if ((n_rec[1] = (s->read_file) (file[B], &data[1], s, &data_set[1], &GMT_io)) < 0) {
					fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file[B]);
					exit (EXIT_FAILURE);
				}

				if (do_project) {	/* Convert all the coordinates */
					for (i = 0; i < n_rec[0]; i++) {
						GMT_geo_to_xy (data[1][s->x_col][i], data[1][s->y_col][i], &xx, &yy);
						data[1][s->x_col][i] = xx;
						data[1][s->y_col][i] = yy;
					}
				}

				dist[1] = x2sys_distances (data[1][s->x_col], data[1][s->y_col], n_rec[1], dist_flag);

				time[1] = (s->t_col < 0) ? x2sys_dummytimes (n_rec[1]) : data[1][s->t_col];

				ylist_B = GMT_init_track (data[1][s->x_col], data[1][s->y_col], n_rec[1]);
			}

			/* Calculate all possible crossover locations */

			nx = GMT_crossover (data[0][s->x_col], data[0][s->y_col], data_set[0].ms_rec, ylist_A, n_rec[0], data[1][s->x_col], data[1][s->y_col], data_set[1].ms_rec, ylist_B, n_rec[1], (A == B), &XC);

			if (nx && xover_locations_only) {	/* Report crossover locations only */
				if (!GMT_io.binary[1]) fprintf (GMT_stdout, "> %s - %s\n", file[A], file[B]);
				for (i = 0; i < nx; i++) {
					out[0] = XC.x[i];
					out[1] = XC.y[i];
					if (geographic) GMT_lon_range_adjust (geodetic, &out[0]);
					GMT_output (GMT_stdout, 2, out);
				}
				GMT_x_free (&XC);
			}
			else if (nx) {	/* Got crossovers, now estimate crossover values */

				first_crossover = TRUE;

				for (i = 0; i < nx; i++) {	/* For each potential crossover */

					memset ((void *)ok, 0, col_alloc);
					n_ok = 0;

					for (k = 0; k < 2; k++) {	/* For each of the two data sets involved */

						/* Get node number to each side of crossover location */

				/*	--o----------o--------o------X-------o-------o----------o--
							      ^      ^       ^
							    left   xover   right			*/

						left[k]  = (int) floor (XC.xnode[k][i]);
						right[k] = (int) ceil  (XC.xnode[k][i]);

						deld = dist[k][right[k]] - dist[k][left[k]];

						/* Check if we have exceeded the maximum point separation */

						if (deld > max_gap) continue;

						delt = time[k][right[k]] - time[k][left[k]];

						/* Check if speed is outside accepted domain */

						if (speed_check && ((speed = (deld/delt)) < lower_speed || speed > upper_speed)) continue;

						/* Linearly estimate the crossover times and distances */

						dt = XC.xnode[k][i] - left[k];
						time_x[k] = time[k][left[k]];
						dist_x[k] = dist[k][left[k]];
						if (dt > 0.0) {
							time_x[k] += dt * delt;
							dist_x[k] += dt * deld;
						}


						for (j = 0; j < n_data_col; j++) {	/* Evaluate each field at the crossover */

							col = col_number[j];
		
							start = t_right = left[k];
							end = t_left = right[k];
							n_left = n_right = 0;

							xdata[k][col] = GMT_d_NaN;	/* In case of nuthin' */

							/* First find the required <window> points to the left of the xover */

							while (start >= 0 && n_left < half_window_width) {
								if (!GMT_is_dnan (data[k][col][start])) {
									n_left++;
									if (t_left > left[k]) t_left = start;
									y[half_window_width-n_left] = data[k][col][start];
									t[half_window_width-n_left] = time[k][start];
								}
								start--;
							}

							if (!n_left || ((time_x[k] - time[k][t_left]) > max_time_separation)) continue;

							/* Ok, that worked.  Now for the right side: */

							while (end < n_rec[k] && n_right < half_window_width) {
								if (!GMT_is_dnan (data[k][col][end])) {
									y[half_window_width+n_right] = data[k][col][end];
									t[half_window_width+n_right] = time[k][end];
									n_right++;
									if (t_right < right[k]) t_right = end;
								}
								end++;
							}

							if (!n_right || ((delt = (time[k][t_right] - time[k][t_left])) > max_time_separation)) continue;

							/* Ok, got enough data to interpolate at xover */

							first = half_window_width - n_left;
							GMT_intpol (&t[first], &y[first], (n_left + n_right), 1, &time_x[k], &xdata[k][col], gmtdefs.interpolant);
							ok[j]++;
							n_ok++;
						}
					}

					/* Only output crossover if there are any data there */

					if (n_ok == 0) continue;
					for (j = n_ok = 0; j < n_data_col; j++) if (ok[j] == 2) n_ok++;
					if (n_ok == 0) continue;

					/* OK, got something to report */

					/* Load the out array */

					out[0] = XC.x[i];	/* Crossover location */
					out[1] = XC.y[i];

					for (k = 0; k < 2; k++) {	/* Get time, distance, and headings */

						/* Get time */

						out[2+k] = time_x[k];

						/* Get cumulative distance at crossover */

						out[k+4] = dist_x[k];

						/* Estimate heading there */

						j = k + 6;
						delx = data[k][s->x_col][right[k]] - data[k][s->x_col][left[k]];
						dely = data[k][s->y_col][right[k]] - data[k][s->y_col][left[k]];
						if (geographic) delx *= cosd (0.5 * (data[k][s->y_col][right[k]] + data[k][s->y_col][left[k]]));
						if (delx == 0.0)
							out[j] = (dely > 0.0) ? 0.0 : 180.0;
						else {
							out[j] = 90.0 - d_atan2 (dely, delx) * R2D;
							if (out[j] < 0.0) out[j] += 360.0;
						}
					}

					/* Calculate crossover and mean value */

					for (k = 0, j = 8; k < n_data_col; k++) {
						if (two_values) {
							col = col_number[k];
							out[j++] = xdata[0][col];
							out[j++] = xdata[1][col];
						}
						else {
							if (ok[k] == 2) {
								col = col_number[k];
								out[j++] = xdata[0][col] - xdata[1][col];
								out[j++] = 0.5 * (xdata[0][col] + xdata[1][col]);
							}
							else {
								out[j] = out[j+1] = X2SYS_NaN;
								j += 2;
							}
						}
					}

					if (first_header) {	/* Write the header record */
						t_or_i = (got_time) ? 't' : 'i';
						fprintf (GMT_stdout, "#x\ty\t%c0\t%c1\tdist_0\tdist_1\thead_0\thead_1", t_or_i, t_or_i);
						for (j = 0; j < n_data_col; j++) {
							col = col_number[j];
							if (two_values)
								fprintf (GMT_stdout, "\t%s_0\t%s_1", s->info[col].name, s->info[col].name);
							else
								fprintf (GMT_stdout, "\t%s_X\t%s_M", s->info[col].name, s->info[col].name);
						}
						fputc ('\n', GMT_stdout);
						first_header = FALSE;
					}

					if (first_crossover) {
						if (GMT_io.binary[1]) {	/* Segment record has x = NaN as flag */
							(*out_record) (GMT_stdout, 1, &GMT_d_NaN);
							sprintf (line, header_format, file[A], data_set[0].year, file[B], data_set[1].year);
							fwrite ((void *)line, sizeof (char), (size_t)n_bin_header, GMT_stdout);
						}
						else
							fprintf (GMT_stdout, header_format, file[A], data_set[0].year, file[B], data_set[1].year);
						first_crossover = FALSE;
					}

					if (geographic) GMT_lon_range_adjust (geodetic, &out[0]);
					(*out_record) (GMT_stdout, n_output, out);
				}

				GMT_x_free (&XC);
			}
			
			if (A != B) {	/* Must free up memory for B */
				x2sys_free_data (data[1], s->n_fields);
				GMT_free ((void *)dist[1]);
				if (s->t_col < 0) GMT_free ((void *)time[1]);
				GMT_free ((void *)ylist_B);
			}
		}

		/* Must free up memory for A */

		x2sys_free_data (data[0], s->n_fields);
		GMT_free ((void *)dist[0]);
		if (s->t_col < 0) GMT_free ((void *)time[0]);
		GMT_free ((void *)ylist_A);
	}

	/* Free up other arrays */

	GMT_free ((void *)xdata[0]);
	GMT_free ((void *)xdata[1]);
	GMT_free ((void *)out);
	if (n_data_col) {
		GMT_free ((void *)t);
		GMT_free ((void *)y);
		GMT_free ((void *)col_number);
		GMT_free ((void *)ok);
	}

	x2sys_free_info (s);

	GMT_end (argc, argv);
}

BOOLEAN combo_ok (char *name_1, char *name_2, struct PAIR *pair, int n_pairs)
{
	int i;

	/* Return TRUE if this particular combination is found in the list of pairs */

	for (i = 0; i < n_pairs; i++) {
		if (!(strcmp (name_1, pair[i].id1) || strcmp (name_2, pair[i].id2))) return (TRUE);
		if (!(strcmp (name_2, pair[i].id1) || strcmp (name_1, pair[i].id2))) return (TRUE);
	}
	return (FALSE);
}
