/*--------------------------------------------------------------------
 *	$Id: grdmath_main.c,v 1.33 2004/08/19 04:13:48 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
 *--------------------------------------------------------------------*/
/*
 * grdmath.c is a reverse polish calculator that operates on grd files 
 * (and constants) and perform basic mathematical operations
 * on them like add, multiply, etc.
 * Some operators only work on one operand (e.g., log, exp)
 *
 * Author:	Paul Wessel
 * Date:	22-AUG-1988
 * Revised:	27-JUL-2000 PW: Added EXTREMA operator
 *		18-AUG-2000 PW: Added LT, LE, EQ, GE, GT, NAN, ISNAN, XOR, MODE, MAD, LMSSCL
 *		20-AUG-2000 PW: Added FDIST, TDIST, CHI2DIST
 *		23-AUG-2000 PW: Added LOWER and UPPER
 *		28-NOV-2001 PW: Added Critical values for Chi2, F, T, and Z distributions
 *		12-NOV-2002 PW: Added Cartesian and Spherical Azimuth operators
 *		27-JAN-2004 PW: Added SINC and NEQ
 *		28-MAR-2004 PW: Added FLIPUD, FLIPLR, ROTX, ROTY, INRANGE, LDIST, PDIST, INSIDE
 *		01-JUL-2004 PW: Added LRAND
 *		17-JUL-2004 PW: Added LOG2
 * Version:	4
 *
 */
 
#include "gmt.h"

#define N_OPERATORS	117
#define STACK_SIZE	50

#define ARG_IS_FILE		0
#define ARG_IS_NUMBER		1
#define ARG_IS_PI		2
#define ARG_IS_E		3
#define ARG_IS_X_MATRIX		4
#define ARG_IS_Y_MATRIX		5
#define ARG_IS_ASCIIFILE	6
#define ARG_IS_SAVE		10
#define ARG_IS_OPERATOR		100
#define N_SPECIAL_ARGS		10
#define N_SKIP_ARGS		5

#include "grdmath_def.h"

PFV call_operator[N_OPERATORS];

int nm = 0, consumed_operands[N_OPERATORS], produced_operands[N_OPERATORS];

char ASCII_file[BUFSIZ];

struct GMT_HASH localhashnode[HASH_SIZE];

float *grd_x, *grd_y;

struct GRD_HEADER header;

int decode_argument(char *txt), get_operator(char *choice);
void grdmath_init(PFV ops[], int n_args[], int n_out[]);

char *special_arg[N_SPECIAL_ARGS] = {
	"PI",
	"pi",
	"E",
	"e",
	"x",
	"X",
	"y",
	"+",
	"-",
	"^"
};

char *skip_args[N_SKIP_ARGS] = {	/* 2nd time we want to skip these */
	"-R",
	"-I",
	"-F",
	"-V",
	"-f"
};

int GMT_mode_selection = 0, GMT_n_multiples = 0;

main (int argc, char **argv)
{
	int i, j, k, arg, op = 0, nstack = 0, new_stack = -1, ok = 1, type, n_items = 0, this_stack;
	
	BOOLEAN constant[STACK_SIZE], error = FALSE, skip_arg, set_r = FALSE, set_inc = FALSE;
	
	float *stack[STACK_SIZE];

	double factor[STACK_SIZE], d, x_noise, y_noise;
	
	char *outfile = CNULL, file[BUFSIZ];
	
	struct GRD_HEADER grd[STACK_SIZE];
	
	argc = GMT_begin (argc, argv);
	
	if (argc == 2 && !strcmp (argv[1], "-")) error = GMT_quick = TRUE;
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdmath %s - Reverse Polish Notation (RPN) calculator for grdfiles (element by element)\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdmath [-R<west/east/south/north> -F -I<xinc[m|c]>[/<yinc>[m|c] -V] A B op C op D op ... = outfile\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tA, B, etc are grdfiles, constants, or symbols (see below)\n");
		fprintf (stderr, "\tThe stack can hold up to %d entries (given enough memory)\n", STACK_SIZE);
		fprintf (stderr, "\tTrigonometric operators expect radians.  The operators are:\n\n");
		fprintf (stderr, "\tName	#args	Returns:\n");
		fprintf (stderr, "\t-----------------------\n");
#include "grdmath_explain.h"
		fprintf (stderr, "\n\tThe special symbols are:\n\n");
		fprintf (stderr, "\t  PI	= 3.1415926...\n");
		fprintf (stderr, "\t  E	= 2.7182818...\n");
		fprintf (stderr, "\t  X	= grid with x-coordinates\n");
		fprintf (stderr, "\t  Y	= grid with y-coordinates\n");
		fprintf (stderr, "\n\tOPTIONS: (only used if no grdfiles are passed as arguments)\n\n");
		fprintf (stderr, "\t-F Set pixel grid registration [Default is gridline orientation]\n");
		fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
		fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or c] to indicate minutes [or seconds].\n");
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}
	
	if (argv[argc-2][0] != '=') {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Usage is <operations> = outfile\n", GMT_program);
		exit (EXIT_FAILURE);
	}

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

	/* Determine if the user have files with stupid names that will cause conflict */

	for (op = 0; op < N_OPERATORS; op++) {
		if (!access (operator[op], R_OK))
			fprintf (stderr, "%s Warning: Your file %s may be confused with a %s operator!\n", GMT_program, operator[op], GMT_program);
	}
	for (i = 0; i < N_SPECIAL_ARGS; i++) {
		if (!access (special_arg[i], R_OK))
			fprintf (stderr, "%s Warning: Your file %s may be confused with a %s operator!\n", GMT_program, special_arg[i], GMT_program);
	}

	GMT_hash_init (localhashnode, operator, HASH_SIZE, N_OPERATORS);

	for (i = 0; i < STACK_SIZE; i++) {
		constant[i] = FALSE;
		factor[i] = 0.0;
		grd[i].nx = grd[i].ny = 0;
		stack[i] = (float *)NULL;
	}
		
	/* Get header from one file so we can allocate space */
	
	GMT_grd_init (&header, argc, argv, FALSE);

	for (arg = 1; nm == 0 && arg < argc; arg++) {
	
		if (argv[arg][0] == '-') continue;	/* Command line option or a negative number */
		if (decode_argument (argv[arg]) == ARG_IS_SAVE) {	/* Output case = */
			arg++;	/* Skip the output file name since the file does not exist yet */
			continue;
		}
		if (decode_argument (argv[arg]) != ARG_IS_FILE) continue;
		if (arg < (argc - 1) && !(strncmp (argv[arg+1], "LDIST", 5) && strncmp (argv[arg+1], "PDIST", 5) && strncmp (argv[arg+1], "INSIDE", 6))) continue;	/* Not a grd file */
		
		strcpy (file, argv[arg]);
		for (j = 0; file[j]; j++) if (file[j] == '=') file[j] = 0;
		if (GMT_read_grd_info (argv[arg], &header)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file);
			exit (EXIT_FAILURE);
		}
		
		nm = header.nx * header.ny;
	}
	
	/* Scan command line for -R, -I, -F, -f, -V */

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

			switch (argv[arg][1]) {

				case 'R':
					set_r = TRUE;
				case 'V':
				case 'f':
					error += GMT_get_common_args (argv[arg], &header.x_min, &header.x_max, &header.y_min, &header.y_max);
					break;

				case 'F':
					header.node_offset = 1;
					break;

				case 'I':
					GMT_getinc (&argv[arg][2], &header.x_inc, &header.y_inc);
					set_inc = (header.x_inc > 0 && header.y_inc > 0);
					break;
			}
		}
	}
	if (nm && set_r && set_inc) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot use -R, -I when grdfiles are specified\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (set_r && set_inc) {
		GMT_grd_RI_verify (&header, 1);
		header.nx = irint ((header.x_max - header.x_min) / header.x_inc) + !header.node_offset;
		header.ny = irint ((header.y_max - header.y_min) / header.y_inc) + !header.node_offset;
		nm = header.nx * header.ny;
	}
	
	if (nm == 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Expression must contain at least one grdfile or -R, -I\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	
	stack[0] = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);

	/* Get x and y vectors */

	grd_x = (float *) GMT_memory (VNULL, (size_t)header.nx, sizeof (float), GMT_program);
	grd_y = (float *) GMT_memory (VNULL, (size_t)header.ny, sizeof (float), GMT_program);
	d = (header.node_offset) ? 0.5 : 0.0;
	for (j = 0; j < header.ny; j++) grd_y[j] = (float)((j == (header.ny-1)) ? header.y_min + d * header.y_inc : header.y_max - (j + d) * header.y_inc);
	for (i = 0; i < header.nx; i++) grd_x[i] = (float)((i == (header.nx-1)) ? header.x_max - d * header.x_inc: header.x_min + (i + d) * header.x_inc);
	x_noise = SMALL * header.x_inc;	y_noise = SMALL * header.y_inc;

	grdmath_init (call_operator, consumed_operands, produced_operands);

	nstack = 0;
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: ", GMT_program);
	for (arg = 1; !error && arg < argc; arg++) {
	
		/* First check if we should skip optional arguments */

		for (k = 0, skip_arg = FALSE; !skip_arg && k < N_SKIP_ARGS; k++) skip_arg = !strncmp (argv[arg], skip_args[k], 2);
		if (skip_arg) continue;

		type = decode_argument (argv[arg]);
		if (type == ARG_IS_SAVE) {	/* Time to save the current stack to output and pop the stack */
			if (nstack <= 0 && n_items) {
				fprintf (stderr, "%s: GMT SYNTAX ERROR:  No items left on the stack for output!\n", GMT_program);
				exit (EXIT_FAILURE);
			}
			arg++;	/* Skip to the file name */
			outfile = argv[arg];
	
			if (gmtdefs.verbose) fprintf (stderr, "= %s", outfile);

			GMT_grd_init (&header, argc, argv, TRUE);
	
			if (n_items && new_stack < 0 && constant[nstack-1]) {	/* Only a constant provided, set grid accordingly */
				for (k = 0; k < nm; k++) stack[nstack-1][k] = (float)factor[nstack-1];
			}
			this_stack = (n_items) ? nstack - 1 : 0;	/* Since giving no args just go get an empty grid is OK */
			if (GMT_write_grd (outfile, &header, stack[this_stack], 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile);
				exit (EXIT_FAILURE);
			}
			if (n_items) nstack--;	/* Pop off the current stack if there is one */
			new_stack = nstack;
			continue;
		}
		if (type != ARG_IS_OPERATOR) {	/* File name or factor */
		
			if (type == ARG_IS_FILE && !(strncmp (argv[arg+1], "LDIST", 5) && strncmp (argv[arg+1], "PDIST", 5) && strncmp (argv[arg+1], "INSIDE", 6))) type = ARG_IS_ASCIIFILE;
			
			if (nstack == STACK_SIZE) {	/* Stack overflow */
				error = TRUE;
				continue;
			}
			n_items++;
			if (type == ARG_IS_NUMBER) {
				constant[nstack] = TRUE;
				ok = sscanf (argv[arg], "%lf", &factor[nstack]);
				error = !ok;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (type == ARG_IS_PI) {
				constant[nstack] = TRUE;
				factor[nstack] = M_PI;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (type == ARG_IS_E) {
				constant[nstack] = TRUE;
				factor[nstack] = M_E;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}

			/* Here we need a matrix */

			GMT_grd_init (&grd[nstack], argc, argv, TRUE);
			if (!stack[nstack]) stack[nstack] = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
			constant[nstack] = FALSE;

			if (type == ARG_IS_X_MATRIX) {		/* Need to set up matrix of x-values */
				if (gmtdefs.verbose) fprintf (stderr, "X ");
				for (j = k = 0; j < header.ny; j++, k += header.nx)
					memcpy ((void *)&stack[nstack][k], (void *)grd_x, (size_t)(header.nx * sizeof (float)));
			}
			else if (type == ARG_IS_Y_MATRIX) {	/* Need to set up matrix of y-values */
				if (gmtdefs.verbose) fprintf (stderr, "Y ");
				for (j = k = 0; j < header.ny; j++) for (i = 0; i < header.nx; i++, k++)
					stack[nstack][k] = grd_y[j];
			}
			else if (type == ARG_IS_ASCIIFILE) {
				strcpy (ASCII_file, argv[arg]);
				if (gmtdefs.verbose) fprintf (stderr, "(%s) ", argv[arg]);
			}
			else if (type == ARG_IS_FILE) {		/* Filename given */
				if (gmtdefs.verbose) fprintf (stderr, "%s ", argv[arg]);
				if (GMT_read_grd_info (argv[arg], &grd[nstack])) {
					fprintf (stderr, "%s: Error opening file %s\n", GMT_program, argv[arg]);
					exit (EXIT_FAILURE);
				}
				if (grd[nstack].nx != header.nx || grd[nstack].ny != header.ny) {
					fprintf (stderr, "%s: grd files not of same size!\n", GMT_program);
					exit (EXIT_FAILURE);
				}
				else if (fabs (grd[nstack].x_min - header.x_min) > x_noise || fabs (grd[nstack].x_max - header.x_max) > x_noise || 
					fabs (grd[nstack].y_min - header.y_min) > y_noise || fabs (grd[nstack].y_max - header.y_max) > y_noise) {
					fprintf (stderr, "%s: grd files do not cover the same area!\n", GMT_program);
					exit (EXIT_FAILURE);
				}
				if (GMT_read_grd (argv[arg], &grd[nstack], stack[nstack], 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
					fprintf (stderr, "%s: Error reading file %s\n", GMT_program, argv[arg]);
					exit (EXIT_FAILURE);
				}

			}
			nstack++;
			continue;
		}
		
		/* Here we have an operator */
		
		if ((op = get_operator (argv[arg])) < 0) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Unrecognized operator %s\n", GMT_program, argv[arg]);
			exit (EXIT_FAILURE);
		}

		if ((new_stack = nstack - consumed_operands[op] + produced_operands[op]) >= STACK_SIZE) {
			error = TRUE;
			continue;
		}

		if (nstack < consumed_operands[op]) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Operation \"%s\" requires %d operands\n", GMT_program, operator[op], consumed_operands[op]);
			exit (EXIT_FAILURE);
		}
		
		n_items++;
		if (gmtdefs.verbose) fprintf (stderr, "%s ", operator[op]);

		for (i = produced_operands[op] - consumed_operands[op]; i > 0; i--) {
			 if (stack[nstack+i-1])	continue;

			/* Must make space for more */

			stack[nstack+i-1] = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
		}

		/* If operators operates on constants only we may have to make space as well */

		for (j = 0, i = nstack - consumed_operands[op]; j < produced_operands[op]; j++, i++) {
			if (constant[i] && !stack[i]) stack[i] = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
		}

		(*call_operator[op]) (stack, constant, factor, nstack - 1);	/* Do it */

		nstack = new_stack;
	
		for (i = 1; i <= produced_operands[op]; i++)
			constant[nstack-i] = FALSE;	/* Now filled with grid */
	}
	
	if (error && !ok) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Unable to decode constant %s (File not found?)\n", GMT_program, argv[i-1]);
		exit (EXIT_FAILURE);
	}
	
	if (error) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Stack overflow (%s)\n", GMT_program, argv[i-1]);
		exit (EXIT_FAILURE);
	}
	
	for (i = 0; i < STACK_SIZE; i++) if (stack[i]) GMT_free ((void *)stack[i]);
	GMT_free ((void *)grd_x);
	GMT_free ((void *)grd_y);
	
	if (gmtdefs.verbose) fprintf (stderr, "\n");

	if (nstack > 1) fprintf (stderr, "%s: Warning: %d more operands left on the stack!\n", GMT_program, nstack-1);

	GMT_end (argc, argv);
}

