/*-----------------------------------------------------------------
 *
 *      Copyright (c) 1999-2024 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
 *      See LICENSE.TXT 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 Lesser General Public License as published by
 *      the Free Software Foundation; version 3 or any later version.
 *
 *      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 Lesser General Public License for more details.
 *
 *      Contact info: www.generic-mapping-tools.org
 *--------------------------------------------------------------------*/
/* 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.
 *		31-MAR-2006: Changed -O to -L to avoid clash with GMT.
 *
 */

#include "gmt_dev.h"
#include "longopt/x2sys_cross_inc.h"
#include "mgd77/mgd77.h"
#include "x2sys.h"

#define THIS_MODULE_CLASSIC_NAME	"x2sys_cross"
#define THIS_MODULE_MODERN_NAME	"x2sys_cross"
#define THIS_MODULE_LIB		"x2sys"
#define THIS_MODULE_PURPOSE	"Calculate crossovers between track data files"
#define THIS_MODULE_KEYS	"<D{,>D}"
#define THIS_MODULE_NEEDS	""
#define THIS_MODULE_OPTIONS "->RVbd"

/* Control structure for x2sys_cross */

#define HHI	0
#define VLO	1
#define	VHI	2

enum x2sys_sets {
	SET_A = 0,
	SET_B = 1
};

enum x2sys_types {
	X2SYS_EXTERNAL = 1,
	X2SYS_INTERNAL = 2,
	X2SYS_BOTH = 3
};

struct X2SYS_CROSS_CTRL {
	struct X2SYS_CROSS_A {	/* -A */
		bool active;
		char *file;
	} A;
	struct X2SYS_CROSS_C {	/* -C */
		bool active;
		char *file;
	} C;
	struct X2SYS_CROSS_D {	/* -Dg|s|n */
		bool active;	/* Force selection if true, else examine */
		int mode;	/* -1 for S pole, +1 for N pole */
	} D;
	struct X2SYS_CROSS_E {	/* -E<limit> */
		bool active;
		double limit;
	} E;
	struct X2SYS_CROSS_I {	/* -Ia|c|l|n */
		bool active;
		int mode;
	} I;
	struct X2SYS_CROSS_S {	/* -S */
		bool active[3];
		double limit[3];
	} S;
	struct X2SYS_CROSS_T {	/* -T */
		bool active;
		char *TAG;
	} T;
	struct X2SYS_CROSS_W {	/* -W */
		bool active;
		unsigned int width;
	} W;
	struct X2SYS_CROSS_Q {	/* -Q */
		bool active;
		int mode;
	} Q;
	struct X2SYS_CROSS_Z {	/* -Z */
		bool active;
	} Z;
};

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

static void *New_Ctrl (struct GMT_CTRL *GMT) {	/* Allocate and initialize a new control structure */
	struct X2SYS_CROSS_CTRL *C;

	C = gmt_M_memory (GMT, NULL, 1, struct X2SYS_CROSS_CTRL);

	/* Initialize values whose defaults are not 0/false/NULL */

	C->S.limit[VHI] = DBL_MAX;	/* Ignore crossovers on segments that implies speed higher than this */
	C->W.width = 3;			/* Number of points on either side in the interpolation */
	return (C);
}

static void Free_Ctrl (struct GMT_CTRL *GMT, struct X2SYS_CROSS_CTRL *C) {	/* Deallocate control structure */
	if (!C) return;
	gmt_M_str_free (C->A.file);
	gmt_M_str_free (C->C.file);
	gmt_M_str_free (C->T.TAG);
	gmt_M_free (GMT, C);
}

static int usage (struct GMTAPI_CTRL *API, int level) {
	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
	GMT_Usage (API, 0, "usage: %s <files> -T<TAG> [-A<pairs>] [-C[<fname>]] [-D[S|N]] [-E<limit>] [-Ia|c|l|n] [-Qe|i] "
		"[%s] [-Sl|h|u<speed>] [%s] [-W<size>] [-Z] [%s] [%s] [%s]\n",
		name, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_do_OPT, GMT_PAR_OPT);

	if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);

	GMT_Usage (API, 1, "Note 1: Output is x y t1 t2 d1 d2 az1 az2 v1 v2 xval1 xmean1 xval2 xmean2 ...");
	GMT_Usage (API, 1, "Note 2: If time is not selected (or present) we use record numbers as proxies i1 i2");

	GMT_Message (API, GMT_TIME_NONE, "\n  REQUIRED ARGUMENTS:\n");
	GMT_Usage (API, 1, "\n<files> is one or more datafiles, or give =<files.lis> for a file with a list of datafiles.");
	GMT_Usage (API, 1, "\n-T<TAG>");
	GMT_Usage (API, -2, "Set the system tag for this compilation.");
	GMT_Message (API, GMT_TIME_NONE, "\n  OPTIONAL ARGUMENTS:\n");
	GMT_Usage (API, 1, "\n-A<pairs>");
	GMT_Usage (API, -2, "Give file with list of track pairs to process [Default processes all combinations].");
	GMT_Usage (API, 1, "\n-C[<fname>]");
	GMT_Usage (API, -2, "Print run time for each pair. Optionally append <fname> to write them to that file.");
	GMT_Usage (API, 1, "\n-D[S|N]");
	GMT_Usage (API, -2, "Control geographic coordinate conversions. By default we automatically convert "
		"lon,lat to polar coordinates if contained within one hemisphere. -D turns this off, while "
		"-DS or -DN forces the conversion using the specified hemisphere [auto-selected].");
	GMT_Usage (API, 1, "\n-E<limit>");
	GMT_Usage (API, -2, "Exclude crossovers from lines with orientation differences less than <limit> [0].");
	GMT_Usage (API, -2, "Select an interpolation mode:");
	GMT_Usage (API, 3, "a: Akima spline interpolation.");
	GMT_Usage (API, 3, "c: Cubic spline interpolation.");
	GMT_Usage (API, 3, "l: Linear interpolation [Default].");
	GMT_Usage (API, 3, "n: Nearest neighbor interpolation.");
	GMT_Usage (API, 1, "\n-Qe|i");
	GMT_Usage (API, -2, "Select a sub-group of crossovers [Default is both]:");
	GMT_Usage (API, 3, "e: Find external crossovers.");
	GMT_Usage (API, 3, "i: Find internal crossovers.");
	GMT_Option (API, "R");
	GMT_Usage (API, 1, "\n-Sl|h|u<speed>");
	GMT_Usage (API, -2, "Set limits on lower and upper speeds (units determined by -Ns):");
	GMT_Usage (API, 3, "l: Set lower speed [Default is 0].");
	GMT_Usage (API, 3, "h: No headings should be computed if velocity drops below this value [0].");
	GMT_Usage (API, 3, "u: Set upper speed [Default is Infinity].");
	GMT_Option (API, "V");
	GMT_Usage (API, 1, "\n-W<size>");
	GMT_Usage (API, -2, "Set maximum points on either side of crossover to use in interpolation [Default is 3].");
	GMT_Usage (API, 1, "\n-Z Return z-values for each track [Default is crossover and mean value].");
	GMT_Option (API, "bo,do,.");

	return (GMT_MODULE_USAGE);
}

