/*--------------------------------------------------------------------
 *	$Id: splitxyz.c,v 1.15 2004/04/25 09:10:46 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * 
 * read a file of lon, lat, zvalue[, distance, azimuth]
 * and split it into profile segments.
 * 
 * Author:	W. H. F. Smith
 * Date:	24 July, 1991-2000
 * Version:	4
 */

#include "gmt.h"

#define F_RES 1000	/* Number of points in filter halfwidth  */

struct  DATA    {
        double  d;
        double  a;
        double  x;
        double  y;
        double  z;
        double  w;
}       *data;

main(int argc, char **argv)
{
	int     i, ndata, n_alloc, nprofiles, begin, end, n_expected_fields, n_out, n_fields;
	int	n_required_fields;
	
	BOOLEAN	error = FALSE, dh_supplied = FALSE, ok, hilow, map_units = FALSE;
	BOOLEAN xy_only = FALSE;

        double  desired_azim = 90.0, tolerance_azim = 360.0;
        double  min_dist = 0.0, xy_filter = 0.0, z_filter = 0.0;
	double	d_gap = 10.0, course_change = 0.0;
	double	dy, dx, last_d, last_c, last_s, csum, ssum, this_c, this_s, dotprod;
	double	mean_azim, fwork[F_RES], *in, out[3];
	
 	char    buffer[BUFSIZ], basename[120], filename[120];
 	
	FILE    *fp = NULL, *fpout;

	void	filterz(struct DATA *data, int ndata, double z_filter, double *fwork, int hilow), filterxy_setup(double *fwork), filterxy(struct DATA *data, int begin, int end, double xy_filter, double *fwork, int hilow);

	argc = GMT_begin (argc, argv);
	basename[0] = '\0';	/* Later, if (basename[0]) we have to write files.  */
	fpout = GMT_stdout;

        for (i = 1; i < argc; i++) {
                if (argv[i][0] == '-') {
                        switch (argv[i][1]) {
 		
				/* Common parameters */
			
				case 'H':
				case 'V':
				case ':':
				case 'b':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;
				
				/* Supplemental parameters */
				
				case 'A':
                                        if ( (sscanf(&argv[i][2], "%lf/%lf", &desired_azim, &tolerance_azim)) != 2) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Can't decipher values\n", GMT_program);
						error = TRUE;
					}
                                        break;
				case 'C':
                                        if ( (sscanf(&argv[i][2], "%lf", &course_change)) != 1) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Can't decipher value\n", GMT_program);
						error = TRUE;
					}
                                        break;
                                case 'D':
                                        if ( (sscanf(&argv[i][2], "%lf", &min_dist)) != 1) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: Can't decipher value\n", GMT_program);
						error = TRUE;
					}
                                        break;
				case 'F':
                                        if ( (sscanf(&argv[i][2], "%lf/%lf", &xy_filter, &z_filter)) != 2) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: Can't decipher values\n", GMT_program);
						error = TRUE;
					}
                                        break;
                                case 'G':
                                        if ( (sscanf(&argv[i][2], "%lf", &d_gap)) != 1) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Can't decipher value\n", GMT_program);
						error = TRUE;
					}
                                        break;
				case 'M':
                                        map_units = TRUE;
					GMT_io.in_col_type[0] = GMT_io.out_col_type[0] = GMT_IS_LON;
					GMT_io.in_col_type[1] = GMT_io.out_col_type[1] = GMT_IS_LAT;
                                        break;
				case 'N':
                                        strcpy(basename, &argv[i][2]);
                                        break;
				case 'S':
                                        dh_supplied = TRUE;
                                        break;
				case 'Z':
                                        xy_only = TRUE;
                                        break;
                                default:
                                        error = TRUE;
					GMT_default_error (argv[i][1]);
                                        break;
                        }
                }
                else {
                        if ( (fp = GMT_fopen(argv[i], GMT_io.r_mode)) == NULL) {
                                fprintf(stderr,"%s:  Cannot open %s\n", GMT_program, argv[i]);
                                error = TRUE;
                        }
                }
        }

	if (argc == 1 || GMT_quick) {	/* Display usage */
 		fprintf (stderr, "splitxyz %s - Split xyz[dh] files into segments\n\n", GMT_VERSION);
		fprintf(stderr,"usage:  splitxyz [<xyz[dh]file>] -C<course_change> -D<minimum_distance>\n");
		fprintf(stderr,"\t[-A<azimuth>/<tolerance>] [-F<xy_filter>/<z_filter>] [-G<gap>] [-H[<nrec>]] [-M]\n");
		fprintf(stderr,"\t[-N<namestem>] [-S] [-V] [-Z] [-:] [-bi[s][<n>]] [-bo[s][<n>]] [-f[i|o]<colinfo>]\n\n");
 		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf(stderr,"\tGive xyz[dh]file name or read stdin.\n");
		fprintf(stderr,"\t-C  Profile ends when change of heading exceeds <course_change>.\n");
		fprintf(stderr,"\t-D Only write profile if length is at least <minimum_distance>.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf(stderr,"\t-A Only write profile if mean direction is w/in +/- <tolerance>\n");
		fprintf(stderr,"\t     of <azimuth>. [Default = All].\n");
		fprintf(stderr,"\t-F Filter the data.  Give full widths of cosine arch filters for xy and z.\n");
		fprintf(stderr,"\t   Defaults are both widths = 0, giving no filtering.\n");
		fprintf(stderr,"\t   Use negative width to highpass.\n");
		fprintf(stderr,"\t-G Do not let profiles have gaps exceeding <gap>. [Default = 10 dist units].\n");
		GMT_explain_option ('H');
		fprintf(stderr,"\t-M Map units TRUE; x,y in degrees, dist units in km.  [Default dist unit = x,y unit].\n");
		fprintf(stderr,"\t-N Write output to separate files named <namestem>.profile#.\n");
		fprintf(stderr,"\t     [Default all to stdout, separated by >].\n");
		fprintf(stderr,"\t-S dh is supplied.  Input is 5 col x,y,z,d,h with d non-decreasing.\n");
		fprintf(stderr,"\t     [Default input is 3 col x,y,z only and computes d,h from the data].\n");
		fprintf(stderr,"\t-Z No z-values.  Input is 2 col x,y only\n");
		GMT_explain_option ('V');
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf(stderr,"\t     Default input columns is set given -S and -Z options.\n");
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
        }

	if (min_dist <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: Minimum segment distance must be positive\n", GMT_program);
		error++;
	}
	if (course_change <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Course change tolerance must be positive\n", GMT_program);
		error++;
	}
	if (tolerance_azim < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Azimuth tolerance must be positive\n", GMT_program);
		error++;
	}
	if (d_gap < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Data gap distance must be positive\n", GMT_program);
		error++;
	}
	if (desired_azim < 0.0 || desired_azim > 360) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: azimuth must be in 0-360 degree range\n", GMT_program);
		error++;
	}
	if (xy_only && dh_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option: Cannot be used with -S option\n", GMT_program);
		error++;
	}
	if (xy_only && z_filter != 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: Cannot specify z-filter while using -Z option\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && gmtdefs.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_OUT] && !basename[0]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary output requires a namestem in -N\n", GMT_program);
		error++;
	}
	n_required_fields = (dh_supplied) ? 5 : ((xy_only) ? 2 : 3);
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN] = n_required_fields;
	if (GMT_io.binary[GMT_IN])
		n_expected_fields = GMT_io.ncol[GMT_IN];
	else
		n_expected_fields = n_required_fields;
	if (GMT_io.binary[GMT_IN] && n_expected_fields < n_required_fields) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data must have at least %d columns\n", GMT_program, n_required_fields);
		error++;
	}

	if (error) exit (EXIT_FAILURE);

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