/* -----------------------------------------------------------------
 *              Definitions of all operator functions
 * -----------------------------------------------------------------*/

void grd_ABS (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = fabs (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : fabs ((double)stack[last][i]));
}

void grd_ACOS (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) > 1.0) fprintf (stderr, "%s: Warning, |operand| > 1 for ACOS!\n", GMT_program);
	if (constant[last]) a = d_acos (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_acos ((double)stack[last][i]));
}

void grd_ACOSH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) > 1.0) fprintf (stderr, "%s: Warning, operand < 1 for ACOSH!\n", GMT_program);
	if (constant[last]) a = acosh (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : acosh ((double)stack[last][i]));
}

void grd_ADD (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a + b);
	}
}

void grd_AND (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)((GMT_is_dnan (a)) ? b : a);
	}
}

void grd_ASIN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) > 1.0) fprintf (stderr, "%s: Warning, |operand| > 1 for ASIN!\n", GMT_program);
	if (constant[last]) a = d_asin (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_asin ((double)stack[last][i]));
}

void grd_ASINH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = asinh (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : asinh ((double)stack[last][i]));
}

void grd_ATAN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = atan (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : atan ((double)stack[last][i]));
}

void grd_ATAN2 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for ATAN2!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for ATAN2!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)d_atan2 (a, b);
	}
}

void grd_ATANH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) >= 1.0) fprintf (stderr, "%s: Warning, |operand| >= 1 for ATANH!\n", GMT_program);
	if (constant[last]) a = atanh (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : atanh ((double)stack[last][i]));
}

void grd_BEI (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_bei (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_bei (fabs((double)stack[last][i])));
}

void grd_BER (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_ber (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_ber (fabs ((double)stack[last][i])));
}

void grd_CAZ (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, prev;
	double x, y, az;
	
	prev = last - 1;
	for (k = j = i = 0; k < nm; k++) {
		x = (constant[prev]) ? factor[prev] : stack[prev][k];
		y = (constant[last]) ? factor[last] : stack[last][k];
		az = 90.0 - (atan2 (y - (double)grd_y[j], x - (double)grd_x[i]) * R2D);
		while (az < -180.0) az += 360.0;
		while (az > +180.0) az -= 360.0;
		stack[prev][k] = (float)az;
		i++;
		if (i == header.nx) i = 0, j++;
	}
}

void grd_CDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, prev;
	double a, b;
	
	prev = last - 1;
	for (k = j = i = 0; k < nm; k++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][k];
		b = (constant[last]) ? factor[last] : stack[last][k];
		stack[prev][k] = (float)hypot (a - (double)grd_x[i], b - (double)grd_y[j]);
		i++;
		if (i == header.nx) i = 0, j++;
	}
}

void grd_CEIL (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = ceil (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : ceil ((double)stack[last][i]));
}

void grd_CHICRIT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for CHICRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for CHICRIT!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)GMT_chi2crit (a, b);
	}
}

void grd_CHIDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b, prob;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for CHIDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for CHIDIST!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		GMT_chi2 (a, b, &prob);
		stack[prev][i] = (float)prob;
	}
}

void grd_COS (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = cos (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : cos ((double)stack[last][i]));
}

void grd_COSD (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = cosd (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : cosd ((double)stack[last][i]));
}

void grd_COSH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = cosh (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : cosh ((double)stack[last][i]));
}

