/*--------------------------------------------------------------------
 *	$Id: originator.c,v 1.3.4.2 2002/02/27 17:41:10 pwessel Exp $
 *
 *   Copyright (c) 2000-2002 by P. Wessel
 *
 *   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/wessel
 *--------------------------------------------------------------------*/
/*
 * originator reads file of seamount locations and tries to match each
 * seamount with a probable hotspot by drawing flowines back in time and
 * keeping track of which hotspot is closest to each flowline.  It then
 * reports the closest hotspot, the stage of the flowline involved, the
 * implied pseudo-age of the seamount, and the minimum distance between
 * the flowline and hotspot (in km).
 *
 * Author:	Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA
 * Date:	29-DEC-1999
 * Version:	1.0
 *
 *-------------------------------------------------------------------------
 * The ASCII Euler file must have following format:
 *
 * 1. Any number of comment lines starting with # in first column
 * 2. Any number of blank lines (just carriage return, no spaces)
 * 2. Any number of stage pole records which each have the format:
 *    lon(deg)  lat(deg)  tstart(Ma)  tstop(Ma)  ccw-angle(deg)
 * 3. stage records must go from oldest to youngest rotation
 * 4. Note tstart is larger (older) that tstop for each record
 * 5. No gaps allowed: tstart must equal the previous records tstop
 *
 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
 *
 * # Time in Ma, angles in degrees
 * # lon  lat	tstart	tend	ccw-angle
 * 165     85	150	100	24.0
 * 284     36	100	74	15.0
 * 265     22	74	65	7.5
 * 253     17	65	42	14.0
 * 285     68	42	0	34.0
 *
 * ASCII seamount location file(s) must have the following format:
 *
 * 1. Any number of comment lines starting with # in first column
 * 2. Any number of blank lines (just carriage return, no spaces)
 * 3. For special header records, see -H
 * 4. Any number of data recordswhich each have the format:
 *    lon lat height radius crustal_age    (or lat lon ..., see -: option).
 *    crustal_age in Ma, height and radius are not used by originator but
 *    are used by hotspotter.
 *
 * Binary files cannot have header records, and data fields must all be
 * either single or double precision (see -bi option).  Output file will
 * be ASCII since it contains a text string (hotspot ID).
 *
 * The file with a list of hotspots must have the following format:
 *
 * 1. Any number of comment lines starting with # in first column
 * 2. Any number of blank lines (just carriage return, no spaces)
 * 2. Any number of hotspot records which each have the format:
 *    lon(deg)  lat(deg)  id  name
 *    the id is a 3-character tag (e.g., HWI), the name is the
 *    full name of the hotspot (e.g., Hawaii).
 *
 * Example: 
 *
 * # Partial list (Pacific) of HotSpots from Table 1 of Yamaji, 1992
 * #Lon		Lat	Abbreviation	Hotspot_name
 * 167		3	CRL	Caroline
 * 230		46	COB	Cobb
 * 205		20	HWI	Hawaii
 * 221.9	-50.9	LSV	Louisville
 * 220		-29	MDN	MacDonald
 * 221		-11	MRQ	Marquesas
 * 231		-27	PTC	Pitcairn
 * 254		-27	SLG	Sala y Gomez
 * 192		-15	SAM	Samoa
 * 212		-18	SOC	Society
 */
 
#include "spotter.h"

struct HOTSPOT {
	double lon, lat;	/* Current location of hot spot */
	double dist;		/* Distance to current seamount */
	double t;		/* Distance to current seamount */
	int id;			/* Hot spot id flag */
	int stage;		/* Stage to which seamount belongs */
	char name[20];		/* Name of hotspot */
	char abbrev[4];		/* 3-char abbreviation of hotspot name */
};

int hotspot_init (char *file, struct HOTSPOT **p);
int comp_hs (const void *p1, const void *p2);