#ifdef SET_IO_MODE
	GMT_setmode (GMT_OUT);
#endif

	tolerance_azim *= D2R;
	/* if (desired_azim > 180.0) desired_azim -= 180.0; */	/* Put in Easterly strikes  */
	desired_azim = D2R * (90.0 - desired_azim);	/* Work in cartesian angle and radians  */
	course_change *= D2R;

	if (fp == NULL) {
		fp = GMT_stdin;
#ifdef SET_IO_MODE
		GMT_setmode (GMT_IN);
#endif
	}
        n_alloc = GMT_CHUNK;
	ndata = 0;
	i = 0;

        data = (struct DATA *) GMT_memory (VNULL, (size_t)n_alloc, sizeof(struct DATA), GMT_program);

	if (!GMT_io.binary[GMT_IN] && gmtdefs.io_header[GMT_IN]) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (buffer, BUFSIZ, fp);

	while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

		if (GMT_io.status & GMT_IO_MISMATCH) {
			fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields,  n_expected_fields, i);
			exit (EXIT_FAILURE);
		
}
		i++;
                if (ndata == n_alloc) {
                        n_alloc += GMT_CHUNK;
                        data = (struct DATA *) GMT_memory ((void *)data, (size_t)n_alloc, sizeof(struct DATA), GMT_program);
                }
		if (!dh_supplied) {
			data[ndata].x = in[0];
			data[ndata].y = in[1];
			if (!xy_only) data[ndata].z = in[2];
			if (ndata) {
				dy = (data[ndata].y - data[ndata-1].y);
				dx = (data[ndata].x - data[ndata-1].x);
				if (map_units) {
					dy *= DEG2KM;
					dx *= (DEG2KM * cos (0.5 * (data[ndata].y + data[ndata-1].y) * D2R) );
				}
				if (dy == 0.0 && dx == 0.0) {
					data[ndata].d = data[ndata-1].d;
					data[ndata].a = data[ndata-1].a;
				}
				else {
					data[ndata].d = data[ndata-1].d + hypot(dx,dy);
					data[ndata].a = d_atan2(dy,dx);
				}
			}
			else {
				data[ndata].d = data[ndata].a = 0.0;
			}
			ndata++;
		}
		else {
			data[ndata].x = in[0];
			data[ndata].y = in[1];
			data[ndata].z = in[2];
			data[ndata].d = in[3];
			data[ndata].a = D2R * (90.0 - in[4]);
			ndata++;
		}
	}
	
	data = (struct DATA *) GMT_memory ((void *)data, (size_t)ndata, sizeof(struct DATA), GMT_program);
	if (!dh_supplied) data[0].a = data[1].a;
	if (fp != GMT_stdin) GMT_fclose(fp);

	/* Now we have read the data and can filter z if necessary.  */

	if (z_filter < 0.0) {
		hilow = TRUE;
		z_filter = fabs(z_filter);
		filterz(data, ndata, z_filter, fwork, hilow);
	}
	else if (z_filter > 0.0) {
		hilow = FALSE;
		filterz(data, ndata, z_filter, fwork, hilow);
	}

	if (xy_filter < 0.0) {
		hilow = TRUE;
		xy_filter = fabs(xy_filter);
		filterxy_setup(fwork);
	}
	else if (xy_filter > 0.0) {
		hilow = FALSE;
		filterxy_setup(fwork);
	}

	/* Now we are ready to search for segments.  */

	n_out = (xy_only) ? 2 : 3;
	nprofiles = 0;
	begin = 0;
	end = begin;
	while (end < ndata-1) {
		last_d = data[begin].d;
		last_c = cos(data[begin].a);
		last_s = sin(data[begin].a);
		csum = last_c;
		ssum = last_s;
		ok = TRUE;
		while (ok && end < ndata-1) {
			end++;
			if (data[end].d - last_d > d_gap) {
				/* Fails due to too much distance gap  */
				ok = FALSE;
				continue;
			}
			this_c = cos(data[end].a);
			this_s = sin(data[end].a);
			dotprod = this_c * last_c + this_s * last_s;
			if (fabs(dotprod) > 1.0) dotprod = copysign(1.0, dotprod);
			if (d_acos(dotprod) > course_change) {
				/* Fails due to too much change in azimuth  */
				ok = FALSE;
				continue;
			}
			/* Get here when this point belongs with last one:  */
			csum += this_c;
			ssum += this_s;
			last_c = this_c;
			last_s = this_s;
			last_d = data[end].d;
		}
		
		/* Get here when we have found a beginning and end  */

		if (ok) {
			/* Last point in input should be included in this segment  */
			end++;
		}

		if (end - begin - 1) {
		
			/* There are at least two points in the list.  */
			
			if ((data[end-1].d - data[begin].d) >= min_dist) {
			
				/* List is long enough.  Check strike. Compute mean_azim in range [-pi/2, pi/2]:  */
				
				/* mean_azim = d_asin(ssum/sqrt(csum*csum + ssum*ssum)); */
				mean_azim = d_atan2 (ssum, csum);
				mean_azim = fabs(mean_azim - desired_azim);
				/* if (mean_azim <= tolerance_azim || mean_azim >= (M_PI-tolerance_azim) ) { */
				if (mean_azim <= tolerance_azim) {
				
					/* List has acceptable strike.  */
					
					if (xy_filter != 0.0) filterxy(data, begin, end, xy_filter, fwork, hilow);
					nprofiles++;
					
					/* If (basename) we write separate files, else GMT_stdout with > marks:  */
					
					if (basename[0]) {
						sprintf(filename, "%s.profile%d", basename, nprofiles);
						if ( (fpout = GMT_fopen(filename, GMT_io.w_mode)) == NULL) {
							fprintf(stderr,"%s:  Cannot create %s\n", GMT_program, argv[i]);
							exit (EXIT_FAILURE);
						}
					}
					else
						fprintf (GMT_stdout, "> Start profile # %d\n", nprofiles);

					for (i = begin; i < end; i++) {
						out[0] = data[i].x;
						out[1] = data[i].y;
						if (!xy_only) out[2] = data[i].z;
						GMT_output (fpout, n_out, out);
					}

					if (fpout != GMT_stdout) {
						GMT_fclose (fpout);
						if (gmtdefs.verbose) fprintf (stderr,"%s: Wrote %d points to file %s\n", GMT_program, end-begin, filename);
					}
				}
			}
		}
		begin = end;
	}

	/* Get here when all profiles have been found and written.  */

	if (gmtdefs.verbose) fprintf(stderr,"%s:  Split %d data into %d files.\n", GMT_program, ndata, nprofiles);
	GMT_free ((void *)data);
	
	GMT_end (argc, argv);
}