void grd_CURV (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, nx;
	double cx, cy;
	BOOLEAN easy;
	float *z;

	/* Curvature (Laplacian) */

	easy = (header.x_inc == header.y_inc);

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to CURV is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	
	z = (float *) (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
 
	nx = header.nx;
	cx = +0.5 / (header.x_inc * header.x_inc);
	cy = -0.5 / (header.y_inc * header.y_inc);	/* Because the loop over j below goes from ymax to ymin we compensate with a minus sign here */

	/* First left/right */

	for (j = 1, k = nx; j < header.ny-1; j++, k += nx)
		z[k] = (float)(cy * (stack[last][k+nx] - 2.0 * stack[last][k] + stack[last][k-nx]));
	for (j = 1, k = 2*nx-1; j < header.ny-1; j++, k += nx)
		z[k] = (float)(cy * (stack[last][k+nx] - 2.0 * stack[last][k] + stack[last][k-nx]));

	/* Then top/bottom */

	for (i = k = 1; i < header.nx - 1; i++, k++)
		z[k] = (float)(cy * (stack[last][k+1] - 2.0 * stack[last][k] + stack[last][k-1]));
	for (i = 1, k = nm - nx + 1; i < header.nx - 1; i++, k++)
		z[k] = (float)(cy * (stack[last][k+1] - 2.0 * stack[last][k] + stack[last][k-1]));

	/* Then inside */

	for (j = 1, k = nx; j < header.ny-1; j++) {
		k++;
		for (i = 1; i < header.nx-1; i++, k++) {
			if (easy)
				z[k] = (float)(cx * (stack[last][k+1] + stack[last][k-1] + stack[last][k+nx] + stack[last][k-nx] - 4.0 * stack[last][k]));
			else
				z[k] = (float)(cx * (stack[last][k+1] - 2.0 * stack[last][k] + stack[last][k-1]) + cy * (stack[last][k+nx] - 2 * stack[last][k] + stack[last][k-nx]));
		}
		k++;
	}

	memcpy ((void *)stack[last], (void *)z, (size_t)(nm * sizeof (float)));
	GMT_free ((void *)z);
}

void grd_D2DX2 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k;
	double c, left, next_left;

	/* Central 2nd difference in x */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to D2DX2 is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}

	c = 1.0 / (header.x_inc * header.x_inc);
	for (j = k = 0; j < header.ny; j++) {
		next_left = stack[last][k];
		stack[last][k++] = (float)0.0;
		for (i = 1; i < header.nx-1; i++, k++) {
			left = next_left;
			next_left = stack[last][k];
			stack[last][k] = (float)(c * (stack[last][k+1] - 2.0 * stack[last][k] + left));
		}
		stack[last][k++] = (float)0.0;
	}
}

void grd_D2DY2 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, nx;
	double c, bottom, next_bottom;

	/* Central 2nd difference in y */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to D2DY2 is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}

	c = -1.0 / (header.y_inc * header.y_inc);	/* Because the loop over j below goes from ymax to ymin we compensate with a minus sign here */
	nx = header.nx;
	for (i = 0; i < header.nx; i++) {
		k = i;
		next_bottom = stack[last][k];
		stack[last][k] = (float)0.0;
		k += nx;
		for (j = 1; j < header.ny - 1; j++, k += nx) {
			bottom = next_bottom;
			next_bottom = stack[last][k];
			stack[last][k] = (float)(c * (stack[last][k+nx] - 2 * stack[last][k] + bottom));
		}
		stack[last][k] = (float)(2.0 * c * (stack[last][k] - bottom));
	}
}

void grd_D2R (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = factor[last] * D2R;
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : (stack[last][i] * D2R));
}

void grd_DDX (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k;
	double c, left, next_left;

	/* Central 1st difference in x */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to DDX is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}

	c = 0.5 / header.x_inc;
	for (j = k = 0; j < header.ny; j++) {
		next_left = 2.0 * stack[last][k] - stack[last][k+1];
		for (i = 0; i < header.nx-1; i++, k++) {
			left = next_left;
			next_left = stack[last][k];
			stack[last][k] = (float)(c * (stack[last][k+1] - left));
		}
		stack[last][k] = (float)(2.0 * c * (stack[last][k] - next_left));
		k++;
	}
}

void grd_DDY (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, nx;
	double c, bottom, next_bottom;

	/* Central 1st difference in y */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to DDY is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}

	c = -0.5 / header.y_inc;	/* Because the loop over j below goes from ymax to ymin we compensate with a minus sign here */
	nx = header.nx;
	for (i = 0; i < header.nx; i++) {
		k = i;
		next_bottom = 2.0 * stack[last][k] - stack[last][k+nx];
		for (j = 0; j < header.ny - 1; j++, k += nx) {
			bottom = next_bottom;
			next_bottom = stack[last][k];
			stack[last][k] = (float)(c * (stack[last][k+nx] - bottom));
		}
		stack[last][k] = (float)(2.0 * c * (stack[last][k] - next_bottom));
	}
}

void grd_DILOG (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_dilog (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_dilog (stack[last][i]));
}

void grd_DIV (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	
	if (constant[last] && factor[last] == 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot divide by zero\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (constant[last]) {	/* Turn divide into multiply */
		a = factor[last];	/* Save original factor */
		factor[last] = 1.0 / factor[last];
		grd_MUL (stack, constant, factor, last);
		factor[last] = a;	/* Restore factor */
		return;
	}
	
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a / b);
	}
}

void grd_DUP (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int next, i;
	
	next = last + 1;
	constant[next] = constant[last];
	factor[next] = factor[last];
	if (constant[last]) {	/* Time to fess up */
		for (i = 0; i < nm; i++) stack[last][i] = (float)factor[last];
	}

	memcpy ((void *)stack[next], (void *)stack[last], (size_t)(nm * sizeof (float)));
}

void grd_ERF (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = erf (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : erf ((double)stack[last][i]));
}

void grd_ERFC (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = erfc (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : erfc ((double)stack[last][i]));
}

void grd_EQ (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a == b);
	}
}

void grd_ERFINV (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_erfinv (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_erfinv ((double)stack[last][i]));
}

void grd_EXCH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		if (constant[prev]) stack[prev][i] = (float)factor[prev];
		if (constant[last]) stack[last][i] = (float)factor[last];
		f_swap (stack[last][i], stack[prev][i]);
	}
	d_swap (factor[last], factor[prev]);
	i_swap (constant[last], constant[prev]);
}

void grd_EXP (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = exp (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : exp ((double)stack[last][i]));
}

void grd_EXTREMA (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, nx1, ny1, dx, dy, diag;
	int do_derivative (float *z, int this_node, int this_pos, int last_i, int last_j, int type);
	int do_diagonal (float *z, int this_node, int this_i, int this_j, int last_i, int last_j, int type);
	float *z;

	/* Find local extrema in grid */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to EXTREMA is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	
	z = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
 
	nx1 = header.nx - 1;
	ny1 = header.ny - 1;

	/* We will visit each node on the grid and determine if there are extrema.  We do this
	 * by looking at the along-x and along-y profiles separately.  If both of them shows the
	 * central node in a min or max location with respect to its two neighbors (one if at the
	 * edge of the grid or in the presence of NaNs) we must do further checking.  If the two
	 * extrema are of different sign we call it a saddle point and we are done.  If both are
	 * min or max then we look at the two diagonal lines to see if they can confirm this.
	 * Here, NaNs are not held agains you - it takes real values to overrule the dx,dy values.
	 *
	 * Min is given -2, Max is given +2, Saddle points are +1 (if max in x) or -1 (if max in y)
	 * Default is 0 which means no extrema.
	 */
	 
	for (j = k = 0; j < header.ny; j++) {	/* k is middle (current) node */
		for (i = 0; i < header.nx; i++, k++) {
		
			if (GMT_is_fnan (stack[last][k])) continue;	/* No extrema if point is NaN */
			
			if ((dx = do_derivative (stack[last], k, i, nx1, header.nx, 0)) == -2) continue;	/* Too many NaNs or flat x-line */
			if ((dy = do_derivative (stack[last], k, j, ny1, header.nx, 1)) == -2) continue;	/* Too many NaNs or flat y-line */
			
			if ((dx * dy) == 0) continue;	/* No min or max possible */
			if ((dx * dy) < 0) {	/* Saddle point - don't need to check diagonals */
				z[k] = (float)((dx > 0) ? +1 : -1);
				continue;
			}
			
			/* Need to examine diagonal trends to verify min or max */
			
			if ((diag = do_diagonal (stack[last], k, i, j, nx1, ny1, 0)) == -2) continue;	/* Sorry, no extrema along diagonal N45E */
			if (diag != 0 && diag != dx) continue;						/* Sorry, extrema of opposite sign along diagonal N45E  */
			if ((diag = do_diagonal (stack[last], k, i, j, nx1, ny1, 1)) == -2) continue;	/* Sorry, no extrema along diagonal N135E */
			if (diag != 0 && diag != dx) continue;						/* Sorry, extrema of opposite sign along diagonal N135E  */
			
			/* OK, we have a min or max point; just use dx to check which kind */
			
			z[k] = (float)((dx > 0) ? +2 : -2);
		}
	}

	memcpy ((void *)stack[last], (void *)z, (size_t)(nm * sizeof (float)));
	GMT_free ((void *)z);
}

/* Subroutines for grd_EXTREMA */