static int parse (struct GMT_CTRL *GMT, struct X2SYS_CROSS_CTRL *Ctrl, struct GMT_OPTION *options) {

	/* This parses the options provided to grdcut and sets parameters in CTRL.
	 * Any GMT common options will override values set previously by other commands.
	 * It also replaces any file names specified as input or output with the data ID
	 * returned when registering these sources/destinations with the API.
	 */

	unsigned int n_errors = 0, n_files[2] = {0, 0};
	struct GMT_OPTION *opt = NULL;
	struct GMTAPI_CTRL *API = GMT->parent;

	for (opt = options; opt; opt = opt->next) {	/* Process all the options given */

		switch (opt->option) {
			/* Common parameters */

			case '<':	/* Skip input files since their paths depend on tag */
				n_files[GMT_IN]++;
				break;
			case '>':	/* Got named output file */
				n_files[GMT_OUT]++;
				break;

			/* Processes program-specific parameters */

			case 'A':	/* Get list of approved filepair combinations to check */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
				n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->A.file));
				break;
			case 'C':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
				if (strlen(opt->arg))
					Ctrl->C.file = strdup (opt->arg);
				break;
			case 'D':	/* Determines if projection should happen for geographic coordinates */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
				switch (opt->arg[0]) {
					case 'S':	Ctrl->D.mode = -1; break;	/* Force projection using S pole */
					case 'N':	Ctrl->D.mode = +1; break;	/* Force projection using N pole */
					case '\0':	Ctrl->D.mode =  0; break;	/* No projection  */
				}
				break;
			case 'E':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
				if (opt->arg[0])
					Ctrl->E.limit = atof (opt->arg);
				else {
					GMT_Report (API, GMT_MSG_ERROR, "Option -E: Argument to set limit is required\n");
					n_errors++;
				}
				break;
			case 'I':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
				switch (opt->arg[0]) {
					case 'l':
						Ctrl->I.mode = GMT_SPLINE_LINEAR;
						break;
					case 'a':
						Ctrl->I.mode = GMT_SPLINE_AKIMA;
						break;
					case 'c':
						Ctrl->I.mode = GMT_SPLINE_CUBIC;
						break;
					case 'n':
						Ctrl->I.mode = GMT_SPLINE_NN;
						break;
					default:
						n_errors++;
						break;
				}
				break;
			case 'S':	/* Speed checks */
				switch (opt->arg[0]) {
					case 'L':
					case 'l':	/* Lower cutoff speed */
						Ctrl->S.limit[VLO] = atof (&opt->arg[1]);
						n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active[VLO]);
						break;
					case 'U':
					case 'u':	/* Upper cutoff speed */
						Ctrl->S.limit[VHI] = atof (&opt->arg[1]);
						n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active[VHI]);
						break;
					case 'H':
					case 'h':	/* Heading calculation cutoff speed */
						Ctrl->S.limit[HHI] = atof (&opt->arg[1]);
						n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active[HHI]);
						break;
					default:
						GMT_Report (API, GMT_MSG_ERROR, "Option -S: Syntax is -S<l|h|u><speed>\n");
						n_errors++;
						break;
				}
				break;
			case 'T':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
				n_errors += gmt_get_required_string (GMT, opt->arg, opt->option, 0, &Ctrl->T.TAG);
				break;
			case 'W':	/* Get new window half-width as number of points */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
				n_errors += gmt_get_required_uint (GMT, opt->arg, opt->option, 0, &Ctrl->W.width);
				break;
			case 'Q':	/* Specify internal or external only */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
				if (opt->arg[0] == 'e') Ctrl->Q.mode = X2SYS_EXTERNAL;
				else if (opt->arg[0] == 'i') Ctrl->Q.mode = X2SYS_INTERNAL;
				else Ctrl->Q.mode = X2SYS_BOTH;
				break;
			case 'Z':	/* Return z1, z1 rather than (z1-z1) and 0.5 * (z1 + z2) */
				n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
				n_errors += gmt_get_no_argument (GMT, opt->arg, opt->option, 0);
				break;
			case 'J':
				if (gmt_M_compat_check (GMT, 6)) {
					GMT_Report (API, GMT_MSG_COMPAT, "Option -J is no longer needed or used in x2sys_cross, ignored\n");
					break;
				}
				/* Intentionally fall through - if not compat mode we fall down here on purpose and fail I think */
			default:	/* Report bad options */
				n_errors += gmt_default_option_error (GMT, opt);
				break;
		}
	}

	n_errors += gmt_M_check_condition (GMT, n_files[GMT_IN] == 0, "No track files given\n");
	n_errors += gmt_M_check_condition (GMT, n_files[GMT_OUT] > 1, "More than one output file given\n");
	n_errors += gmt_M_check_condition (GMT, !Ctrl->T.active || !Ctrl->T.TAG, "Option -T must be used to set the TAG\n");
	n_errors += gmt_M_check_condition (GMT, Ctrl->W.width < 1, "Option -W: window must be at least 1\n");
	n_errors += gmt_M_check_condition (GMT, Ctrl->S.limit[VLO] > Ctrl->S.limit[VHI], "Option -S: lower speed cutoff higher than upper cutoff!\n");
	n_errors += gmt_M_check_condition (GMT, Ctrl->Q.mode == X2SYS_BOTH, "Option -Q: Only one of -Qe -Qi can be specified!\n");


	return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}

GMT_LOCAL int x2syscross_combo_ok (char *name_1, char *name_2, struct X2SYS_CROSS_PAIR *pair, uint64_t n_pairs) {
	uint64_t 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);
}

GMT_LOCAL void x2syscross_free_pairs (struct X2SYS_CROSS_PAIR *pair, uint64_t n_pairs) {
	/* Free the strings in the array of pairs */
	uint64_t k;
	for (k = 0; k < n_pairs; k++) {
		gmt_M_str_free (pair[k].id1);
		gmt_M_str_free (pair[k].id2);
	}
}