main (int argc, char **argv)
{
	int i, j, k, n, kk, ns, nh, nc, np, x, y, n_out, fno, n_files = 0, n_args;
	int n_fields, n_expected_fields, n_read, n_skipped = 0, n_best_hs = 1;
		
	BOOLEAN	error = FALSE, done, nofile = TRUE, truncate_ages = FALSE, first = TRUE;
	BOOLEAN finite = FALSE;		/* TRUE if stage pole file contains finite rotation poles instead */
	
	double	x_smt, y_smt, z_smt, r_smt, t_smt, *c, *in, d_km = 5.0, upper_age = 180.0, dist;
	double max_dist = 1.0e100;
	
	char *hsfile = CNULL, *efile = CNULL, age[32];
	
	FILE *fp;
	
	struct EULER *p;
	struct HOTSPOT *hotspot, *hot;

	argc = GMT_begin (argc, argv);
		
	/* Check command line arguments */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
			
				/* Common parameters */
			
				case 'H':
				case 'V':
				case ':':
				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 'C':	/* Use finite rotation poles */
					finite = TRUE;
					break;

				case 'D':
					d_km = atof (&argv[i][2]);
					break;
				case 'E':
					efile = &argv[i][2];
					break;
				case 'F':
					hsfile = &argv[i][2];
					break;
				case 'N':
					upper_age = atof (&argv[i][2]);
					break;
				case 'S':
					n_best_hs = atoi (&argv[i][2]);
					break;
				case 'T':
					truncate_ages = TRUE;
					break;
				case 'W':
					max_dist = atof (&argv[i][2]);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}

	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] = 5;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < 5) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 5 columns\n", GMT_program);
		error++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "%s %s - Associate seamounts with hotspot point sources\n\n", GMT_program, GMT_VERSION);
		fprintf (stderr, "usage: %s [<xyfiles>] -E<euler_file> -F<hotspot_file> [-C] [-D<d_km>]\n", GMT_program);
		fprintf (stderr, "	[-H] [-N<upper_age>] [-S<n_hs>] [-T] [-V] [-W<maxdist>] [-:]\n\n");
		
		if (GMT_quick) exit (-1);
		
		fprintf (stderr, "\txyfiles is one or more seamount (x,y,z,r,t) files\n");
		fprintf (stderr, "\t-E Specify file name for Euler poles.\n");
		fprintf (stderr, "\t-F Specify file name for hotspot locations.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C The file given with -E contains finite rotation poles [Default is stage poles]\n");
		fprintf (stderr, "\t-D set sampling interval in km along tracks [5].\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-N set age (in m.y.) for seafloor where age == NaN [180].\n");
		fprintf (stderr, "\t-S Report the <n_hs> closest hotSpots [1].\n");
		fprintf (stderr, "\t-T Truncate seamount ages exceeding the upper age set with -N [no truncation] \n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W Only report seamounts whose closest encounter to a hotspot is less than <maxdist> km\n");
		fprintf (stderr, "\t   [Default reports for all seamounts] \n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "	   Default is 5 input columns\n");
		exit (EXIT_FAILURE);
	}

	if (!hsfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -F:  Must specify hotspot file\n", GMT_program);
		error = TRUE;
	}
	if (!efile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E:  Must specify Euler pole file\n", GMT_program);
		error = TRUE;
	}
	if (d_km <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -D:  Must specify a positive interval\n", GMT_program);
		error = TRUE;
	}
	if (max_dist <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -W:  Must specify a positive distance in km\n", GMT_program);
		error = TRUE;
	}
	
	if (error) exit (EXIT_FAILURE);

	nh = hotspot_init (hsfile, &hotspot);
	if (n_best_hs <= 0 || n_best_hs > nh) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S:  Give value between 1 and %d\n", GMT_program, nh);
		exit (EXIT_FAILURE);
	}
	n_out = MIN (n_best_hs, nh);

	ns = spotter_init (efile, &p, TRUE, finite, &upper_age);

	hot = (struct HOTSPOT *) GMT_memory (VNULL, nh, sizeof (struct HOTSPOT), GMT_program);

	x = (gmtdefs.xy_toggle) ? 1 : 0;	y = 1 - x;		/* Set up which columns have x and y */

	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
	
	done = FALSE;
	
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 5;
	n = 0;
	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over all input files */

		if (!nofile && argv[fno][0] == '-') continue;
		if (nofile) {
			fp = stdin;
			done = TRUE;
		}
		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]);
		
		n_read = 0;
		if (gmtdefs.io_header) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (GMT_io.segment_header, BUFSIZ, fp);
				if (first && !GMT_io.binary[1]) fprintf (GMT_stdout, "%s", GMT_io.segment_header);
				n_read++;
			}
			first = FALSE;
		}
		
		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */
			n_read++;
			while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) {
				GMT_write_segmentheader (GMT_stdout, n_expected_fields);
				n_fields = GMT_input (fp, &n_expected_fields, &in);
				n_read++;
			}
			if (GMT_io.status & GMT_IO_EOF) continue;

			if (GMT_io.status & GMT_IO_MISMATCH) {
				fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d (skipped)\n", GMT_program, n_fields, n_expected_fields, n_read);
				continue;
			}
			if (GMT_is_dnan (in[4])) {	/* Age is NaN, assign value */
				t_smt = upper_age;
			}
			else {			/* Assign given value, truncate if necessary */
				t_smt = in[4];
				if (t_smt > upper_age) {
					if (truncate_ages) {
						t_smt = upper_age;
					}
					else {
						fprintf (stderr, "%s: Seamounts near line %d has age (%lg) > oldest stage (%lg) (skipped)\n", GMT_program, n_read, t_smt, upper_age);
						continue;
					}
				}
			}
			if (t_smt < 0.0) {	/* Negative ages are flags for points to be skipped */
				n_skipped++;
				continue;
			}
			x_smt = in[0] * D2R;
			y_smt = in[1] * D2R;
			z_smt = in[2];
			r_smt = in[3];

			if (gmtdefs.verbose && !(n % 10)) fprintf (stderr, "%s: Working on seamount # %5d\r", GMT_program, n);

			nc = spotter_forthtrack (&x_smt, &y_smt, &t_smt, 1, p, ns, d_km, 0.0, TRUE, &c);

			np = (int) c[0];

			memcpy ((void *)hot, (void *)hotspot, nh * sizeof (struct HOTSPOT));

			for (kk = 0, k = 1; kk < np; kk++, k += 3) {	/* For this seamounts track */
				for (j = 0; j < nh; j++) {	/* For all hotspots */
					dist = GMT_great_circle_dist (hot[j].lon, hot[j].lat, R2D * c[k], R2D * c[k+1]) * KM_PR_DEG;
					if (dist < hot[j].dist) {
						hot[j].dist = dist;
						hot[j].t = c[k+2];
					}
				}
			}
			for (j = 0; j < nh; j++) {	/* Assign stage id */
				k = 0;
				while (k < ns && hot[j].t <= p[k].t_stop) k++;
				hot[j].stage = ns - k;
				if (hot[j].stage == 0) hot[j].stage++;
			}
			
			qsort ((void *)hot, nh, sizeof(struct HOTSPOT), comp_hs);
			
			if (hot[0].dist < max_dist) {
				if (t_smt == 180.0)
					strcpy (age, "NaN");
				else
					sprintf (age, "%lg", t_smt);
				
				fprintf (GMT_stdout, "%lg\t%lg\t%lg\t%lg\t%s", in[x], in[y], z_smt, r_smt, age);
				for (j = 0; j < n_out; j++) fprintf (GMT_stdout, "\t%s\t%d\t%lg\t%lg", hot[j].abbrev, hot[j].stage, hot[j].t, hot[j].dist);
				fprintf (GMT_stdout, "\n");
			}

			GMT_free ((void *)c);
			n++;
		}
		
		if (fp != stdin) fclose (fp);
	}
		
	if (gmtdefs.verbose) fprintf (stderr, "%s: Working on seamount # %5d\n", GMT_program, n);

	GMT_free ((void *)hotspot);
	GMT_free ((void *)hot);
	GMT_free ((void *)p);
	
	GMT_end (argc, argv);
}