int do_derivative (float *z, int this_node, int this_pos, int last, int nx, int type)
{
	/* z is the data matrix */
	/* this_node is the index that gives us the center node */
	/* this_pos is either the x index (i) or the y index (j), depending on type */
	/* last_i is the last x index (nx-1) */
	/* last_j is the last y index (ny-1) */
	/* type = 0 means x-derivative, type = 1 means y-derivative */
	
	int off, next_node, prev_node;
	
	off = (type == 0) ? 1 : nx;
	
	if (this_pos == 0) {	/* Node at left column (x) or top row (y) - hence only one neighbor  - hence only one neighbor to the right (x) or below (y) */
		next_node = this_node + off;
		if (GMT_is_fnan (z[next_node])) return (-2);		/* One of the two points is a NaN */
		if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
		if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
		return (+1);						/* Else it must be a local maximum */
	}
	else if (this_pos == last) {	/* Node at right column (x) or bottom row (y) - hence only one neighbor to the left (x) or above (y) */
		prev_node = this_node - off;
		if (GMT_is_fnan (z[prev_node])) return (-2);		/* One of the two points is a NaN */
		if (z[this_node] == z[prev_node]) return (-2);		/* Flat line, no extrema possible */
		if (z[this_node] < z[prev_node]) return (-1);		/* A local minimum */
		return (+1);						/* Else it must be a local maximum */
	}
	else {	/* Node has neighbors on either side */
		prev_node = this_node - off;	/* Node to the left */
		next_node = this_node + off;	/* Node to the right */
		if (GMT_is_fnan (z[prev_node])) {			/* At least one of the two neighbor points is a NaN */
			if (GMT_is_fnan (z[next_node])) return (-2);	/* Both points are NaN */
			if (z[this_node] == z[next_node]) return (-2);	/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);	/* A local minimum */
			return (+1);					/* Else it must be a local maximum */
		}
		else if (GMT_is_fnan (z[next_node])) {			/* One of the two neighbor points is a NaN */
			if (z[this_node] == z[prev_node]) return (-2);	/* Flat line, no extrema possible */
			if (z[this_node] < z[prev_node]) return (-1);	/* A local minimum */
			return (+1);					/* Else it must be a local maximum */
		}
		else if (z[this_node] == z[prev_node] && z[this_node] == z[next_node])
			return (-2);					/* Flat line, no extrema possible */
		else if (z[this_node] < z[prev_node] && z[this_node] < z[next_node])
			return (-1);					/* A local minimum */
		else if (z[this_node] > z[prev_node] && z[this_node] > z[next_node])
			return (+1);					/* A local maximum */
		else
			return (0);					/* Nuthin' */
	}

	/* Should never get here... */
	
	fprintf (stderr, "GMT: Error in do_derivative - notify GMT developers\n");
	return (0);
}

int do_diagonal (float *z, int this_node, int this_i, int this_j, int last_i, int last_j, int type)
{
	/* Need to take d/dn in either the N45E direction (type = 0) or N135E (type = 1) direction */
	
	int off, next_node, prev_node;

	off = last_i + 1; 
	if (type == 0) off = -off;
	
	/* diagonal check is only fatal if there is actual data to rule out a min or max */
	/* If there are NaNs we do not rule out the min/max */
	
	if (this_j == 0) {	/* Node is somewhere along top row so limited checking only */
		
		if (this_i == 0) {	/* Node is upper left corner - can only do limited checking on one half-diagonal */
			next_node = this_node + last_i + 2;
			if (GMT_is_fnan (z[next_node])) return (0);		/* The only neighbor points is a NaN */
			if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
		else if (this_i == last_i) {	/* Node is upper right corner - can only do limited checking on one half-diagonal */
			prev_node = this_node + last_i;
			if (GMT_is_fnan (z[prev_node])) return (0);		/* The only neighbor points is a NaN */
			if (z[this_node] == z[prev_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[prev_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
		else {	/* Node is along top row but not at the corners - must check one of the two half-diagonals below */
			next_node = this_node + ((type) ? last_i + 2 : last_i);
			if (GMT_is_fnan (z[next_node])) return (0);		/* One of the two neighbor points is a NaN */
			if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
	}
	else if (this_j == last_j) {	/* Node is somewhere along bottom row so limited checking only */

		if (this_i == 0) {	/* Node is lower left corner - can only do limited checking on one half-diagonal */
			next_node = this_node - last_i;
			if (GMT_is_fnan (z[next_node])) return (0);		/* The only neighbor points is a NaN */
			if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
		else if (this_i == last_i) {	/* Node is lower right corner - can only do limited checking on one half-diagonal */
			prev_node = this_node - last_i - 2;
			if (GMT_is_fnan (z[prev_node])) return (0);		/* The only neighbor points is a NaN */
			if (z[this_node] == z[prev_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[prev_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
		else {	/* Node is along bottom row but not at the corners - must check one of the two half-diagonals below */
			next_node = this_node - ((type) ? last_i + 2 : last_i);
			if (GMT_is_fnan (z[next_node])) return (0);		/* One of the two neighbor points is a NaN */
			if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
	}
	else {	/* Node is along one of the intermediate rows - checking is limited by x-position only */
		
		if (this_i == 0) {	/* Node is on left row - can only do limited checking on two half-diagonals */
			next_node = this_node + off + 1;
			if (GMT_is_fnan (z[next_node])) return (0);		/* One of the two neighbor points is a NaN */
			if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
		else if (this_i == last_i) {	/* Node is on right row - can only do limited checking on two half-diagonals */
			next_node = this_node + off - 1;
			if (GMT_is_fnan (z[next_node])) return (0);		/* One of the two neighbor points is a NaN */
			if (z[this_node] == z[next_node]) return (-2);		/* Flat line, no extrema possible */
			if (z[this_node] < z[next_node]) return (-1);		/* A local minimum */
			return (+1);						/* Else it must be a local maximum */
		}
		else {	/* Node has all required neighbors */
			next_node = this_node + off + 1;
			prev_node = this_node - off - 1;
			if (GMT_is_fnan (z[prev_node])) {			/* At least one of the two neighbor points is a NaN */
				if (GMT_is_fnan (z[next_node])) return (0);	/* Both points are NaN */
				if (z[this_node] == z[next_node]) return (-2);	/* Flat line, no extrema possible */
				if (z[this_node] < z[next_node]) return (-1);	/* A local minimum */
				return (+1);					/* Else it must be a local maximum */
			}
			else if (GMT_is_fnan (z[next_node])) {			/* One of the two neighbor points is a NaN */
				if (z[this_node] == z[prev_node]) return (-2);	/* Flat line, no extrema possible */
				if (z[this_node] < z[prev_node]) return (-1);	/* A local minimum */
				return (+1);					/* Else it must be a local maximum */
			}
			else if (z[this_node] == z[prev_node] && z[this_node] == z[next_node])
				return (-2);					/* Flat line, no extrema possible */
			else if (z[this_node] < z[prev_node] && z[this_node] < z[next_node])
				return (-1);					/* A local minimum */
			else if (z[this_node] > z[prev_node] && z[this_node] > z[next_node])
				return (+1);					/* A local maximum */
			else
				return (0);					/* Nuthin' */
		}
	}

	/* Should never get here... */
	
	fprintf (stderr, "GMT: Error in do_diagonal - notify GMT developers\n");
	return (0);
}
		
void grd_FCRIT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, nu1, nu2, prev1, prev2;
	double alpha;
	
	prev1 = last - 1;
	prev2 = last - 2;
	if (gmtdefs.verbose && constant[prev2] && factor[prev2] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for FCRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev1] && factor[prev1] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for FCRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand three == 0 for FCRIT!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		alpha = (constant[prev2]) ? factor[prev2] : stack[prev2][i];
		nu1 = irint ((double)((constant[prev1]) ? factor[prev1] : stack[prev1][i]));
		nu2= irint ((double)((constant[last]) ? factor[last] : stack[last][i]));
		stack[prev2][i] = (float)GMT_Fcrit (alpha, (double)nu1, (double)nu2);
	}
}

void grd_FDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, nu1, nu2, prev1, prev2;
	double F, chisq1, chisq2 = 1.0, prob;
	
	prev1 = last - 1;
	prev2 = last - 2;
	if (gmtdefs.verbose && constant[prev1] && factor[prev1] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand three == 0 for FDIST!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		F = (constant[prev2]) ? factor[prev2] : stack[prev2][i];
		nu1 = irint ((double)((constant[prev1]) ? factor[prev1] : stack[prev1][i]));
		nu2 = irint ((double)((constant[last]) ? factor[last] : stack[last][i]));
		/* Since GMT_f_q needs chisq1 and chisq2, we set chisq2 = 1 and solve for chisq1 */
		chisq1 = F * nu1 / nu2;
		(void) GMT_f_q (chisq1, nu1, chisq2, nu2, &prob);
		stack[prev2][i] = (float)prob;
	}
}

void grd_FLIPLR (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, k0, nx1, nx_half;

	/* Reverse order of all rows */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to FLIPLR is constant!\n", GMT_program);
		return;
	}

	nx_half = header.nx / 2;
	nx1 = header.nx - 1;
	for (j = k0 = 0; j < header.ny; j++, k0 += header.nx) {	/* Do this to all rows */
		for (i = 0, k = nx1; i < nx_half; i++, k--) {
			f_swap (stack[last][k0+i], stack[last][k0+k]);
		}
	}
}

void grd_FLIPUD (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, ny1, ny_half;

	/* Reverse order of all columns */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to FLIPLR is constant!\n", GMT_program);
		return;
	}

	ny_half = header.ny / 2;
	ny1 = header.ny - 1;
	for (i = 0; i < header.nx; i++) {
		for (j = 0, k = ny1; j < ny_half; j++, k--) {	/* Do this to all rows */
			f_swap (stack[last][j*header.nx+i], stack[last][k*header.nx+i]);
		}
	}
}

void grd_FLOOR (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = floor (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : floor ((double)stack[last][i]));
}

void grd_FMOD (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)fmod (a, b);
	}
}

void grd_GE (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a >= b);
	}
}

void grd_GT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a > b);
	}
}

void grd_HYPOT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)hypot (a, b);
	}
}

void grd_I0 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;

	if (constant[last]) a = GMT_i0 (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_i0 ((double)stack[last][i]));
}

void grd_I1 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;

	if (constant[last]) a = GMT_i1 (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_i1 ((double)stack[last][i]));
}

void grd_IN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for IN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for IN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = GMT_in (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < nm; i++) {
		if (simple)
			stack[prev][i] = (float)b;
		else {
			if (!constant[last]) order = irint (fabs ((double)stack[last][i]));
			stack[last][i] = (float)GMT_in (order, fabs ((double)stack[prev][i]));
		}
	}
}