void	filterz(struct DATA *data, int ndata, double z_filter, double *fwork, int hilow)
{
	
	double	half_width, sum, dt, tmp;
	int	i, j, k, istart, istop;

	half_width = 0.5 * z_filter;
	sum = 0.0;
	dt = F_RES/half_width;
	tmp = M_PI / F_RES;
	for (i = 0; i < F_RES; i++) {
		fwork[i] = 1.0 + cos(i*tmp);
		sum += fwork[i];
	}
	for (i = 1; i < F_RES; i++) {
		fwork[i] /= sum;
	}

	j = 0;
	istart = 0;
	istop = 0;
	while(j < ndata) {
		while(istart < ndata && data[istart].d - data[j].d <= -half_width) istart++;
		while(istop < ndata && data[istop].d - data[j].d < half_width) istop++;
		sum = 0.0;
		data[j].w = 0.0;
		for(i = istart; i < istop; i++) {
			k = (int)floor(dt*fabs(data[i].d - data[j].d));
			sum += fwork[k];
			data[j].w += (data[i].z * fwork[k]);
		}
		data[j].w /= sum;
		j++;
	}
	if (hilow) {
		for (i = 0; i < ndata; i++) data[i].z = data[i].z - data[i].w;
	}
	else {
		for (i = 0; i < ndata; i++) data[i].z = data[i].w;
	}
	return;
}

