/*--------------------------------------------------------------------
 *	$Id: grdvector.c,v 1.4.4.2 2002/01/17 17:42:51 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 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
 *--------------------------------------------------------------------*/
/*
   grdvector reads 2 grdfiles that contains the 2 components of a vector
   field (cartesian or polar) and plots vectors at the grid positions.
   This is basically a short-hand for using grd2xyz | psxy -SV and is
   more convenient for such plots on a grid.

   Author:	Paul Wessel
   Date:	12-JUN-1995
   Revised:	15-FEB-2000
   Version:	3.4.1
 
 */


#include "gmt.h"

float *r, *theta;

main (int argc, char **argv)
{

	int	i, j, n = 0, nm, nx, ny, ij, i0, j0, di, dj, off;
	
	BOOLEAN convert_angles = FALSE, get_rgb = FALSE, cartesian = TRUE, shrink = FALSE, set_fill = FALSE;
	BOOLEAN error = FALSE, center = FALSE, outline = FALSE, azimuth = FALSE, stick_plot = TRUE, inc_set = FALSE;
	BOOLEAN clip = TRUE, got_fix_length = FALSE;
	
	char *file[2], *cpt, txt_a[32], txt_b[32], txt_c[32], unit = 0;
	
	double dx2, dy2, v_width = 0.03, h_length = 0.12, h_width = 0.1;
	double v_w, h_l, h_w, v_shrink, v_norm = 0.0, tmp, x, y, plot_x, plot_y, x_off, y_off;
	double west, east, south, north, x2, y2, scale = 1.0, fix_length = 0.0;
	double data_west, data_east, data_south, data_north, value, c, s;
	
	struct GRD_HEADER h[2];
	struct GMT_FILL fill;
	struct GMT_PEN pen;

	GMT_init_pen (&pen, GMT_PENWIDTH);
	GMT_init_fill (&fill, -1, -1, -1);
	west = east = south = north = 0.0;
	di = dj = 1;
	i0 = j0 = 0;
	dx2 = dy2 = 0.0;
	
	argc = GMT_begin (argc, argv);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* Common parameters */
			
				case 'B':
				case 'J':
				case 'K':
				case 'O':
				case 'P':
				case 'R':
				case 'U':
				case 'V':
				case 'X':
				case 'x':
				case 'Y':
				case 'y':
				case 'c':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
			
				case 'A':
					cartesian = FALSE;
					break;
				case 'C':	/* Vary symbol color with z */
					cpt = &argv[i][2];
					get_rgb = TRUE;
					break;
				case 'E':
					center = TRUE;
					break;
				case 'G':		/* Set Gray shade for polygon */
					if (GMT_getfill (&argv[i][2], &fill)) {
						GMT_fill_syntax ('G');
						error++;
					}
					set_fill = TRUE;
					break;
				case 'I':	/* Only use gridnodes dx2,dy2 apart */
					GMT_getinc (&argv[i][2], &dx2, &dy2);
					inc_set = TRUE;
					break;
				case 'N':	/* Do not clip at border */
					clip = FALSE;
					break;
				case 'Q':
					if (argv[i][2] && argv[i][3] != 'n') {
						if (sscanf (&argv[i][2], "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c) != 3) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -Q option:  Could not decode arrowwidth/headlength/headwidth\n", GMT_program);
							error++;
						}
						else {
							v_width  = GMT_convert_units (txt_a, GMT_INCH);
							h_length = GMT_convert_units (txt_b, GMT_INCH);
							h_width  = GMT_convert_units (txt_c, GMT_INCH);
						}
					}
					for (j = 2; argv[i][j] && argv[i][j] != 'n'; j++);
					if (argv[i][j]) {	/* Normalize option used */
						v_norm = atof (&argv[i][j+1]);
						if (v_norm > 0.0) {
							v_shrink = 1.0 / v_norm;
							shrink = TRUE;
						}
						else {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -Qn option:  No reference length given\n", GMT_program);
							error++;
						}
					}
					stick_plot = FALSE;
					break;
				case 'S':
					j = strlen (argv[i]) - 1;
					if (strchr ("cimpCIMP", (int)argv[i][j])) unit = argv[i][j];
					if (argv[i][2] == 'l' || argv[i][2] == 'L') {
						got_fix_length = TRUE;
						fix_length = atof (&argv[i][3]);
					}
					else
						scale = atof (&argv[i][2]);
					break;
				case 'T':
					convert_angles = TRUE;
					break;
				case 'W':		/* Set line attributes */
					if (argv[i][2] && GMT_getpen (&argv[i][2], &pen)) {
						GMT_pen_syntax ('W');
						error++;
					}
					outline = TRUE;
					break;
				case 'Z':
					azimuth = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else if (n < 2)
			file[n++] = argv[i];
		else
			n++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdvector %s - Plot vector fields from grdfiles\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdvector compx.grd compy.grd -J<params> -R<west/east/south/north> [-A]\n");
		fprintf (stderr, "\t[-B<tickinfo>] [-C<cpt>] [-E] [-G<fill>] [-I<dx/dy>] [-K] [-O] [-P] [-Q<params>] [-N] [-S[l]<scale>] [-T]\n");
		fprintf (stderr, "\t[-U[<label>]] [-V] [-W<pen>] [-X<x_shift>] [-Y<y_shift>] [-Z] [-c<ncopies>] [-:]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tcompx & compy are grdfiles with the 2 vector components.\n");
		GMT_explain_option ('j');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A means grdfiles have polar (r, theta) components [Default is Cartesian]\n");
		GMT_explain_option ('b');
		fprintf (stderr, "\t-C Use cpt-file to assign colors based on vector length\n");
		fprintf (stderr, "\t-E cEnter vectors on grid nodes [Default draws from grid node]\n");
		fprintf (stderr, "\t-G Specify color. fill can be either <r/g/b> (each 0-255) for color or (0-255) for gray [0].\n");
		fprintf (stderr, "\t   Default is no fill (vector outlines only)\n");
		fprintf (stderr, "\t-I plots only those nodes that are <dx/dy> apart [Default is all nodes]\n");
		GMT_explain_option ('K');
		fprintf (stderr, "\t-N Do Not clip vectors that exceed the map boundaries [Default will clip]\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr, "\t-Q Select vector plot [Default is stick-plot].\n");
		fprintf (stderr, "\t   Optionally, specify vector parameters\n");
		fprintf (stderr, "\t   <params> are arrowwidth/headlength/headwidth [Default is 0.03i/0.12i/0.09i]\n");
		fprintf (stderr, "\t   Append n<size> which will cause vectors shorter than <size> to be\n");
		fprintf (stderr, "\t     scaled down\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\t-S sets scale for vector length in data units per %s [1]\n", GMT_unit_names[gmtdefs.measure_unit]);
		fprintf (stderr, "\t   Append c, i, m, or p to indicate cm, inch, m, or points as the distance unit.\n");
		fprintf (stderr, "\t   Alternatively, prepend l to indicate a fixed length for all vectors.\n");
		fprintf (stderr, "\t-T means azimuth should be converted to angles based on map projection\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets pen attributes [width = %lgp, color = (%d/%d/%d), texture = solid line].\n", 
			pen.width, pen.rgb[0], pen.rgb[1], pen.rgb[2]);
		GMT_explain_option ('X');
		fprintf (stderr, "\t-Z means the angles provided are azimuths rather than direction\n");
		GMT_explain_option ('c');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	if (inc_set && (dx2 <= 0.0 || dy2 <= 0.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify positive increments\n", GMT_program);
		error++;
	}
	if (scale == 0.0 && !got_fix_length) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Scale must be nonzero\n", GMT_program);
		error++;
	}
	if (fix_length <= 0.0 && got_fix_length) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Sl option:  Length must be positive\n", GMT_program);
		error++;
	}
	if (got_fix_length && shrink) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Sl, -Q options:  Cannot use -Q..n<size> with -Sl\n", GMT_program);
		error++;
	}
	if (azimuth && cartesian) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  Azimuths not valid input for Cartesian data\n", GMT_program);
		error++;
	}
	if (get_rgb && !cpt) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Must specify a color palette table\n", GMT_program);
		error++;
	}
	if (!(set_fill || outline || get_rgb)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one of -G, -W, -C\n", GMT_program);
		error++;
	}
	if (n != 2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify two input grdfiles\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */

	if (get_rgb) GMT_read_cpt (cpt);

	if (!(strcmp (file[0], "=") || strcmp (file[1], "="))) {
		fprintf (stderr, "%s: Piping of grdfiles not supported!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (GMT_read_grd_info (file[0], &h[0])) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file[0]);
		exit (EXIT_FAILURE);
	}
	
	if (GMT_read_grd_info (file[1], &h[1])) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file[1]);
		exit (EXIT_FAILURE);
	}
	
	if (!(h[0].nx == h[1].nx && h[0].ny == h[1].ny && h[0].x_min == h[1].x_min && h[0].y_min == h[1].y_min 
		&& h[0].x_inc == h[1].x_inc && h[0].y_inc == h[1].y_inc)) {
		fprintf (stderr, "%s: files %s and %s does not match!\n", GMT_program, file[0], file[1]);
		exit (EXIT_FAILURE);
	}
	off = (h[0].node_offset) ? 0 : 1;
		
	/* Determine what wesn to pass to map_setup */

	if (!project_info.region_supplied) {
		west = h[0].x_min;
		east = h[0].x_max;
		south = h[0].y_min;
		north = h[0].y_max;
	}

	GMT_map_setup (west, east, south, north);

	/* Determine the wesn to be used to read the grdfile */

	GMT_grd_setregion (&h[0], &data_west, &data_east, &data_south, &data_north);

	/* Read data */

	nx = irint ( (data_east - data_west) / h[0].x_inc) + off;
	ny = irint ( (data_north - data_south) / h[0].y_inc) + off;
	nm = nx * ny;
	r = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
	theta = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
	if (GMT_read_grd (file[0], &h[0], r, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file[0]);
		exit (EXIT_FAILURE);
	}
	if (GMT_read_grd (file[1], &h[1], theta, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file[1]);
		exit (EXIT_FAILURE);
	}
	if (!got_fix_length) scale = 1.0 / scale;

	switch (unit) {	/* Adjust for possible unit selection */
		case 'C':
		case 'c':
			scale *= GMT_u2u[GMT_CM][GMT_INCH];
			fix_length *= GMT_u2u[GMT_CM][GMT_INCH];
			break;
		case 'I':
		case 'i':
			scale *= GMT_u2u[GMT_INCH][GMT_INCH];
			fix_length *= GMT_u2u[GMT_INCH][GMT_INCH];
			break;
		case 'M':
		case 'm':
			scale *= GMT_u2u[GMT_M][GMT_INCH];
			fix_length *= GMT_u2u[GMT_M][GMT_INCH];
			break;
		case 'P':
		case 'p':
			scale *= GMT_u2u[GMT_PT][GMT_INCH];
			fix_length *= GMT_u2u[GMT_PT][GMT_INCH];
			break;
		default:
			scale *= GMT_u2u[gmtdefs.measure_unit][GMT_INCH];
			fix_length *= GMT_u2u[gmtdefs.measure_unit][GMT_INCH];
			break;
	}
	
	ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
		gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
		gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0]));
	GMT_echo_command (argc, argv);
	if (gmtdefs.unix_time) GMT_timestamp (argc, argv);

	GMT_setpen (&pen);

        if (clip) GMT_map_clip_on (GMT_no_rgb, 3);

	if (dx2 != 0.0 && dy2 != 0.0) {
		dj = irint (dy2 / h[0].y_inc);
		di = irint (dx2 / h[0].x_inc);
		tmp = ceil (h[0].y_max / dy2) * dy2;
		if (tmp > h[0].y_max) tmp -= dy2;
		j0 = irint ((h[0].y_max - tmp) / h[0].y_inc);
		tmp = floor (h[0].x_min / dx2) * dx2;
		if (tmp < h[0].x_min) tmp += dx2;
		i0 = irint ((tmp - h[0].x_min) / h[0].x_inc);
	}
	
	for (j = j0; j < h[1].ny; j += dj) {
		y = h[0].y_max - j *h[0].y_inc;
		for (i = i0; i < h[1].nx; i += di) {
		
			ij = j * h[0].nx + i;
			if (GMT_is_fnan (r[ij]) || GMT_is_fnan (theta[ij])) continue;	/* Cannot plot NaN-vectors */
			
			value = r[ij];
			
			if (cartesian) {
				value = hypot (theta[ij], r[ij]);
				if (value == 0.0) continue;
				theta[ij] = (float)(R2D * atan2 (theta[ij], r[ij]));
				r[ij] = (float)value;
			}
			else if (r[ij] < 0.0) {
				r[ij] = -r[ij];
				theta[ij] += 180.0;
			}
			else if (r[ij] == 0.0) continue;
			
			if (get_rgb) GMT_get_rgb24 (value, fill.rgb);
			
			x = h[0].x_min + i * h[0].x_inc;
			GMT_geo_to_xy (x, y, &plot_x, &plot_y);
			
			if (convert_angles) {
				if (!azimuth) theta[ij] = (float)90.0 - theta[ij];
				GMT_azim_to_angle (x, y, 0.1, (double)theta[ij], &tmp);
				theta[ij] = (float)(tmp * D2R);
			}
			else
				theta[ij] *= (float)D2R;
			
			if (got_fix_length)
				r[ij] = (float)fix_length;
			else
				r[ij] *= (float)scale;
			
			sincos (theta[ij], &s, &c);
			x2 = plot_x + r[ij] * c;
			y2 = plot_y + r[ij] * s;
			
			if (center) {
				x_off = 0.5 * (x2 - plot_x);
				y_off = 0.5 * (y2 - plot_y);
				plot_x -= x_off;
				plot_y -= y_off;
				x2 -= x_off;
				y2 -= y_off;
			}
			
			if (stick_plot) {
				if (get_rgb) ps_setpaint (fill.rgb);
				ps_plot (plot_x, plot_y, 3);
				ps_plot (x2, y2, 2);
				continue;
			}
			
			if (shrink && r[ij] < v_norm) {	/* Scale arrow attributes down with length */
				v_w = v_width * r[ij] * v_shrink;
				h_l = h_length * r[ij] * v_shrink;
				h_w = h_width * r[ij] * v_shrink;
				ps_vector (plot_x, plot_y, x2, y2, v_w, h_l, h_w, gmtdefs.vector_shape, fill.rgb, outline);
			}
			else	/* Leave as specified */
				ps_vector (plot_x, plot_y, x2, y2, v_width, h_length, h_width, gmtdefs.vector_shape, fill.rgb, outline);
		}
	}
	
        if (clip) GMT_map_clip_off ();

	if (frame_info.plot) GMT_map_basemap ();

	ps_plotend (gmtdefs.last_page);

	GMT_free ((void *)r);
	GMT_free ((void *)theta);
	
	GMT_end (argc, argv);
}