void grd_INRANGE (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev1, prev2;
	BOOLEAN inrange;
	float a, b, c;
	
	/* last is C */
	prev1 = last - 1;	/* This is B */
	prev2 = last - 2;	/* This is A */

	/* Set to 1 where B <= A <= C, 0 elsewhere, except where
	 * A, B, or C = NaN, in which case we set answer to NaN */
	 
	if (constant[prev2]) a = (float)factor[prev2];
	if (constant[prev1]) b = (float)factor[prev1];
	if (constant[last])  c = (float)factor[last];
	
	for (i = 0; i < nm; i++) {
		if (!constant[prev2]) a = stack[prev2][i];
		if (!constant[prev1]) b = stack[prev1][i];
		if (!constant[last])  c = stack[last][i];
	
		if (GMT_is_fnan (a) || GMT_is_fnan (b) || GMT_is_fnan (c)) {
			stack[prev2][i] = GMT_f_NaN;
			continue;
		}
		
		inrange = (b <= a && a <= c);
		stack[prev2][i] = (float)inrange;
	}
}

void grd_INSIDE (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, n_pols, P, sign_lat, inside;
	double xp, yp, r, sa, ca, *sin_x, *cos_x;
	BOOLEAN greenwich, polar;
	struct GMT_LINES *pol;
	
	if (GMT_io.in_col_type[0] & GMT_IS_GEO) {	/* Geographic data */
		greenwich = (header.x_min < 0.0 && header.x_max > 0.0);
	}
	else {
		greenwich = FALSE;
	}
	
	n_pols  = GMT_lines_init (ASCII_file, &pol, -9999.0, greenwich);	/* -9999 means closed polygons */

	for (P = 0, polar = FALSE; P < n_pols; P++) if (pol[P].polar) {
		/* Polygon is a polar cap, convert to simple polar (r-theta) coordinates instead */
		polar = TRUE;
		for (j = 0; j < pol[P].np; j++) {
			r = 90.0 - fabs (pol[P].lat[j]);
			sincos (pol[P].lon[j] * D2R, &sa, &ca);
			pol[P].lon[j] = r * ca;
			pol[P].lat[j] = r * sa;
		}
	}
	if (polar) {	/* We will need sincos of longitudes */
		sin_x = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof (double), GMT_program);
		cos_x = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof (double), GMT_program);
		for (i = 0; i < header.nx; i++) sincos (grd_x[i] * D2R, &sin_x[i], &cos_x[i]);
	}
	
	sign_lat = irint (copysign (1.0, grd_y[0]));
	for (i = j = k = 0; k < nm; k++) {	/* Visit each node */
		P = 0;
		inside = FALSE;
		while (P < n_pols && !inside) {
			if (pol[P].polar && sign_lat == pol[P].polar) {	/* Polygon encloses the pole on same hemisphere as point */
				r = 90.0 - grd_y[j];
				xp = r * cos_x[i];
				yp = r * sin_x[i];
			}
			else {	/* Regular polygon not containing a pole */
				xp = grd_x[i];
				yp = grd_y[j];
			}
			inside = (GMT_non_zero_winding (xp, yp, pol[P].lon, pol[P].lat, pol[P].np) >= 1);
			P++;
		}
		stack[last][k] = (float)((inside) ? 1.0 : 0.0);
		i++;
		if (i == header.nx) {	/* Go to next row */
			i = 0, j++;
			sign_lat = irint (copysign (1.0, grd_y[j]));
		}
	}
	
	/* Free memory used for pol */
	
	for (P = 0; P < n_pols; P++) {
		GMT_free ((void *)pol[P].lon);
		GMT_free ((void *)pol[P].lat);
	}
	GMT_free ((void *)pol);
	if (polar) {
		GMT_free ((void *)sin_x);
		GMT_free ((void *)cos_x);
	}
}

void grd_INV (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last] && factor[last] == 0.0) {
		fprintf (stderr, "%s: Error, Cannot take inverse of zero!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (constant[last]) factor[last] = 1.0 / factor[last];
	for (i = 0; i < nm; i++) {
		a = (constant[last]) ? factor[last] : 1.0 / stack[last][i];
		stack[last][i] = (float)a;
	}
}

void grd_ISNAN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_is_dnan (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_is_fnan (stack[last][i]));
}

void grd_J0 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;

	if (constant[last]) a = j0 (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : j0 ((double)stack[last][i]));
}

void grd_J1 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = j1 (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : j1 (fabs ((double)stack[last][i])));
}

void grd_JN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for JN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for JN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = jn (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < nm; i++) {
		if (simple)
			stack[prev][i] = (float)b;
		else {
			if (!constant[last]) order = irint (fabs ((double)stack[last][i]));
			stack[last][i] = (float)jn (order, fabs ((double)stack[prev][i]));
		}
	}
}

void grd_K0 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;

	if (constant[last]) a = GMT_k0 (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_k0 ((double)stack[last][i]));
}

void grd_K1 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;

	if (constant[last]) a = GMT_k1 (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_k1 ((double)stack[last][i]));
}

void grd_KEI (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_kei (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_kei (fabs ((double)stack[last][i])));
}

void grd_KER (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_ker (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_ker (fabs ((double)stack[last][i])));
}

void grd_KN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for KN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for KN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = GMT_kn (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < nm; i++) {
		if (simple)
			stack[prev][i] = (float)b;
		else {
			if (!constant[last]) order = irint (fabs ((double)stack[last][i]));
			stack[last][i] = (float)GMT_kn (order, fabs ((double)stack[prev][i]));
		}
	}
}

void grd_LDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, n_lines;
	double d;
	BOOLEAN greenwich;
	PFI near_a_line;
	struct GMT_LINES *line;
	
	if (GMT_io.in_col_type[0] & GMT_IS_GEO) {	/* Geographic data */
		GMT_distance_func = (PFD) GMT_km_dist;
		near_a_line = (PFI) GMT_near_a_line_spherical;
		greenwich = (header.x_min < 0.0 && header.x_max > 0.0);
	}
	else {
		GMT_distance_func = (PFD) GMT_cartesian_dist;
		near_a_line = (PFI) GMT_near_a_line_cartesian;
		greenwich = FALSE;
	}
	
	n_lines  = GMT_lines_init (ASCII_file, &line, 0.0, greenwich);

	for (i = j = k = 0; k < nm; k++) {	/* Visit each node */
		(void) near_a_line (grd_x[i], grd_y[j], line, n_lines, TRUE, &d);
		stack[last][k] = (float)d;
		i++;
		if (i == header.nx) i = 0, j++;
	}
	
	/* Free memory used for line */
	
	for (k = 0; k < n_lines; k++) {
		GMT_free ((void *)line[k].lon);
		GMT_free ((void *)line[k].lat);
	}
	GMT_free ((void *)line);
}

void grd_LE (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a <= b);
	}
}

void grd_LOG (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, argument to log = 0\n", GMT_program);

	if (constant[last]) a = d_log (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_log (fabs ((double)stack[last][i])));
}

void grd_LOG10 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, argument to log10 = 0\n", GMT_program);

	if (constant[last]) a = d_log10 (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_log10 (fabs ((double)stack[last][i])));
}

void grd_LOG1P (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, argument to log1p < 0\n", GMT_program);

	if (constant[last]) a = d_log1p (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_log1p (fabs ((double)stack[last][i])));
}

void grd_LOG2 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, argument to log2 = 0\n", GMT_program);

	if (constant[last]) a = d_log (fabs (factor[last])) * M_LN2_INV;
	for (i = 0; i < nm; i++) stack[last][i] = (float)(((constant[last]) ? a : d_log (fabs ((double)stack[last][i]))) * M_LN2_INV);
}

void grd_LMSSCL (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double mode, lmsscl;
	float lmsscl_f;
	
	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	
	/* Sort will put any NaNs to the end - we then count to find the real data */
	
	qsort ((void *)stack[last], (size_t)nm, sizeof (float), GMT_comp_float_asc);
	for (i = nm; GMT_is_fnan (stack[last][i-1]) && i > 1; i--);
	if (i) {
		GMT_mode_f (stack[last], i, i/2, 0, GMT_mode_selection, &GMT_n_multiples, &mode);
		GMT_getmad_f (stack[last], i, mode, &lmsscl);
		lmsscl_f = (float)lmsscl;
	}
	else
		lmsscl_f = GMT_f_NaN;

	for (i = 0; i < nm; i++) stack[last][i] = lmsscl_f;
	if (GMT_n_multiples > 0) {
		fprintf (stderr, "%s: WARNING: %d Multiple modes found\n", GMT_program, GMT_n_multiples);
		GMT_n_multiples = 0;
	}
}

void grd_LOWER (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	float low;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < nm; i++) stack[last][i] = (float)factor[last];
		return;
	}
	
	for (i = 0, low = FLT_MAX; i < nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		if (stack[last][i] < low) low = stack[last][i];
	}
	for (i = 0; i < nm; i++) if (!GMT_is_fnan (stack[last][i])) stack[last][i] = low;
}

void grd_LRAND (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < nm; i++) {
		if (!constant[prev]) a = (double)stack[prev][i];
		if (!constant[last]) b = (double)stack[last][i];
		stack[prev][i] = (float)(a + b * GMT_lrand ());
	}
}

void grd_LT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a < b);
	}
}

void grd_MAD (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double mad, med;
	float mad_f;
	
	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	
	/* Sort will put any NaNs to the end - we then count to find the real data */
	
	qsort ((void *)stack[last], (size_t)nm, sizeof (float), GMT_comp_float_asc);
	for (i = nm; GMT_is_fnan (stack[last][i-1]) && i > 1; i--);
	if (i) {
		med = (i%2) ? stack[last][i/2] : (float)(0.5 * (stack[last][(i-1)/2] + stack[last][i/2]));
		GMT_getmad_f (stack[last], i, med, &mad);
		mad_f = (float)mad;
	}
	else
		mad_f = GMT_f_NaN;

	for (i = 0; i < nm; i++) stack[last][i] = mad_f;
}