void	filterxy_setup(double *fwork)
{
	int	i;
	double	tmp, sum = 0.0;

	tmp = M_PI / F_RES;
	for (i = 0; i < F_RES; i++) {
		fwork[i] = 1.0 + cos(i*tmp);
		sum += fwork[i];
	}
	for (i = 1; i < F_RES; i++) {
		fwork[i] /= sum;
	}
	return;
}

void	filterxy(struct DATA *data, int begin, int end, double xy_filter, double *fwork, int hilow)
{
	int	i, j, k, istart, istop;
	double	half_width, dt, sum;

	half_width = 0.5 * fabs(xy_filter);
	dt = F_RES/half_width;

	j = begin;
	istart = begin;
	istop = begin;
	while(j < end) {
		while(istart < end && data[istart].d - data[j].d <= -half_width) istart++;
		while(istop < end && data[istop].d - data[j].d < half_width) istop++;
		sum = 0.0;
		data[j].w = 0.0;
		for(i = istart; i < istop; i++) {
			k = (int)floor(dt*fabs(data[i].d - data[j].d));
			sum += fwork[k];
			data[j].w += (data[i].x * fwork[k]);
		}
		data[j].w /= sum;
		j++;
	}
	if (hilow) {
		for (i = begin; i < end; i++) data[i].x = data[i].x - data[i].w;
	}
	else {
		for (i = begin; i < end; i++) data[i].x = data[i].w;
	}

	j = begin;
	istart = begin;
	istop = begin;
	while(j < end) {
		while(istart < end && data[istart].d - data[j].d <= -half_width) istart++;
		while(istop < end && data[istop].d - data[j].d < half_width) istop++;
		sum = 0.0;
		data[j].w = 0.0;
		for(i = istart; i < istop; i++) {
			k = (int)floor(dt*fabs(data[i].d - data[j].d));
			sum += fwork[k];
			data[j].w += (data[i].y * fwork[k]);
		}
		data[j].w /= sum;
		j++;
	}
	if (hilow) {
		for (i = begin; i < end; i++) data[i].y = data[i].y - data[i].w;
	}
	else {
		for (i = begin; i < end; i++) data[i].y = data[i].w;
	}

	return;
}