GMT_LOCAL bool x2syscross_is_outside_region (struct GMT_CTRL *GMT, double lon, double lat, bool geo) {
	if (!GMT->common.R.active[RSET]) return false;	/* Requires -R to be outside */
	if (lat < GMT->common.R.wesn[YLO] || lat > GMT->common.R.wesn[YHI]) return true;	/* Outside y range */
	if (geo) {	/* Must wind the longitude properly */
		if (lon > GMT->common.R.wesn[XHI])    lon -= 360.0;
		while (lon < GMT->common.R.wesn[XLO]) lon += 360.0;
		if (lon > GMT->common.R.wesn[XHI]) return true;	/* Outside longitude range */
	}
	else if (lon < GMT->common.R.wesn[XLO] || lon > GMT->common.R.wesn[XHI]) return true;	/* Outside Cartesian x range */
	return false;	/* Inside */
}

GMT_LOCAL void x2syscross_local_geo_to_xy (double *lon, double *lat, double plat) {
	/* Project lon,lat to polar coordinates (r,theta) then Cartesian x,y */
	double r, s, c;
	r = fabs (*lat - plat);		/* Distance from the pole plat in degrees */
	sincosd (*lon, &s, &c);		/* Longitude is our theta */
	*lon = r * c;	*lat = r * s;	/* Cartesian coordinates */
}

GMT_LOCAL void x2syscross_local_xy_to_geo (double *x, double *y, int plat) {
	/* Project Cartesian x,y to polar coordinates (r,theta) then lon,lat */
	double r, lon, lat;
	r = hypot (*x, *y);				/* Distance from our pole in degrees */
	lat = (plat == 1) ? 90.0 - r : r - 90.0;	/* Latitude depends on the which pole, N or S */
	lon = d_atan2d (*y, *x);			/* Determine theta which is our longitude */
	*x = lon;	*y = lat;			/* Pass back the geographic coordinates */
}

#define bailout(code) {gmt_M_free_options (mode); return (code);}
#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
#define Crashout(code) {gmt_M_free (GMT, duplicate); x2syscross_free_pairs (pair, n_pairs); gmt_M_free (GMT, pair); x2sys_free_list (GMT, trk_name, n_tracks); if (fpC) fclose (fpC); x2sys_end (GMT, s); Return (error);}

EXTERN_MSC int GMT_x2sys_cross (void *V_API, int mode, void *args) {
	char **trk_name = NULL;			/* Name of tracks */
	char line[GMT_BUFSIZ] = {""};		/* buffer */
	char item[GMT_BUFSIZ] = {""};		/* buffer */
	char t_or_i;				/* t = time, i = dummy node time */
	char name1[80] = {""}, name2[80] = {""};		/* Name of two files to be examined */
	char *x2sys_header = "%s %d %s %d %s";

	uint64_t n_rec[2];			/* Number of data records for both files */
	uint64_t window_width;			/* Max number of points to use in the interpolation */
	uint64_t n_tracks = 0;			/* Total number of data sets to compare */
	uint64_t nx;				/* Number of crossovers found for this pair */
	uint64_t *col_number = NULL;		/* Array with the column numbers of the data fields */
	unsigned int n_output;			/* Number of columns on output */
	uint64_t n_pairs = 0;			/* Number of acceptable combinations */
	uint64_t A, B, i, j, col, k, start, n_bad;	/* Misc. counters and local variables */
	uint64_t end, first, n_ok;
	uint64_t n_data_col, left[2], t_left;
	uint64_t n_left, right[2], t_right, n_right;
	uint64_t n_duplicates, n_errors;
	uint64_t add_chunk;
	int scol;
	int error = 0;				/* nonzero for invalid arguments */
	int iplat = 0;				/* Sign of pole for reprojections */
	unsigned int *ok = NULL;

	bool xover_locations_only = false;	/* true if only x,y (and possible indices) to be output */
	bool internal = true;			/* false if only external xovers are needed */
	bool external = true;			/* false if only internal xovers are needed */
	bool do_project = false;		/* true if we must project first */
	bool do_examine = false;		/* true if we should examine y-range to pick closest pole*/
	bool is_geographic = false;		/* true if we must project to polar Cartesian */
	bool got_time = false;			/* true if there is a time column */
	bool first_header = true;		/* true for very first crossover */
	bool first_crossover;			/* true for first crossover between two data sets */
	bool same = false;			/* true when the two tracks we compare have the same name */
	bool has_time[2];			/* true for each tracks that actually has a time column */
	bool *duplicate = NULL;			/* Array, true for any track that is already listed */
	bool cmdline_files = false;		/* true if files where given directly on the command line */

	size_t n_alloc = 1;

	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 deld, delt;			/* Differences in dist and time across xpoint */
	double speed[2];			/* speed across the xpoint ( = deld/delt) */
	double **data[2] = {NULL, NULL};	/* Data matrices for the two data sets to be checked */
	double *xdata[2] = {NULL, NULL};	/* Data vectors with estimated values at crossover points */
	double *dist[2] = {NULL, NULL};		/* Data vectors with along-track distances */
	double *time[2] = {NULL, NULL};		/* Data vectors with along-track times (or dummy node indices) */
	double *t = NULL, *y = NULL;		/* Interpolation y(t) arrays */
	double *out = NULL;			/* Output record array */
	double X2SYS_NaN;			/* Value to write out when result is NaN */
	double dist_scale;			/* Scale to give selected distance units */
	double vel_scale;			/* Scale to give selected velocity units */
	double t_scale;				/* Scale to give time in seconds */
	double plat[2] = {0.0, 0.0};		/* Pole latitude for polar reprojections */
	double ymin[2] = {0.0, 0.0}, ymax[2] = {0.0, 0.0};	/* Latitude range of each file */
	double delta_orientation;	/* Angle between two intersecting tracks */

	clock_t tic = 0, toc = 0;

	struct X2SYS_INFO *s = NULL;			/* Data format information  */
	struct GMT_XSEGMENT *ylist[2] = {NULL, NULL};	/* y-indices sorted in increasing order */
	struct GMT_XOVER XC;				/* Structure with resulting crossovers */
	struct X2SYS_FILE_INFO data_set[2];		/* File information */
	struct X2SYS_BIX Bix;
	struct X2SYS_CROSS_PAIR *pair = NULL;			/* Used with -Akombinations.lis option */
	FILE *fp = NULL, *fpC = NULL;
	struct GMT_RECORD *Out = NULL;
	struct X2SYS_CROSS_CTRL *Ctrl = NULL;
	struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
	struct GMT_OPTION *options = NULL;
	struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API);	/* Cast from void to GMTAPI_CTRL pointer */