void grd_MAX (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)MAX (a, b);
	}
}

void grd_MEAN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, n_a = 0;
	double sum_a = 0.0;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < nm; i++) stack[last][i] = (float)factor[last];
		return;
	}
	
	for (i = 0; i < nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		sum_a += stack[last][i];
		n_a++;
	}
	sum_a = (n_a) ? sum_a / n_a : 0.0;
	for (i = 0; i < nm; i++) stack[last][i] = (float)sum_a;
}

void grd_MED (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	float med;
	
	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < nm; i++) stack[last][i] = (float)factor[last];
		return;
	}

	qsort ((void *)stack[last], (size_t)nm, sizeof (float), GMT_comp_float_asc);
	for (i = nm; GMT_is_fnan (stack[last][i-1]) && i > 1; i--);
	if (i)
		med = (i%2) ? stack[last][i/2] : (float)(0.5 * (stack[last][(i-1)/2] + stack[last][i/2]));
	else
		med = GMT_f_NaN;

	for (i = 0; i < nm; i++) stack[last][i] = med;
}

void grd_MIN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)MIN (a, b);
	}
}

void grd_MODE (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double mode = 0.0;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < nm; i++) stack[last][i] = (float)factor[last];
		return;
	}
	
	qsort ((void *)stack[last], (size_t)nm, sizeof (float), GMT_comp_float_asc);
	for (i = nm; GMT_is_fnan (stack[last][i-1]) && i > 1; i--);
	if (i)
		GMT_mode_f (stack[last], i, i/2, 0, GMT_mode_selection, &GMT_n_multiples, &mode);
	else
		mode = GMT_f_NaN;

	for (i = 0; i < nm; i++) stack[last][i] = (float)mode;
	if (GMT_n_multiples > 0) {
		fprintf (stderr, "%s: WARNING: %d Multiple modes found\n", GMT_program, GMT_n_multiples);
		GMT_n_multiples = 0;
	}
}

void grd_MUL (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a * b);
	}
}

void grd_NAN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < nm; i++) {
		if (!constant[prev]) a = stack[prev][i];
		if (!constant[last]) b = stack[last][i];
		stack[prev][i] = ((a == b) ? GMT_f_NaN : (float)a);
	}
}

void grd_NEG (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = -factor[last];
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : -stack[last][i]);
}

void grd_NEQ (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a != b);
	}
}

void grd_NRAND (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < nm; i++) {
		if (!constant[prev]) a = (double)stack[prev][i];
		if (!constant[last]) b = (double)stack[last][i];
		stack[prev][i] = (float)(a + b * GMT_nrand ());
	}
}

void grd_OR (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)((GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_f_NaN : a);
	}
}

void grd_PDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, np, dummy;
	double *xx, *yy, *dd, d;
	BOOLEAN greenwich;
	
	if (GMT_io.in_col_type[0] & GMT_IS_GEO) {	/* Geographic data */
		GMT_distance_func = (PFD) GMT_km_dist;
		greenwich = (header.x_min < 0.0 && header.x_max > 0.0);
	}
	else {
		GMT_distance_func = (PFD) GMT_cartesian_dist;
		greenwich = FALSE;
	}
	
	np = GMT_points_init (ASCII_file, &xx, &yy, &dd, 0.0, greenwich);
	for (k = 0; k < np; k++) dd[k] = 0.0;	/* In order to get min dist returned */
	
	for (i = j = k = 0; k < nm; k++) {	/* Visit each node */
		d = GMT_dist_to_point (grd_x[i], grd_y[j], xx, yy, np, &dummy);
		stack[last][k] = (float)d;
		i++;
		if (i == header.nx) i = 0, j++;
	}
	
	/* Free memory used for points */
	
	GMT_free ((void *)xx);
	GMT_free ((void *)yy);
	GMT_free ((void *)dd);
}

void grd_POP (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	
	/* Dummy routine that does nothing but consume the top element of stack */
}

void grd_PLM (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev, first, L, M;
	double a;
				/* last holds the order M */
	prev  = last - 1;	/* prev holds the degree L */
	first = prev - 1;	/* first holds the argument x = cos(colat) */

	if (!(constant[prev] && constant[last])) {
		fprintf (stderr, "%s: L and M must be constants in PLM!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	L = irint (factor[prev]);
	M = irint (factor[last]);
	if (gmtdefs.verbose && constant[first] && (factor[first] < -1.0 || factor[first] > 1.0)) fprintf (stderr, "%s: Warning, argument to PLM outside domain!\n", GMT_program);

	if (constant[first]) a = GMT_plm (L, M, factor[first]);
	for (i = 0; i < nm; i++) stack[first][i] = (float)((constant[first]) ? a : GMT_plm (L, M, stack[first][i]));
}

void grd_POW (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;

	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)pow (a, b);
	}
}

void grd_R2 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	if (constant[prev]) factor[prev] *= factor[prev];
	if (constant[last]) factor[last] *= factor[last];
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i] * stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i] * stack[last][i];
		stack[prev][i] = (float)(a + b);
	}
}

void grd_R2D (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = R2D * factor[last];
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : R2D * stack[last][i]);
}

void grd_RAND (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < nm; i++) {
		if (!constant[prev]) a = (double)stack[prev][i];
		if (!constant[last]) b = (double)stack[last][i];
		stack[prev][i] = (float)(a + GMT_rand () * (b - a));
	}
}

void grd_RINT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = rint (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : rint ((double)stack[last][i]));
}

void grd_ROTX (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, prev, shift, *node;
	float *z;
	
	/* Shift grid A by the x-shift B.  B must be a constant */
	
	prev = last - 1;
	if (!constant[last]) {
		fprintf (stderr, "%s: DX shift (B) must be a constant in ROTX!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	shift = irint (factor[last] / header.x_inc);	/* Shift of nodes */
	
	if (constant[prev] || !shift) return;	/* Trivial since A is a constant or shift is zero */
	if (shift < 0) shift += header.nx;	/* Same thing */
	/* Set up permutation vector */
	
	node = (int *) GMT_memory (VNULL, (size_t)header.nx, sizeof (int), GMT_program);
	z =  (float *) GMT_memory (VNULL, (size_t)header.nx, sizeof (float), GMT_program);
	for (i = 0; i < header.nx; i++) node[i] = (i + shift) % header.nx;	/* Move by shift but rotate around */
	for (j = k = 0; j < header.ny; j++, k += header.nx) {	/* For each row */
		for (i = 0; i < header.nx; i++) z[node[i]] = stack[prev][k+i];
		memcpy ((void *)&stack[prev][k], (void *)z, (size_t)(header.nx * sizeof (float)));
	}
	GMT_free ((void *)z);
	GMT_free ((void *)node);
}

void grd_ROTY (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev, shift, *node;
	float *z;
	
	/* Shift grid A by the y-shift B.  B must be a constant */
	
	prev = last - 1;
	if (!constant[last]) {
		fprintf (stderr, "%s: DY shift (B) must be a constant in ROTY!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	shift = irint (factor[last] / header.y_inc);	/* Shift of nodes */
	
	if (constant[prev] || !shift) return;	/* Trivial since A is a constant or shift is zero */
	if (shift < 0) shift += header.ny;	/* Same thing */
	/* Set up permutation vector */
	
	node = (int *) GMT_memory (VNULL, (size_t)header.ny, sizeof (int), GMT_program);
	z =  (float *) GMT_memory (VNULL, (size_t)header.ny, sizeof (float), GMT_program);
	for (j = 0; j < header.ny; j++) node[j] = (j + header.ny - shift) % header.ny;	/* Move by shift but rotate around */
	for (i = 0; i < header.nx; i++) {	/* For each column */
		for (j = 0; j < header.ny; j++) z[node[j]] = stack[prev][j*header.nx+i];
		for (j = 0; j < header.ny; j++) stack[prev][j*header.nx+i] = z[j];
	}
	GMT_free ((void *)z);
	GMT_free ((void *)node);
}

void grd_SDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, prev;
	double a, b;
	
	prev = last - 1;
	for (k = j = i = 0; k < nm; k++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][k];
		b = (constant[last]) ? factor[last] : stack[last][k];
		stack[prev][k] = (float)GMT_great_circle_dist (a, b, (double)grd_x[i], (double)grd_y[j]);
		i++;
		if (i == header.nx) i = 0, j++;
	}
}

void grd_SAZ (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, prev;
	double x0, y0, sin_y0, cos_y0, sin_dlon, cos_dlon, az, *cy, *sy;
	
	prev = last - 1;
	if (constant[prev]) x0 = factor[prev];
	if (constant[last]) {
		y0 = factor[last];
		sincos (y0 * D2R, &sin_y0, &cos_y0);
	}
	cy = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof (double), GMT_program);
	sy = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof (double), GMT_program);
	for (j = 0; j < header.ny; j++) sincos (grd_y[j] * D2R, &sy[j], &cy[j]);
	for (k = j = i = 0; k < nm; k++) {
		if (!constant[prev]) x0 = (double)stack[prev][i];
		if (!constant[last]) {
			y0 = (double)stack[last][i];
			sincos (y0 * D2R, &sin_y0, &cos_y0);
		}
		sincos ((x0 - grd_x[i]) * D2R, &sin_dlon, &cos_dlon);
		az = (float)(atan2 (cos_y0 * sin_dlon, cy[j] * sin_y0 - sy[j] * cos_y0 * cos_dlon) * R2D);
		while (az < -180.0) az += 360.0;
		while (az > +180.0) az -= 360.0;
		stack[prev][k] = (float)az;
		i++;
		if (i == header.nx) i = 0, j++;
	}
	GMT_free ((void *)cy);
	GMT_free ((void *)sy);
}