int comp_hs (const void *p1, const void *p2)
{
	struct HOTSPOT *a, *b;

	a = (struct HOTSPOT *) p1;
	b = (struct HOTSPOT *) p2;
	if (a->dist < b->dist) return (-1);
	if (a->dist > b->dist) return (1);
	return (0);
}

int hotspot_init (char *file, struct HOTSPOT **p)
{
	FILE *fp;
	struct HOTSPOT *e;
	char buffer[BUFSIZ];
	int i = 0;
	size_t n_alloc = GMT_CHUNK;

	if ((fp = fopen (file, "r")) == NULL) {
		fprintf (stderr, "%s: Cannot open file %s - aborts\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	e = (struct HOTSPOT *) GMT_memory (VNULL, n_alloc, sizeof (struct HOTSPOT), GMT_program);

	while (fgets (buffer, 512, fp) != NULL) {
		if (buffer[0] == '#' || buffer[0] == '\n') continue;
		sscanf (buffer, "%lf %lf %s %s", &e[i].lon, &e[i].lat, e[i].abbrev, e[i].name);
		e[i].id = i;
		e[i].dist = 1.0e100;
		i++;
		if ((size_t)i == n_alloc) {
			n_alloc += GMT_CHUNK;
			e = (struct HOTSPOT *) GMT_memory ((void *)e, n_alloc, sizeof (struct HOTSPOT), GMT_program);
		}
	}
	fclose (fp);
	e = (struct HOTSPOT *) GMT_memory ((void *)e, (size_t)i, sizeof (struct HOTSPOT), GMT_program);
	*p = e;

	return (i);
}