/*----------------------------------END OF VARIABLE DECLARATIONS---------------------------------------------*/

	/*----------------------- Standard module initialization and parsing ----------------------*/

	if (API == NULL) return (GMT_NOT_A_SESSION);
	if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE));	/* Return the purpose of program */
	options = GMT_Create_Options (API, mode, args);	if (API->error) return (API->error);	/* Set or get option list */

	if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error);	/* Give usage if requested */

	/* Parse the command-line arguments */

	if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, module_kw, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
	if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
	Ctrl = New_Ctrl (GMT);	/* Allocate and initialize a new control structure */
	if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);

	/*---------------------------- This is the x2sys_cross main code ----------------------------*/

	if (x2sys_err_fail (GMT, x2sys_set_system (GMT, Ctrl->T.TAG, &s, &Bix, &GMT->current.io), Ctrl->T.TAG))
		Return (GMT_RUNTIME_ERROR);
	if (!s->geographic) {
		gmt_set_column_type (GMT, GMT_IO, GMT_X, GMT_IS_UNKNOWN);
		gmt_set_column_type (GMT, GMT_IO, GMT_Y, GMT_IS_UNKNOWN);
	}

	if (s->x_col == -1 || s->y_col == -1) {
		GMT_Report (API, GMT_MSG_ERROR, "lon,lat or x,y are not among data columns!\n");
		x2sys_end (GMT, s);
		Return (GMT_RUNTIME_ERROR);
	}

	if ((error = x2sys_get_tracknames (GMT, options, &trk_name, &cmdline_files)) == 0) {
		GMT_Report (API, GMT_MSG_ERROR, "Must give at least one data set!\n");
		x2sys_end (GMT, s);
		Return (GMT_RUNTIME_ERROR);
	}
	n_tracks = (uint64_t)error;

	GMT->current.setting.interpolant = Ctrl->I.mode;
	if (Ctrl->Q.active) {
		if (Ctrl->Q.mode == X2SYS_EXTERNAL) internal = false;
		if (Ctrl->Q.mode == X2SYS_INTERNAL) external = false;
	}

	GMT_Report (API, GMT_MSG_INFORMATION, "Files found: %" PRIu64 "\n", n_tracks);

	duplicate = gmt_M_memory (GMT, NULL, n_tracks, bool);

	GMT_Report (API, GMT_MSG_INFORMATION, "Checking for duplicates : ");
	/* Make sure there are no duplicates */
	for (A = n_duplicates = 0; A < n_tracks; A++) {	/* Loop over all files */
		if (duplicate[A]) continue;
		for (B = A + 1; B < n_tracks; B++) {
			if (duplicate[B]) continue;
			same = !strcmp (trk_name[A], trk_name[B]);
			if (same) {
				GMT_Report (API, GMT_MSG_ERROR, "File %s repeated on command line - skipped\n", trk_name[A]);
				duplicate[B] = true;
				n_duplicates++;
			}
		}
	}
	GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " found\n", n_duplicates);

	if (Ctrl->A.active) {	/* Read list of acceptable trk_name combinations */

		GMT_Report (API, GMT_MSG_INFORMATION, "Explicit combinations found: ");
		if ((fp = fopen (Ctrl->A.file, "r")) == NULL) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not open combinations file %s!\n", Ctrl->A.file);
			Crashout (GMT_ERROR_ON_FOPEN);
		}

		n_alloc = add_chunk = GMT_CHUNK;
		pair = gmt_M_memory (GMT, NULL, n_alloc, struct X2SYS_CROSS_PAIR);

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

			if (line[0] == '#' || line[0] == '\n') continue;	/* Skip comments and blanks */
			gmt_chop (line);	/* Get rid of CR, LF stuff */

			if (sscanf (line, "%s %s", name1, name2) != 2) {
				GMT_Report (API, GMT_MSG_ERROR, "Unable to decode combinations file for pair %" PRIu64 "!\n", n_pairs);
				fclose (fp);
				Crashout (GMT_RUNTIME_ERROR);
			}
			pair[n_pairs].id1 = strdup (name1);
			pair[n_pairs].id2 = strdup (name2);
			n_pairs++;
			if (n_pairs == n_alloc) {
				size_t old_n_alloc = n_alloc;
				add_chunk *= 2;
				n_alloc += add_chunk;
				pair = gmt_M_memory (GMT, pair, n_alloc, struct X2SYS_CROSS_PAIR);
				gmt_M_memset (&(pair[old_n_alloc]), n_alloc - old_n_alloc, struct X2SYS_CROSS_PAIR);
			}
		}
		fclose (fp);

		if (!n_pairs) {
			GMT_Report (API, GMT_MSG_ERROR, "No combinations found in file %s!\n", Ctrl->A.file);
			Crashout (GMT_RUNTIME_ERROR);
		}
		if (n_pairs < n_alloc) pair = gmt_M_memory (GMT, pair, n_pairs, struct X2SYS_CROSS_PAIR);
		GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 "\n", n_pairs);
	}

	if (Ctrl->C.file) {	/* Open file to store the per pair run time */
		if ((fpC = fopen (Ctrl->C.file, "w")) == NULL) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not open save times file %s!\n", Ctrl->C.file);
			gmt_M_str_free (Ctrl->C.file);
		}
	}

	X2SYS_NaN = GMT->session.d_NaN;

	if (GMT->current.setting.interpolant == GMT_SPLINE_LINEAR) Ctrl->W.width = 1;
	window_width = 2 * Ctrl->W.width;
	n_data_col = x2sys_n_data_cols (GMT, s);
	got_time = (s->t_col >= 0);
	if (!got_time) Ctrl->S.active[VLO] = false;	/* Cannot check speed if there is no time */

	n_output = (unsigned int)(10 + 2 * n_data_col);
	gmt_set_column_type (GMT, GMT_OUT, GMT_X, (!strcmp (s->info[s->x_col].name, "lon")) ? GMT_IS_LON : GMT_IS_FLOAT);
	gmt_set_column_type (GMT, GMT_OUT, GMT_Y, (!strcmp (s->info[s->y_col].name, "lat")) ? GMT_IS_LAT : GMT_IS_FLOAT);
	gmt_set_column_type (GMT, GMT_OUT, GMT_Z, (got_time) ? GMT_IS_ABSTIME : GMT_IS_FLOAT);
	gmt_set_column_type (GMT, GMT_OUT, 3, (got_time) ? GMT_IS_ABSTIME : GMT_IS_FLOAT);

	for (i = 0; i < n_data_col+2; i++) {
		gmt_set_column_type (GMT, GMT_OUT, 4+2*(unsigned int)i, GMT_IS_FLOAT);
		gmt_set_column_type (GMT, GMT_OUT, 5+2*(unsigned int)i, GMT_IS_FLOAT);
	}

	if (n_data_col == 0) {
		xover_locations_only = true;
		n_output = 2;
	}
	else {	/* Set the actual column numbers with data fields */
		t = gmt_M_memory (GMT, NULL, window_width, double);
		y = gmt_M_memory (GMT, NULL, window_width, double);
		col_number = gmt_M_memory (GMT, NULL, n_data_col, uint64_t);
		ok = gmt_M_memory (GMT, NULL, n_data_col, unsigned int);
		for (col = k = scol = 0; col < s->n_out_columns; col++, scol++) {
			if (scol == s->x_col || scol == s->y_col || scol == s->t_col) continue;
			col_number[k++] = col;
		}
		if (s->t_col < 0) GMT_Report (API, GMT_MSG_WARNING, "No time column, use dummy times\n");
	}

	out = gmt_M_memory (GMT, NULL, n_output, double);
	xdata[SET_A] = gmt_M_memory (GMT, NULL, s->n_out_columns, double);
	xdata[SET_B] = gmt_M_memory (GMT, NULL, s->n_out_columns, double);
	Out = gmt_new_record (GMT, out, NULL);	/* Since we only need to worry about numerics in this module */

	gmt_set_segmentheader (GMT, GMT_OUT, true);	/* Turn on segment headers on output */
	gmt_set_tableheader (GMT, GMT_OUT, true);	/* Turn on -ho explicitly */

	/* PW: Cannot use gmt_proj_setup since all that stuff assumes lon,lat are now in cols
	 * 0 and 1.  But that is not true in x2sys since it honors the *.def file layout.
	 * Since gmt_crossover deals with periodic longitudes then we really do not need to
	 * project to Cartesian anyway.  If -DS|N is given (or no -D) and geographic then we internally
	 * project via a r-theta schene to x,y get the crossings and project back afterwards.
	 * We examine the latitude range to select the most suitable pole for the projection (unless set in -D),
	 * This may have issues if data cover both hemispheres to high latitudes. The projection we use is
	 *
	 * (lon,lat) --> (x,y), via x = r * cos(lon), y = r * sin(lon), with r = distance to the pole,
	 *	so r is either 90 - lat for N pole or lat - 90 for S pole.  After computing crossovers we
	 * recover the geographic coordinates via lon = atan2d (y, x) and r = hypot (x,y); the latitude
	 * is then given as 90 - r or r - 90 depending on pole.
	 */

	do_examine = do_project = !(Ctrl->D.active || !s->geographic);	/* If -D is given then no examination, else examine if geographic */
	is_geographic = s->geographic;
	if (s->geographic && Ctrl->D.active && Ctrl->D.mode) {	/* Must project if geographic and user gave -Ds or -Dn */
		do_project = true;		/* Convert lon,lat to a cylindrical, polar projection */
		is_geographic = false;		/* Once we project we have Cartesian data */
		iplat = Ctrl->D.mode;		/* User specified which pole to project to */
		plat[SET_A] = plat[SET_B] = iplat * 90.0;	/* Corresponding latitude of pole */
		GMT_Report (API, GMT_MSG_INFORMATION, "Based on -D setting we will polar project all data using pole latitude %g\n", plat[SET_A]);
	}
	if (gmt_init_distaz (GMT, s->dist_flag ? GMT_MAP_DIST_UNIT : 'X', s->dist_flag, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE)
		Crashout (GMT_NOT_A_VALID_TYPE);

	MGD77_Set_Unit (GMT, s->unit[X2SYS_DIST_SELECTION], &dist_scale, -1);	/* Gets scale which multiplies meters to chosen distance unit */
	MGD77_Set_Unit (GMT, s->unit[X2SYS_SPEED_SELECTION], &vel_scale, -1);	/* Sets output scale for distances using in velocities */
	switch (s->unit[X2SYS_SPEED_SELECTION][0]) {
		case 'c':
			vel_scale = 1.0;
			break;
		case 'e':
			vel_scale /= dist_scale;			/* Must counteract any distance scaling to get meters. dt is in sec so we get m/s */
			break;
		case 'f':
			vel_scale /= (METERS_IN_A_FOOT * dist_scale);	/* Must counteract any distance scaling to get feet. dt is in sec so we get ft/s */
			break;
		case 'k':
			vel_scale *= (3600.0 / dist_scale);		/* Must counteract any distance scaling to get km. dt is in sec so 3600 gives km/hr */
			break;
		case 'm':
			vel_scale *= (3600.0 / dist_scale);		/* Must counteract any distance scaling to get miles. dt is in sec so 3600 gives miles/hr */
			break;
		case 'n':
			vel_scale *= (3600.0 / dist_scale);		/* Must counteract any distance scaling to get miles. dt is in sec so 3600 gives miles/hr */
			break;
		case 'u':
			vel_scale /= (METERS_IN_A_SURVEY_FOOT * dist_scale);	/* Must counteract any distance scaling to get survey feet. dt is in sec so we get ft/s */
			break;
		default:	/*Cartesian */
			break;
	}
	t_scale = GMT->current.setting.time_system.scale;	/* Convert user's TIME_UNIT to seconds */

	if ((error = GMT_Set_Columns (API, GMT_OUT, n_output, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
		Crashout (error);
	}
	if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {	/* Registers default output destination, unless already set */
		Crashout (error);
	}
	if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) {	/* Enables data output and sets access mode */
		Crashout (error);
	}
	if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) {	/* Sets output geometry */
		Crashout (error);
	}

	for (A = 0; A < n_tracks; A++) {	/* Loop over all files */
		if (duplicate[A]) continue;

		if (s->x_col < 0 || s->y_col < 0) {
			GMT_Report (API, GMT_MSG_ERROR, "x and/or y column not found for track %s!\n", trk_name[A]);
			Crashout (error);
		}

		if (x2sys_err_fail (GMT, (s->read_file) (GMT, trk_name[A], &data[SET_A], s, &data_set[SET_A], &GMT->current.io, &n_rec[SET_A]), trk_name[A]))
			Return (GMT_RUNTIME_ERROR);

		if (n_rec[SET_A] == 0) {	/* No data in track A */
			x2sys_free_data (GMT, data[SET_A], s->n_out_columns, &data_set[SET_A]);
			continue;
		}

		has_time[SET_A] = false;
		if (got_time) {	/* Check to make sure we do in fact have time */
			for (i = n_bad = 0; i < n_rec[SET_A]; i++) n_bad += gmt_M_is_dnan (data[SET_A][s->t_col][i]);
			if (n_bad < n_rec[SET_A]) has_time[SET_A] = true;
		}

		if ((dist[SET_A] = gmt_dist_array_2 (GMT, data[SET_A][s->x_col], data[SET_A][s->y_col], n_rec[SET_A], dist_scale, s->dist_flag)) == NULL) {
			error = gmt_M_err_fail (GMT, GMT_MAP_BAD_DIST_FLAG, "");
			Crashout (error);
		}

		if (do_examine) {	/* Check all the coordinates and find suitable pole */
			ymin[SET_A] = ymax[SET_A] = data[SET_A][s->y_col][0];
			for (i = 1; i < n_rec[SET_A]; i++) {
				if (data[SET_A][s->y_col][i] < ymin[SET_A]) ymin[SET_A] = data[SET_A][s->y_col][i];
				if (data[SET_A][s->y_col][i] > ymax[SET_A]) ymax[SET_A] = data[SET_A][s->y_col][i];
			}
			plat[SET_A] = ((90.0 - ymax[SET_A]) <= (ymin[SET_A] + 90.0)) ? 90.0 : -90.0;	/* Pick closest pole */
			iplat = (plat[SET_A] > 0.0) ? 1 : -1;
		}
		if (do_project) {	/* Convert coordinates */
			for (i = 0; i < n_rec[SET_A]; i++)
				x2syscross_local_geo_to_xy (&data[SET_A][s->x_col][i], &data[SET_A][s->y_col][i], plat[SET_A]);
		}

		time[SET_A] = (has_time[SET_A]) ? data[SET_A][s->t_col] : x2sys_dummytimes (GMT, n_rec[SET_A]) ;

		gmt_init_track (GMT, data[SET_A][s->y_col], n_rec[SET_A], &ylist[SET_A]);

		for (B = A; B < n_tracks; B++) {
			if (duplicate[B]) continue;

			same = !strcmp (trk_name[A], trk_name[B]);
			if (same && !(A == B)) {
				GMT_Report (API, GMT_MSG_ERROR, "File %s repeated on command line - skipped\n", trk_name[A]);
				continue;
			}
			if (!internal &&  same) continue;	/* Only do external errors */
			if (!external && !same) continue;	/* Only do internal errors */

			if (Ctrl->A.active && !x2syscross_combo_ok (trk_name[A], trk_name[B], pair, n_pairs)) continue;	/* Do not want this combo */

			if (Ctrl->C.active) tic = clock();	/* To report execution time from this pair */

			if (same) {	/* Just set pointers */
				data[SET_B] = data[SET_A];
				dist[SET_B] = dist[SET_A];
				time[SET_B] = time[SET_A];
				has_time[SET_B] = has_time[SET_A];
				n_rec[SET_B] = n_rec[SET_A];
				ylist[SET_B] = ylist[SET_A];
				data_set[SET_B] = data_set[SET_A];
			}
			else {	/* Must read a second file */

				if (x2sys_err_fail (GMT, (s->read_file) (GMT, trk_name[B], &data[SET_B], s, &data_set[SET_B], &GMT->current.io, &n_rec[SET_B]), trk_name[B]))
					Return (GMT_RUNTIME_ERROR);

				if (n_rec[SET_B] == 0) {	/* No data in track B */
					x2sys_free_data (GMT, data[SET_B], s->n_out_columns, &data_set[SET_B]);
					continue;
				}
				has_time[SET_B] = false;
				if (got_time) {	/* Check to make sure we do in fact have time */
					for (i = n_bad = 0; i < n_rec[SET_B]; i++) n_bad += gmt_M_is_dnan (data[SET_B][s->t_col][i]);
					if (n_bad < n_rec[SET_B]) has_time[SET_B] = true;
				}

				if ((dist[SET_B] = gmt_dist_array_2 (GMT, data[SET_B][s->x_col], data[SET_B][s->y_col], n_rec[SET_B], dist_scale, s->dist_flag)) == NULL) {
					error = gmt_M_err_fail (GMT, GMT_MAP_BAD_DIST_FLAG, "");
					Crashout (error);
				}

				if (do_examine) {	/* Check the coordinates and find suitable pole */
					ymin[SET_B] = ymax[SET_B] = data[SET_B][s->y_col][0];
					for (i = 1; i < n_rec[SET_B]; i++) {
						if (data[SET_B][s->y_col][i] < ymin[SET_B]) ymin[SET_B] = data[SET_B][s->y_col][i];
						if (data[SET_B][s->y_col][i] > ymax[SET_B]) ymax[SET_B] = data[SET_B][s->y_col][i];
					}
					plat[SET_B] = ((90.0 - ymax[SET_B]) <= (ymin[SET_B] + 90.0)) ? 90.0 : -90.0;	/* Pick closest pole for dataset B */
					if ((plat[SET_A] * plat[SET_B]) < 0.0) {	/* Auto-determined to be on separate hemispheres; check if this could be problematic */
						/* Arbitrary consider data more than 45 from the Equator to possibly be at high latitude */
						fprintf (stderr, "A (%d): %g/%g  B (%d): %g/%g\n", (int)A, ymin[SET_A], ymax[SET_A], (int)B, ymin[SET_B], ymax[SET_B]);
						if (MAX (fabs (ymin[SET_A]), fabs (ymax[SET_A])) < 45.0 && MAX (fabs (ymin[SET_B]), fabs (ymax[SET_B])) < 45.0)
							GMT_Report (API, GMT_MSG_WARNING, "The two files occupy opposite hemispheres at low latitudes - we select use pole at %g (see -D for discussion)\n", plat[SET_A]);
						else
							GMT_Report (API, GMT_MSG_WARNING, "The two files occupy opposite hemispheres at high latitudes - we select use pole at %g (see -D for discussion)\n", plat[SET_A]);
					}
				}
				if (do_project) {	/* Convert coordinates the same way as A */
					for (i = 0; i < n_rec[SET_B]; i++)
						x2syscross_local_geo_to_xy (&data[SET_B][s->x_col][i], &data[SET_B][s->y_col][i], plat[SET_A]);
				}

				time[SET_B] = (has_time[SET_B]) ? data[SET_B][s->t_col] : x2sys_dummytimes (GMT, n_rec[SET_B]);

				gmt_init_track (GMT, data[SET_B][s->y_col], n_rec[SET_B], &ylist[SET_B]);
			}

			/* Calculate all possible crossover locations */

			nx = gmt_crossover (GMT, data[SET_A][s->x_col], data[SET_A][s->y_col], data_set[SET_A].ms_rec, ylist[SET_A], n_rec[SET_A], data[SET_B][s->x_col], data[SET_B][s->y_col], data_set[SET_B].ms_rec, ylist[SET_B], n_rec[SET_B], (A == B), is_geographic, &XC);

			if (nx && xover_locations_only) {	/* Report crossover locations only */
				double az[2];
				first_crossover = true;
				for (i = 0; i < nx; i++) {
					if (Ctrl->E.active) {
						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-- ----> time
							                          ^    ^       ^
							                        left xover   right			*/

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

							if (left[k] == right[k]) {	/* Crosses exactly on a node; move left or right so interpolation will work */
								if (left[k] > 0)
									left[k]--;	/* Move back so cross occurs at right[k] */
								else
									right[k]++;	/* Move forward so cross occurs at left[k] */
							}
							az[k] = (*GMT->current.map.azimuth_func) (GMT, data[k][s->x_col][right[k]], data[k][s->y_col][right[k]], data[k][s->x_col][left[k]], data[k][s->y_col][left[k]], false);
						}
						/* Ensure azimuths are positive (adding 360) then double to get orientations and then delete by 2 and subtract to get final absolute difference */
						delta_orientation = fabs (fmod (2.0 * (az[0] + 360.0), 360.0) / 2.0 - fmod (2.0 * (az[1] + 360.0), 360.0) / 2.0);
						if (delta_orientation < Ctrl->E.limit)
							continue;	/* Skip this crossover */
					}
					if (first_crossover) {
						sprintf (line, "%s - %s", trk_name[A], trk_name[B]);
						GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, line);
						first_crossover = false;
					}
					out[GMT_X] = XC.x[i];
					out[GMT_Y] = XC.y[i];
					if (do_project)
						x2syscross_local_xy_to_geo (&out[GMT_X], &out[GMT_Y], iplat);
					if (x2syscross_is_outside_region (GMT, out[GMT_X], out[GMT_Y], s->geographic))
						continue;	/* Outside our area of interest */
					if (s->geographic) gmt_lon_range_adjust (s->geodetic, &out[GMT_X]);
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);	/* Write this to output */
				}
				gmt_x_free (GMT, &XC);
			}
			else if (nx) {	/* Got crossovers, now estimate crossover values */
				int64_t start_s;
				first_crossover = true;

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

					out[GMT_X] = XC.x[i];	/* Crossover location */
					out[GMT_Y] = XC.y[i];
					if (do_project)
						x2syscross_local_xy_to_geo (&out[GMT_X], &out[GMT_Y], iplat);
					if (x2syscross_is_outside_region (GMT, out[GMT_X], out[GMT_Y], s->geographic))
						continue;	/* Outside our area of interest */

					gmt_M_memset (ok, n_data_col, unsigned int);
					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-- ----> time
				                          ^    ^       ^
				                        left xover   right			*/

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

						if (left[k] == right[k]) {	/* Crosses exactly on a node; move left or right so interpolation will work */
							if (left[k] > 0)
								left[k]--;	/* Move back so cross occurs at right[k] */
							else
								right[k]++;	/* Move forward so cross occurs at left[k] */
						}

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

						/* Check if speed is outside accepted domain; units were set via x2sys_init -Ns */

						speed[k] = (delt == 0.0) ? GMT->session.d_NaN : vel_scale * (deld / (delt * t_scale));
						if (Ctrl->S.active[VLO] && !gmt_M_is_dnan (speed[k]) && (speed[k] < Ctrl->S.limit[VLO] || speed[k] > Ctrl->S.limit[VHI])) continue;

						/* Linearly estimate the crossover times and distances using fractional point counter */

						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->session.d_NaN;	/* In case of nuthin' */

							/* First find the required <window> points to the left of the xover */
							start_s = start;
							while (start_s >= 0 && n_left < Ctrl->W.width) {
								if (!gmt_M_is_dnan (data[k][col][start])) {
									n_left++;
									if (t_left > left[k]) t_left = start;
									y[Ctrl->W.width-n_left] = data[k][col][start];
									t[Ctrl->W.width-n_left] = time[k][start];
								}
								start--;
								start_s--;
							}

							if (!n_left) continue;
							if (got_time && ((time_x[k] - time[k][t_left]) > Bix.time_gap)) continue;
							if ((dist_x[k] - dist[k][t_left]) > Bix.dist_gap) continue;

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

							while (end < n_rec[k] && n_right < Ctrl->W.width) {
								if (!gmt_M_is_dnan (data[k][col][end])) {
									y[Ctrl->W.width+n_right] = data[k][col][end];
									t[Ctrl->W.width+n_right] = time[k][end];
									n_right++;
									if (t_right < right[k]) t_right = end;
								}
								end++;
							}

							if (!n_right) continue;
							/* See if we pass any gap criteria */
							if (got_time && ((time[k][t_right] - time_x[k]) > Bix.time_gap)) continue;	/* Exceeded time gap */
							if ((dist[k][t_right] - dist_x[k]) > Bix.dist_gap) continue;			/* Exceeded distance gap */

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

							first = Ctrl->W.width - n_left;
							n_errors = gmt_intpol (GMT, &t[first], &y[first], NULL, (n_left + n_right), 1, &time_x[k], &xdata[k][col], 0.0, GMT->current.setting.interpolant);
							if (n_errors == 0) {	/* OK */
								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 */

					for (k = 0; k < 2; k++) {	/* Get times, distances, headings, and velocities */

						/* Get time */

						out[2+k] = (got_time && !has_time[k]) ? X2SYS_NaN : time_x[k];

						/* Get cumulative distance at crossover */

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

						/* Estimate heading there */

						j = k + 6;
						out[j] = (!gmt_M_is_dnan (speed[k]) && (!Ctrl->S.active[HHI] || speed[k] > Ctrl->S.limit[HHI])) ? (*GMT->current.map.azimuth_func) (GMT, data[k][s->x_col][right[k]], data[k][s->y_col][right[k]], data[k][s->x_col][left[k]], data[k][s->y_col][left[k]], false) : X2SYS_NaN;

						/* Estimate velocities there */

						j = k + 8;
						out[j] = (has_time[k]) ? speed[k] : X2SYS_NaN;
					}

					if (Ctrl->E.active) {	/* Avoid grazing crossovers */
						/* Ensure azimuths are positive (adding 360) then double to get orientations and then delete by 2 and subtract to get final absolute difference */
						double delta_orientation = fabs (fmod (2.0 * (out[6] + 360.0), 360.0) / 2.0 - fmod (2.0 * (out[7] + 360.0), 360.0) / 2.0);
						if (delta_orientation < Ctrl->E.limit)
							continue;	/* Skip this crossover */
					}

					/* Calculate crossover and mean value */

					for (k = 0, j = 10; k < n_data_col; k++) {
						if (Ctrl->Z.active) {
							col = col_number[k];
							out[j++] = xdata[SET_A][col];
							out[j++] = xdata[SET_B][col];
						}
						else {
							if (ok[k] == 2) {
								col = col_number[k];
								out[j++] = xdata[SET_A][col] - xdata[SET_B][col];
								out[j++] = 0.5 * (xdata[SET_A][col] + xdata[SET_B][col]);
							}
							else {
								out[j] = out[j+1] = X2SYS_NaN;
								j += 2;
							}
						}
					}

					if (first_header) {	/* Write the header record */
						char *cmd = NULL, *c = GMT->current.setting.io_col_separator;
						t_or_i = (got_time) ? 't' : 'i';
						sprintf (line, "Tag: %s", Ctrl->T.TAG);
						GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, line);
						cmd = GMT_Create_Cmd (API, options);
						sprintf (line, "Command: %s %s", THIS_MODULE_CLASSIC_NAME, cmd);	/* Build command line argument string */
						gmt_M_free (GMT, cmd);
						GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, line);
						sprintf (line, "%s%s%s%s%c_1%s%c_2%sdist_1%sdist_2%shead_1%shead_2%svel_1%svel_2",
							s->info[s->out_order[s->x_col]].name, c, s->info[s->out_order[s->y_col]].name, c, t_or_i, c, t_or_i, c, c, c, c, c, c);
						for (j = 0; j < n_data_col; j++) {
							col = col_number[j];
							if (Ctrl->Z.active)
								sprintf (item, "%s%s_1%s%s_2", c, s->info[s->out_order[col]].name, c, s->info[s->out_order[col]].name);
							else
								sprintf (item, "%s%s_X%s%s_M", c, s->info[s->out_order[col]].name, c, s->info[s->out_order[col]].name);
							strcat (line, item);
						}
						GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, line);
						first_header = false;
					}

					if (first_crossover) {
						char info[GMT_BUFSIZ] = {""}, l_start[2][GMT_LEN64], stop[2][GMT_LEN64];
						for (k = 0; k < 2; k++) {
							if (has_time[k]) {	/* Find first and last record times */
								for (j = 0; j < n_rec[k] && gmt_M_is_dnan (time[k][j]); j++);	/* Find first non-NaN time */
								gmt_ascii_format_col (GMT, l_start[k], time[k][j], GMT_OUT, 2);
								for (j = n_rec[k]-1; j > 0 && gmt_M_is_dnan (time[k][j]); j--);	/* Find last non-NaN time */
								gmt_ascii_format_col (GMT, stop[k], time[k][j], GMT_OUT, 3);
							}
							else {
								strcpy (l_start[k], "NaN");
								strcpy (stop[k], "NaN");
							}
						}
						sprintf (info, "%s/%s/%g %s/%s/%g", l_start[SET_A], stop[SET_A], dist[SET_A][n_rec[SET_A]-1], l_start[SET_B], stop[SET_B], dist[SET_B][n_rec[SET_B]-1]);
						sprintf (line, x2sys_header, trk_name[A], data_set[SET_A].year, trk_name[B], data_set[SET_B].year, info);
						GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, line);
						first_crossover = false;
					}

					if (s->geographic) gmt_lon_range_adjust (s->geodetic, &out[GMT_X]);
					GMT_Put_Record (API, GMT_WRITE_DATA, Out);	/* Write this to output */
				}

				gmt_x_free (GMT, &XC);
			}

			if (!same) {	/* Must free up memory for B */
				x2sys_free_data (GMT, data[SET_B], s->n_out_columns, &data_set[SET_B]);
				gmt_M_free (GMT, dist[SET_B]);
				if (!got_time) gmt_M_free (GMT, time[SET_B]);
				gmt_M_free (GMT, ylist[SET_B]);
			}
			if (!Ctrl->C.active)
				GMT_Report (API, GMT_MSG_INFORMATION, "Processing %s - %s : %" PRIu64 "\n", trk_name[A], trk_name[B], nx);
			else {
				toc = clock();
				GMT_Report (API, GMT_MSG_INFORMATION, "Processing %s - %s : %" PRIu64 "\t%.3f sec\n", trk_name[A], trk_name[B], nx, (double)(toc - tic)/1000);
				if (fpC)	/* Save also the run time in file */
					fprintf (fpC, "%s\t%s\t%d\t%.3f\n", trk_name[A], trk_name[B], (int)nx, (double)(toc - tic)/1000);
			}
		}

		/* Must free up memory for A */

		x2sys_free_data (GMT, data[SET_A], s->n_out_columns, &data_set[SET_A]);
		gmt_M_free (GMT, dist[SET_A]);
		if (!got_time) gmt_M_free (GMT, time[SET_A]);
		gmt_M_free (GMT, ylist[SET_A]);
	}

	if (fpC) fclose (fpC);

	if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) {	/* Disables further data output */
		Return (API->error);
	}

	/* Free up other arrays */

	if (Ctrl->A.active) {	/* Free strings in pairs, then pairs itself */
		x2syscross_free_pairs (pair, n_pairs);
		gmt_M_free (GMT, pair);
	}

	gmt_M_free (GMT, Out);
	gmt_M_free (GMT, xdata[SET_A]);
	gmt_M_free (GMT, xdata[SET_B]);
	gmt_M_free (GMT, out);
	gmt_M_free (GMT, duplicate);
	if (n_data_col) {
		gmt_M_free (GMT, t);
		gmt_M_free (GMT, y);
		gmt_M_free (GMT, col_number);
		gmt_M_free (GMT, ok);
	}
	x2sys_free_list (GMT, trk_name, n_tracks);

	x2sys_end (GMT, s);

	Return (GMT_NOERROR);
}