void grd_SIGN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = copysign (1.0, factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : copysign (1.0, (double)stack[last][i]));
}

void grd_SIN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = sin (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : sin ((double)stack[last][i]));
}

void grd_SINC (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_sinc (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_sinc ((double)stack[last][i]));
}

void grd_SIND (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = sind (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : sind ((double)stack[last][i]));
}

void grd_SINH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = sinh (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : sinh ((double)stack[last][i]));
}

void grd_SQRT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (constant[last]) a = sqrt (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : sqrt ((double)stack[last][i]));
}

void grd_STD (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, n_a = 0;
	double sum_a = 0.0, sum_a2 = 0.0, std;
	
	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}

	for (i = 0; i < nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		sum_a += stack[last][i];
		sum_a2 += (stack[last][i] * stack[last][i]);
		n_a++;
	}
	if (n_a > 1)
		std = sqrt ((n_a * sum_a2 - sum_a * sum_a) / (n_a * (n_a - 1.0)));
	else
		std = 0.0;
	for (i = 0; i < nm; i++) stack[last][i] = (float)std;
}

void grd_STEP (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = factor[last];
	for (i = 0; i < nm; i++) {
		if (!constant[last]) a = (double)stack[last][i];
		if (a == 0.0)
			stack[last][i] = (float)0.5;
		else
			stack[last][i] = (float)((a < 0.0) ? 0.0 : 1.0);
	}
}

void grd_STEPX (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	for (i = 0; i < nm; i++) {
		a = grd_x[i%header.nx] - ((constant[last]) ? factor[last] : stack[last][i]);
		if (a == 0.0)
			stack[last][i] = (float)0.5;
		else
			stack[last][i] = (float)((a < 0.0) ? 0.0 : 1.0);
	}
}

void grd_STEPY (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	for (i = 0; i < nm; i++) {
		a = grd_y[i/header.nx] - ((constant[last]) ? factor[last] : stack[last][i]);
		if (a == 0.0)
			stack[last][i] = (float)0.5;
		else
			stack[last][i] = (float)((a < 0.0) ? 0.0 : 1.0);
	}
}

void grd_SUB (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (float)(a - b);
	}
}

void grd_TAN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = tan (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : tan ((double)stack[last][i]));
}

void grd_TAND (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = tand (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : tand ((double)stack[last][i]));
}

void grd_TANH (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = tanh (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : tanh ((double)stack[last][i]));
}

void grd_TCRIT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, b, prev;
	double a;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for TCRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for TCRIT!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = irint ((double)((constant[last]) ? factor[last] : stack[last][i]));
		stack[prev][i] = (float)GMT_tcrit (a, b);
	}
}

void grd_TDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, b, prev;
	double a, prob;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for TDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for TDIST!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = irint ((double)((constant[last]) ? factor[last] : stack[last][i]));
		(void) GMT_student_t_a (a, b, &prob);
		stack[prev][i] = (float)prob;
	}
}

void grd_UPPER (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	float high;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < nm; i++) stack[last][i] = (float)factor[last];
		return;
	}
	
	for (i = 0, high = -FLT_MAX; i < nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		if (stack[last][i] > high) high = stack[last][i];
	}
	for (i = 0; i < nm; i++) if (!GMT_is_fnan (stack[last][i])) stack[last][i] = high;
}

void grd_XOR (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < nm; i++) {
		if (!constant[prev]) a = stack[prev][i];
		if (!constant[last]) b = stack[last][i];
		stack[prev][i] = (float)((GMT_is_dnan (a)) ? b : a);
	}
}

void grd_Y0 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = y0 (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : y0 (fabs ((double)stack[last][i])));
}

void grd_Y1 (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = y1 (fabs (factor[last]));
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : y1 (fabs ((double)stack[last][i])));
}

void grd_YLM (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, prev, L, M;
	double x, z, normal, P, *C, *S, GMT_factorial(int n);

	/* Fully normalized surface harmonics on the -R grid given L and M
	 * -R assumed to be in degrees, returns both Re and Im components
	 * as two separate grids */

	prev = last - 1;
	if (!(constant[prev] && constant[last])) {
		fprintf (stderr, "%s: YLM takes constants L and M only!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	L = irint (factor[prev]);
	M = irint (factor[last]);

	if (M > 0) {	/* Need to make the Sin/Cos expression */

		C = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof (double), GMT_program);
		S = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof (double), GMT_program);
		z = M * D2R;
		for (i = 0; i < header.nx; i++) {
			x = z * grd_x[i];
			sincos (x, &S[i], &C[i]);
		}
		normal = pow (-1.0, (double) M) * d_sqrt ((2.0 * L + 1.0) * GMT_factorial (L - M) / (4.0 * M_PI * GMT_factorial (L + M)));
	}
	else
		normal = d_sqrt ((2.0 * L + 1.0) / (4.0 * M_PI ));


	for (k = j = 0; j < header.ny; j++) {	/* For each latitude */

		x = cosd (90.0 - grd_y[j]);	/* plm takes cos(colatitude) */
		P = normal * GMT_plm (L, M, x);
		if (M > 0) {
			for (i = 0; i < header.nx; i++, k++) {
				stack[prev][k] = (float)(P * C[i]);
				stack[last][k] = (float)(P * S[i]);
			}
		}
		else {
			for (i = 0; i < header.nx; i++, k++) {
				stack[prev][k] = stack[last][k] = (float)P;
			}
		}
	}

	if (M > 0) {
		GMT_free ((void *)C);
		GMT_free ((void *)S);
	}
}

void grd_YN (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, prev, order;
	double b;
	BOOLEAN simple = FALSE;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for YN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && (rint(factor[last]) != factor[last])) fprintf (stderr, "%s: Warning, order not an integer for YN!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, argument = 0 for YN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = yn (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < nm; i++) {
		if (simple)
			stack[prev][i] = (float)b;
		else {
			if (!constant[last]) order = irint (fabs ((double)stack[last][i]));
			stack[last][i] = (float)yn (order, fabs ((double)stack[prev][i]));
		}
	}
}

void grd_ZCRIT (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_zcrit (factor[last]);
	for (i = 0; i < nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_zcrit (a));
}

/* ---------------------- end operator functions --------------------- */

int decode_argument (txt)
char *txt; {
	int j, number, minus, plus, exp, dec, n_digits = 0;
	char *s, file[BUFSIZ];
	
	if (!strcmp (txt, "=")) return ARG_IS_SAVE;	/* Time to save stack; next arg is filename */
	if (!strcmp (txt, "/")) return ARG_IS_OPERATOR;	/* Special test since / is also a directory */

	/* First see if argument it can be opened as a file */
	
	strcpy (file, txt);
	for (j = 0; file[j]; j++) if (file[j] == '=') file[j] = 0;
	if (!access (file, R_OK)) return ARG_IS_FILE; 	/* returns FALSE if file exists and can be read */

	/* Next look for symbols with special meaning */

	if (!(strcmp (txt, "PI") && strcmp (txt, "pi"))) return ARG_IS_PI;
	if (!(strcmp (txt, "E") && strcmp (txt, "e"))) return ARG_IS_E;
	if (!strcmp (txt, "X")) return ARG_IS_X_MATRIX;
	if (!(strcmp (txt, "Y") && strcmp (txt, "y"))) return ARG_IS_Y_MATRIX;

	/* Here we must check if argument is a numerical value */

	s = txt;
	if (*s == '-' || *s == '+') s++;	/* Skip leading sign */
	
	minus = plus = exp = dec = 0;
	number = TRUE;

	while (number && *s) {
		if (isdigit ((int)*s))
			n_digits++;
		else {
			switch (*s) {
				case '-':
					minus++;
					break;
				case '+':
					plus++;
					break;
				case 'E':
				case 'e':
					exp++;
					break;
				case '.':
					dec++;
					break;
				default:
					number = FALSE;
					break;
			}
		}
		if (minus > 1 || exp > 1 || dec > 1) number = FALSE;
		s++;
	}
	if (txt[0] == '-' && !number) {	/* Probably a bad commandline option */
		fprintf (stderr, "%s: ERROR: Option %s not recognized\n", GMT_program, txt);
		exit (EXIT_FAILURE);
	}
	
	return ((number && n_digits > 0) ? ARG_IS_NUMBER : ARG_IS_OPERATOR);
	
}

int get_operator (choice)
char *choice; {
	int op;
	
	/* Returns -1 if not a registered operator */

	if (!strncmp (choice, "GDIST", 5)) choice[0] = 'S';	/* Changing GDIST to SDIST for backwards compatibility */
	
	op = GMT_hash_lookup (choice, localhashnode, HASH_SIZE, HASH_SIZE);

	if (op < 0 && strlen (choice) == 1) {	/* Check for old-style operators */
	
		switch (choice[0]) {
			case '+':
				op = ADD;
				break;
			case '-':
				op = SUB;
				break;
			case 'x':
				op = MUL;
				break;
			case '/':
				op = DIV;
				break;
			case '^':
				op = RAISE;
				break;
		}
	}

	return (op);
}

void grdmath_init (PFV ops[], int n_args[], int n_out[])
{

	/* Operator function		# of operands  		# of outputs */

	ops[0]=grd_ABS;		n_args[0]=1;		n_out[0]=1;
	ops[1]=grd_ACOS;		n_args[1]=1;		n_out[1]=1;
	ops[2]=grd_ACOSH;		n_args[2]=1;		n_out[2]=1;
	ops[3]=grd_ADD;		n_args[3]=2;		n_out[3]=1;
	ops[4]=grd_AND;		n_args[4]=2;		n_out[4]=1;
	ops[5]=grd_ASIN;		n_args[5]=1;		n_out[5]=1;
	ops[6]=grd_ASINH;		n_args[6]=1;		n_out[6]=1;
	ops[7]=grd_ATAN;		n_args[7]=1;		n_out[7]=1;
	ops[8]=grd_ATAN2;		n_args[8]=2;		n_out[8]=1;
	ops[9]=grd_ATANH;		n_args[9]=1;		n_out[9]=1;
	ops[10]=grd_BEI;		n_args[10]=1;		n_out[10]=1;
	ops[11]=grd_BER;		n_args[11]=1;		n_out[11]=1;
	ops[12]=grd_CAZ;		n_args[12]=2;		n_out[12]=1;
	ops[13]=grd_CDIST;		n_args[13]=2;		n_out[13]=1;
	ops[14]=grd_CEIL;		n_args[14]=1;		n_out[14]=1;
	ops[15]=grd_CHICRIT;		n_args[15]=2;		n_out[15]=1;
	ops[16]=grd_CHIDIST;		n_args[16]=2;		n_out[16]=1;
	ops[17]=grd_COS;		n_args[17]=1;		n_out[17]=1;
	ops[18]=grd_COSD;		n_args[18]=1;		n_out[18]=1;
	ops[19]=grd_COSH;		n_args[19]=1;		n_out[19]=1;
	ops[20]=grd_CURV;		n_args[20]=1;		n_out[20]=1;
	ops[21]=grd_D2DX2;		n_args[21]=1;		n_out[21]=1;
	ops[22]=grd_D2DY2;		n_args[22]=1;		n_out[22]=1;
	ops[23]=grd_D2R;		n_args[23]=1;		n_out[23]=1;
	ops[24]=grd_DDX;		n_args[24]=1;		n_out[24]=1;
	ops[25]=grd_DDY;		n_args[25]=1;		n_out[25]=1;
	ops[26]=grd_DILOG;		n_args[26]=1;		n_out[26]=1;
	ops[27]=grd_DIV;		n_args[27]=2;		n_out[27]=1;
	ops[28]=grd_DUP;		n_args[28]=1;		n_out[28]=2;
	ops[29]=grd_ERF;		n_args[29]=1;		n_out[29]=1;
	ops[30]=grd_ERFC;		n_args[30]=1;		n_out[30]=1;
	ops[31]=grd_ERFINV;		n_args[31]=1;		n_out[31]=1;
	ops[32]=grd_EQ;		n_args[32]=2;		n_out[32]=1;
	ops[33]=grd_EXCH;		n_args[33]=2;		n_out[33]=2;
	ops[34]=grd_EXP;		n_args[34]=1;		n_out[34]=1;
	ops[35]=grd_EXTREMA;		n_args[35]=1;		n_out[35]=1;
	ops[36]=grd_FCRIT;		n_args[36]=3;		n_out[36]=1;
	ops[37]=grd_FDIST;		n_args[37]=3;		n_out[37]=1;
	ops[38]=grd_FLIPLR;		n_args[38]=1;		n_out[38]=1;
	ops[39]=grd_FLIPUD;		n_args[39]=1;		n_out[39]=1;
	ops[40]=grd_FLOOR;		n_args[40]=1;		n_out[40]=1;
	ops[41]=grd_FMOD;		n_args[41]=2;		n_out[41]=1;
	ops[42]=grd_GE;		n_args[42]=2;		n_out[42]=1;
	ops[43]=grd_GT;		n_args[43]=2;		n_out[43]=1;
	ops[44]=grd_HYPOT;		n_args[44]=2;		n_out[44]=1;
	ops[45]=grd_I0;		n_args[45]=1;		n_out[45]=1;
	ops[46]=grd_I1;		n_args[46]=1;		n_out[46]=1;
	ops[47]=grd_IN;		n_args[47]=2;		n_out[47]=1;
	ops[48]=grd_INRANGE;		n_args[48]=3;		n_out[48]=1;
	ops[49]=grd_INSIDE;		n_args[49]=1;		n_out[49]=1;
	ops[50]=grd_INV;		n_args[50]=1;		n_out[50]=1;
	ops[51]=grd_ISNAN;		n_args[51]=1;		n_out[51]=1;
	ops[52]=grd_J0;		n_args[52]=1;		n_out[52]=1;
	ops[53]=grd_J1;		n_args[53]=1;		n_out[53]=1;
	ops[54]=grd_JN;		n_args[54]=2;		n_out[54]=1;
	ops[55]=grd_K0;		n_args[55]=1;		n_out[55]=1;
	ops[56]=grd_K1;		n_args[56]=1;		n_out[56]=1;
	ops[57]=grd_KN;		n_args[57]=2;		n_out[57]=1;
	ops[58]=grd_KEI;		n_args[58]=1;		n_out[58]=1;
	ops[59]=grd_KER;		n_args[59]=1;		n_out[59]=1;
	ops[60]=grd_LDIST;		n_args[60]=1;		n_out[60]=1;
	ops[61]=grd_LE;		n_args[61]=2;		n_out[61]=1;
	ops[62]=grd_LMSSCL;		n_args[62]=1;		n_out[62]=1;
	ops[63]=grd_LOG;		n_args[63]=1;		n_out[63]=1;
	ops[64]=grd_LOG10;		n_args[64]=1;		n_out[64]=1;
	ops[65]=grd_LOG1P;		n_args[65]=1;		n_out[65]=1;
	ops[66]=grd_LOG2;		n_args[66]=1;		n_out[66]=1;
	ops[67]=grd_LOWER;		n_args[67]=1;		n_out[67]=1;
	ops[68]=grd_LRAND;		n_args[68]=2;		n_out[68]=1;
	ops[69]=grd_LT;		n_args[69]=2;		n_out[69]=1;
	ops[70]=grd_MAD;		n_args[70]=1;		n_out[70]=1;
	ops[71]=grd_MAX;		n_args[71]=2;		n_out[71]=1;
	ops[72]=grd_MEAN;		n_args[72]=1;		n_out[72]=1;
	ops[73]=grd_MED;		n_args[73]=1;		n_out[73]=1;
	ops[74]=grd_MIN;		n_args[74]=2;		n_out[74]=1;
	ops[75]=grd_MODE;		n_args[75]=1;		n_out[75]=1;
	ops[76]=grd_MUL;		n_args[76]=2;		n_out[76]=1;
	ops[77]=grd_NAN;		n_args[77]=2;		n_out[77]=1;
	ops[78]=grd_NEG;		n_args[78]=1;		n_out[78]=1;
	ops[79]=grd_NEQ;		n_args[79]=2;		n_out[79]=1;
	ops[80]=grd_NRAND;		n_args[80]=2;		n_out[80]=1;
	ops[81]=grd_OR;		n_args[81]=2;		n_out[81]=1;
	ops[82]=grd_PDIST;		n_args[82]=1;		n_out[82]=1;
	ops[83]=grd_PLM;		n_args[83]=3;		n_out[83]=1;
	ops[84]=grd_POP;		n_args[84]=1;		n_out[84]=0;
	ops[85]=grd_POW;		n_args[85]=2;		n_out[85]=1;
	ops[86]=grd_R2;		n_args[86]=2;		n_out[86]=1;
	ops[87]=grd_R2D;		n_args[87]=1;		n_out[87]=1;
	ops[88]=grd_RAND;		n_args[88]=2;		n_out[88]=1;
	ops[89]=grd_RINT;		n_args[89]=1;		n_out[89]=1;
	ops[90]=grd_ROTX;		n_args[90]=2;		n_out[90]=1;
	ops[91]=grd_ROTY;		n_args[91]=2;		n_out[91]=1;
	ops[92]=grd_SAZ;		n_args[92]=2;		n_out[92]=1;
	ops[93]=grd_SDIST;		n_args[93]=2;		n_out[93]=1;
	ops[94]=grd_SIGN;		n_args[94]=1;		n_out[94]=1;
	ops[95]=grd_SIN;		n_args[95]=1;		n_out[95]=1;
	ops[96]=grd_SINC;		n_args[96]=1;		n_out[96]=1;
	ops[97]=grd_SIND;		n_args[97]=1;		n_out[97]=1;
	ops[98]=grd_SINH;		n_args[98]=1;		n_out[98]=1;
	ops[99]=grd_SQRT;		n_args[99]=1;		n_out[99]=1;
	ops[100]=grd_STD;		n_args[100]=1;		n_out[100]=1;
	ops[101]=grd_STEP;		n_args[101]=1;		n_out[101]=1;
	ops[102]=grd_STEPX;		n_args[102]=1;		n_out[102]=1;
	ops[103]=grd_STEPY;		n_args[103]=1;		n_out[103]=1;
	ops[104]=grd_SUB;		n_args[104]=2;		n_out[104]=1;
	ops[105]=grd_TAN;		n_args[105]=1;		n_out[105]=1;
	ops[106]=grd_TAND;		n_args[106]=1;		n_out[106]=1;
	ops[107]=grd_TANH;		n_args[107]=1;		n_out[107]=1;
	ops[108]=grd_TCRIT;		n_args[108]=2;		n_out[108]=1;
	ops[109]=grd_TDIST;		n_args[109]=2;		n_out[109]=1;
	ops[110]=grd_UPPER;		n_args[110]=1;		n_out[110]=1;
	ops[111]=grd_XOR;		n_args[111]=2;		n_out[111]=1;
	ops[112]=grd_Y0;		n_args[112]=1;		n_out[112]=1;
	ops[113]=grd_Y1;		n_args[113]=1;		n_out[113]=1;
	ops[114]=grd_YLM;		n_args[114]=2;		n_out[114]=2;
	ops[115]=grd_YN;		n_args[115]=2;		n_out[115]=1;
	ops[116]=grd_ZCRIT;		n_args[116]=1;		n_out[116]=1;
}
