/*--------------------------------------------------------------------
 *	$Id: grdmath_main.c,v 1.88 2006/05/08 23:31:05 pwessel Exp $
 *
 *	Copyright (c) 1991-2006 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
 *		27-JUL-2005 Pw: Added TN (Chebyshev)
 *		10-AUG-2005 PW: Added lower case x and y for normalized -1/+1 coordinates
 *		07-SEP-2005 PW: Added CORRCOEFF
 *		12-JAN-2006 PW: Use size_t for counters to allow > 2 Gb arrays
 *		      PW 03/22/06 Added CPOISS
 *		      PW 03/24/06 Added ZDIST
 * Version:	4.1
 *
 */
 
#include "gmt.h"

struct GRDMATH_CTRL {	/* All control options for this program (except common args) */
	/* active is TRUE if the option has been activated */
	struct F {	/* -F */
		BOOLEAN active;
	} F;
	struct I {	/* -Idx[/dy] */
		BOOLEAN active;
		double xinc, yinc;
	} I;
};

#define GRDMATH_N_OPERATORS	121
#define GRDMATH_STACK_SIZE	100

#define GRDMATH_ARG_IS_OPERATOR		 0
#define GRDMATH_ARG_IS_FILE		-1
#define GRDMATH_ARG_IS_NUMBER		-2
#define GRDMATH_ARG_IS_PI		-3
#define GRDMATH_ARG_IS_E		-4
#define GRDMATH_ARG_IS_X_MATRIX		-5
#define GRDMATH_ARG_IS_x_MATRIX		-6
#define GRDMATH_ARG_IS_Y_MATRIX		-7
#define GRDMATH_ARG_IS_y_MATRIX		-8
#define GRDMATH_ARG_IS_ASCIIFILE	-9
#define GRDMATH_ARG_IS_SAVE		-10
#define GRDMATH_N_SKIP_ARGS		 6

struct GRDMATH_INFO {
	size_t nm;
	char ASCII_file[BUFSIZ];
	float *grd_x, *grd_y;
	float *grd_xn, *grd_yn;
	struct GRD_HEADER header;
};

#include "grdmath_def.h"

int main (int argc, char **argv)
{
	int i, j, k, arg, op = 0, nstack = 0, new_stack = -1, ok = 1, n_items = 0, this_stack;
	int consumed_operands[GRDMATH_N_OPERATORS], produced_operands[GRDMATH_N_OPERATORS];

	size_t i64;

	BOOLEAN constant[GRDMATH_STACK_SIZE], error = FALSE, skip_arg, set_r = FALSE;

	float *stack[GRDMATH_STACK_SIZE];

	double factor[GRDMATH_STACK_SIZE], x_noise, y_noise, off, scale;

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

	struct GRD_HEADER grd[GRDMATH_STACK_SIZE];
	struct GMT_HASH *p, *current, localhashnode[GMT_HASH_SIZE];
	struct GRDMATH_INFO info;
	struct GRDMATH_CTRL *Ctrl;

	PFV call_operator[GRDMATH_N_OPERATORS];

	int decode_argument(char *txt, struct GMT_HASH *H);
	void grdmath_init (PFV ops[], int n_args[], int n_out[]);
	void *New_Grdmath_Ctrl (), Free_Grdmath_Ctrl (struct GRDMATH_CTRL *C);
	
	argc = GMT_begin (argc, argv);

	Ctrl = (struct GRDMATH_CTRL *) New_Grdmath_Ctrl ();		/* Allocate and initialize defaults in a new control structure */

	memset ((void *)&info, 0, sizeof (struct GRDMATH_INFO));

	if (argc == 2 && !strcmp (argv[1], "-")) error = GMT_give_synopsis_and_exit = TRUE;

	if (argc == 1 || GMT_give_synopsis_and_exit) {
		fprintf (stderr, "grdmath %s - Reverse Polish Notation (RPN) calculator for grdfiles (element by element)\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdmath [%s] [-F] [%s] [-V] %s\n", GMT_Rgeo_OPT, GMT_I_OPT, GMT_f_OPT);
		fprintf (stderr, "	A B op C op D op ... = outfile\n\n");

		if (GMT_give_synopsis_and_exit) 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", GRDMATH_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, "\t  Xn	= grid with normalized [-1|+1] x-coordinates\n");
		fprintf (stderr, "\t  Yn	= grid with normalized [-1|+1] 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]. Requires -R and -I.\n");
		GMT_inc_syntax ('I', 0);
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		GMT_explain_option ('f');
		exit (EXIT_FAILURE);
	}

	for (i = 1, error = TRUE; error && i < argc; i++) if (argv[i][0] == '=' && argv[i][1] == '\0') error = FALSE;
	if (error) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Usage is <operations> = outfile\n", GMT_program);
		exit (EXIT_FAILURE);
	}

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

	GMT_hash_init (localhashnode, operator, GMT_HASH_SIZE, GRDMATH_N_OPERATORS);

	for (i = 0; i < GRDMATH_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 (&info.header, argc, argv, FALSE);

	for (arg = 1; info.nm == 0 && arg < argc; arg++) {

		if (argv[arg][0] == '-') continue;	/* Command line option or a negative number */
		if (decode_argument (argv[arg], localhashnode) == GRDMATH_ARG_IS_SAVE) {	/* Output case = */
			arg++;	/* Skip the output file name since the file does not exist yet */
			while (argv[arg][0] == '-' && arg < argc) arg++;	/* In case someone placed an option bewteen = and output gridname */
			continue;
		}
		if (decode_argument (argv[arg], localhashnode) != GRDMATH_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 */

		if (GMT_read_grd_info (argv[arg], &info.header)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, info.header.name);
			exit (EXIT_FAILURE);
		}

		info.nm = ((size_t)info.header.nx) * ((size_t)info.header.ny);
	}

	/* Scan command line for -R, -I, -F, -f, -b, -V */

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

			switch (argv[arg][1]) {

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

				case 'F':
					Ctrl->F.active = TRUE;
					break;

				case 'I':
					Ctrl->I.active = TRUE;
					if (GMT_getinc (&argv[arg][2], &Ctrl->I.xinc, &Ctrl->I.yinc)) {
						GMT_inc_syntax ('I', 1);
						error = TRUE;
					}
					break;
			}
		}
	}
	if (info.nm && set_r && Ctrl->I.active) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot use -R, -I when grdfiles are specified\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (info.nm && Ctrl->F.active) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot use -F when grdfiles are specified\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if ((set_r + Ctrl->I.active)%2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  -R and -I must both be specified\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (Ctrl->I.active && (Ctrl->I.xinc <= 0.0 || Ctrl->I.yinc <= 0.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (set_r && Ctrl->I.active) {
		info.header.node_offset = Ctrl->F.active;
		info.header.x_inc = Ctrl->I.xinc;
		info.header.y_inc = Ctrl->I.yinc;
		GMT_RI_prepare (&info.header);	/* Ensure -R -I consistency and set nx, ny */
		GMT_grd_RI_verify (&info.header, 1);
		info.nm = ((size_t)info.header.nx) * ((size_t)info.header.ny);
		info.header.xy_off = 0.5 * info.header.node_offset;
	}

	if (info.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, info.nm, sizeof (float), GMT_program);

	/* Get x and y vectors */

	info.grd_x = (float *) GMT_memory (VNULL, (size_t)info.header.nx, sizeof (float), GMT_program);
	info.grd_y = (float *) GMT_memory (VNULL, (size_t)info.header.ny, sizeof (float), GMT_program);
	info.grd_xn = (float *) GMT_memory (VNULL, (size_t)info.header.nx, sizeof (float), GMT_program);
	info.grd_yn = (float *) GMT_memory (VNULL, (size_t)info.header.ny, sizeof (float), GMT_program);
	for (j = 0; j < info.header.ny; j++) info.grd_y[j] = (float)GMT_j_to_y (j, info.header.y_min, info.header.y_max, info.header.y_inc, info.header.xy_off, info.header.ny);
	for (i = 0; i < info.header.nx; i++) info.grd_x[i] = (float)GMT_i_to_x (i, info.header.x_min, info.header.x_max, info.header.x_inc, info.header.xy_off, info.header.nx);
	off = 0.5 * (info.header.x_max + info.header.x_min);
	scale = 2.0 / (info.header.x_max - info.header.x_min);
	for (i = 0; i < info.header.nx; i++) info.grd_xn[i] = (float)((info.grd_x[i] - off) * scale);
	off = 0.5 * (info.header.y_max + info.header.y_min);
	scale = 2.0 / (info.header.y_max - info.header.y_min);
	for (j = 0; j < info.header.ny; j++) info.grd_yn[j] = (float)((info.grd_y[j] - off) * scale);
	x_noise = GMT_SMALL * info.header.x_inc;	y_noise = GMT_SMALL * info.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 < GRDMATH_N_SKIP_ARGS; k++) skip_arg = !strncmp (argv[arg], skip_args[k], 2);
		if (skip_arg) continue;

		op = decode_argument (argv[arg], localhashnode);

		if (op != GRDMATH_ARG_IS_FILE && !GMT_access(argv[arg], R_OK)) fprintf (stderr, "%s Warning: The number or operator %s may be confused with an existing file %s!\n", GMT_program, argv[arg], argv[arg]);

		if (op == GRDMATH_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 (&info.header, argc, argv, TRUE);

			if (n_items && new_stack < 0 && constant[nstack-1]) {	/* Only a constant provided, set grid accordingly */
				for (i64 = 0; i64 < info.nm; i64++) stack[nstack-1][i64] = (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, &info.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 (op < GRDMATH_ARG_IS_OPERATOR) {	/* File name or factor */

			if (op == GRDMATH_ARG_IS_FILE && !(strncmp (argv[arg+1], "LDIST", 5) && strncmp (argv[arg+1], "PDIST", 5) && strncmp (argv[arg+1], "INSIDE", 6))) op = GRDMATH_ARG_IS_ASCIIFILE;

			if (nstack == GRDMATH_STACK_SIZE) {	/* Stack overflow */
				error = TRUE;
				continue;
			}
			n_items++;
			if (op == GRDMATH_ARG_IS_NUMBER) {
				constant[nstack] = TRUE;
				ok = sscanf (argv[arg], "%lf", &factor[nstack]);
				error = !ok;
				if (gmtdefs.verbose) fprintf (stderr, "%g ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (op == GRDMATH_ARG_IS_PI) {
				constant[nstack] = TRUE;
				factor[nstack] = M_PI;
				if (gmtdefs.verbose) fprintf (stderr, "%g ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (op == GRDMATH_ARG_IS_E) {
				constant[nstack] = TRUE;
				factor[nstack] = M_E;
				if (gmtdefs.verbose) fprintf (stderr, "%g ", 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, info.nm, sizeof (float), GMT_program);
			constant[nstack] = FALSE;

			if (op == GRDMATH_ARG_IS_X_MATRIX) {		/* Need to set up matrix of x-values */
				if (gmtdefs.verbose) fprintf (stderr, "X ");
				for (j = k = 0; j < info.header.ny; j++, k += info.header.nx)
					memcpy ((void *)&stack[nstack][k], (void *)info.grd_x, (size_t)(info.header.nx * sizeof (float)));
			}
			else if (op == GRDMATH_ARG_IS_x_MATRIX) {		/* Need to set up matrix of normalized x-values */
				if (gmtdefs.verbose) fprintf (stderr, "Xn ");
				for (j = k = 0; j < info.header.ny; j++, k += info.header.nx)
					memcpy ((void *)&stack[nstack][k], (void *)info.grd_xn, (size_t)(info.header.nx * sizeof (float)));
			}
			else if (op == GRDMATH_ARG_IS_Y_MATRIX) {	/* Need to set up matrix of y-values */
				if (gmtdefs.verbose) fprintf (stderr, "Y ");
				for (j = k = 0; j < info.header.ny; j++) for (i = 0; i < info.header.nx; i++, k++)
					stack[nstack][k] = info.grd_y[j];
			}
			else if (op == GRDMATH_ARG_IS_y_MATRIX) {	/* Need to set up matrix of normalized y-values */
				if (gmtdefs.verbose) fprintf (stderr, "Yn ");
				for (j = k = 0; j < info.header.ny; j++) for (i = 0; i < info.header.nx; i++, k++)
					stack[nstack][k] = info.grd_yn[j];
			}
			else if (op == GRDMATH_ARG_IS_ASCIIFILE) {
				strcpy (info.ASCII_file, argv[arg]);
				if (gmtdefs.verbose) fprintf (stderr, "(%s) ", argv[arg]);
			}
			else if (op == GRDMATH_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 != info.header.nx || grd[nstack].ny != info.header.ny) {
					fprintf (stderr, "%s: grd files not of same size!\n", GMT_program);
					exit (EXIT_FAILURE);
				}
				else if (fabs (grd[nstack].x_min - info.header.x_min) > x_noise || fabs (grd[nstack].x_max - info.header.x_max) > x_noise || 
					fabs (grd[nstack].y_min - info.header.y_min) > y_noise || fabs (grd[nstack].y_max - info.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 ((new_stack = nstack - consumed_operands[op] + produced_operands[op]) >= GRDMATH_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, info.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, info.nm, sizeof (float), GMT_program);
		}

		(*call_operator[op]) (&info, 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 < GRDMATH_STACK_SIZE; i++) if (stack[i]) GMT_free ((void *)stack[i]);
	GMT_free ((void *)info.grd_x);
	GMT_free ((void *)info.grd_y);
	GMT_free ((void *)info.grd_xn);
	GMT_free ((void *)info.grd_yn);
	for (i = 0; i < GMT_HASH_SIZE; i++) {
		p = localhashnode[i].next;
		while ((current = p)) {
			p = p->next;
			GMT_free ((void *)current);
		}
	}

	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);

	Free_Grdmath_Ctrl (Ctrl);	/* Deallocate control structure */

	GMT_end (argc, argv);

	exit (EXIT_SUCCESS);
}

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

void grd_ABS (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : fabs ((double)stack[last][i]));
}

void grd_ACOS (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_acos ((double)stack[last][i]));
}

void grd_ACOSH (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : acosh ((double)stack[last][i]));
}

void grd_ADD (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_asin ((double)stack[last][i]));
}

void grd_ASINH (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_ATAN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_ATAN2 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : atanh ((double)stack[last][i]));
}

void grd_BEI (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_BER (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_CAZ (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev;
	size_t k;
	double x, y, az;

	prev = last - 1;
	for (k = j = i = 0; k < info->nm; k++) {
		x = (constant[prev]) ? factor[prev] : stack[prev][k];
		y = (constant[last]) ? factor[last] : stack[last][k];
		az = 90.0 - (atan2 (y - (double)info->grd_y[j], x - (double)info->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 == info->header.nx) i = 0, j++;
	}
}

void grd_CDIST (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev;
	size_t k;
	double a, b;

	prev = last - 1;
	for (k = j = i = 0; k < info->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)info->grd_x[i], b - (double)info->grd_y[j]);
		i++;
		if (i == info->header.nx) i = 0, j++;
	}
}

void grd_CEIL (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_CHICRIT (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	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 < info->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_CORRCOEFF (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	double coeff;

	prev = last - 1;
	if (constant[prev] || constant[last]) {
		fprintf (stderr, "%s: Error, cannot have constant operands for CORRCOEFF!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	coeff = GMT_corrcoeff_f (stack[prev], stack[last], info->nm, 0);
	for (i = 0; i < info->nm; i++) stack[prev][i] = (float)coeff;
}

void grd_COS (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_COSD (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_COSH (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_CPOISS (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	double a, b, prob;

	prev = last - 1;
	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 < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		GMT_cumpoisson (a, b, &prob);
		stack[prev][i] = (float)prob;
	}
}

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

	/* Curvature (Laplacian) */

	easy = (info->header.x_inc == info->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, info->nm * sizeof (float));
		return;
	}

	z = (float *) (float *) GMT_memory (VNULL, info->nm, sizeof (float), GMT_program);
 
	nx = (size_t)info->header.nx;
	cx = +0.5 / (info->header.x_inc * info->header.x_inc);
	cy = -0.5 / (info->header.y_inc * info->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 < info->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 < info->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 < info->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 = info->nm - nx + 1; i < info->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 < info->header.ny-1; j++) {
		k++;
		for (i = 1; i < info->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, info->nm * sizeof (float));
	GMT_free ((void *)z);
}

void grd_D2DX2 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j;
	size_t 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, info->nm * sizeof (float));
		return;
	}

	c = 1.0 / (info->header.x_inc * info->header.x_inc);
	for (j = k = 0; j < info->header.ny; j++) {
		next_left = stack[last][k];
		stack[last][k++] = (float)0.0;
		for (i = 1; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j;
	size_t 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, info->nm * sizeof (float));
		return;
	}

	c = 1.0 / (info->header.y_inc * info->header.y_inc);
	nx = (size_t)info->header.nx;
	for (i = 0; i < info->header.nx; i++) {
		k = i;
		next_bottom = stack[last][k];
		stack[last][k] = (float)0.0;	/* Natural BC has curvature = 0 at edges */
		k += nx;
		for (j = 1; j < info->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)0.0;	/* Natural BC has curvature = 0 at edges */
	}
}

void grd_D2R (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = factor[last] * D2R;
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : (stack[last][i] * D2R));
}

void grd_DDX (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j;
	size_t 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, info->nm * sizeof (float));
		return;
	}

	c = 0.5 / info->header.x_inc;
	for (j = k = 0; j < info->header.ny; j++) {
		next_left = 2.0 * stack[last][k] - stack[last][k+1];
		for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j;
	size_t 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, info->nm * sizeof (float));
		return;
	}

	c = -0.5 / info->header.y_inc;	/* Because the loop over j below goes from ymax to ymin we compensate with a minus sign here */
	nx = (size_t)info->header.nx;
	for (i = 0; i < info->header.nx; i++) {
		k = i;
		next_bottom = 2.0 * stack[last][k] - stack[last][k+nx];
		for (j = 0; j < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = GMT_dilog (factor[last]);
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_dilog (stack[last][i]));
}

void grd_DIV (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	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 (info, stack, constant, factor, last);
		factor[last] = a;	/* Restore factor */
		return;
	}

	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int next;
	size_t i;

	next = last + 1;
	constant[next] = constant[last];
	factor[next] = factor[last];
	if (constant[last]) {	/* Time to fess up */
		for (i = 0; i < info->nm; i++) stack[last][i] = (float)factor[last];
	}

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

void grd_ERF (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_ERFC (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_EQ (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_f_NaN : (float)(a == b);
	}
}

void grd_ERFINV (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_EXCH (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev;
	size_t i;

	prev = last - 1;
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : exp ((double)stack[last][i]));
}

void grd_EXTREMA (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, nx1, ny1, dx, dy, diag;
	size_t k;
	int do_derivative (float *z, size_t this_node, int this_pos, int last_i, int last_j, int type);
	int do_diagonal (float *z, size_t 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, info->nm * sizeof (float));
		return;
	}

	z = (float *) GMT_memory (VNULL, info->nm, sizeof (float), GMT_program);
 
	nx1 = info->header.nx - 1;
	ny1 = info->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 < info->header.ny; j++) {	/* k is middle (current) node */
		for (i = 0; i < info->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, info->header.nx, 0)) == -2) continue;	/* Too many NaNs or flat x-line */
			if ((dy = do_derivative (stack[last], k, j, ny1, info->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, info->nm * sizeof (float));
	GMT_free ((void *)z);
}

/* Subroutines for grd_EXTREMA */

int do_derivative (float *z, size_t 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;
	size_t 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, size_t 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 */

	size_t next_node, prev_node;
	int off;

	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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t k, k0, nx1;
	int i, j, 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 = info->header.nx / 2;
	nx1 = (size_t)info->header.nx - 1;
	for (j = k0 = 0; j < info->header.ny; j++, k0 += (size_t)info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t  k, ny1, nx;
	int i, j, 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 = info->header.ny / 2;
	ny1 = (size_t)info->header.ny - 1;
	nx = (size_t)info->header.nx;
	for (i = 0; i < info->header.nx; i++) {
		for (j = 0, k = ny1; j < ny_half; j++, k--) {	/* Do this to all rows */
			f_swap (stack[last][j*nx+i], stack[last][k*nx+i]);
		}
	}
}

void grd_FLOOR (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_FMOD (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_f_NaN : (float)(a >= b);
	}
}

void grd_GT (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_f_NaN : (float)(a > b);
	}
}

void grd_HYPOT (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_I1 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_IN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev, order = 0;
	BOOLEAN simple = FALSE;
	double b = 0.0;

	prev = last - 1;
	if (constant[last]) {
		if (gmtdefs.verbose && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for IN!\n", GMT_program);
		if (gmtdefs.verbose && fabs (rint(factor[last]) - factor[last]) > GMT_SMALL) fprintf (stderr, "%s: Warning, order not an integer for IN!\n", GMT_program);
		order = irint (fabs (factor[last]));
		if (constant[prev]) {
			b = GMT_in (order, fabs (factor[prev]));
			simple = TRUE;
		}
	}
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev1, prev2;
	BOOLEAN inrange;
	float a = 0.0, b = 0.0, c = 0.0;

	/* 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	void grd_INSIDE_GEO (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last);
	void grd_INSIDE_XY (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last);

	if (GMT_io.in_col_type[0] & GMT_IS_GEO)	/* Geographic data */
		grd_INSIDE_GEO (info, stack, constant, factor, last);
	else					/* Cartesian data */
		grd_INSIDE_XY (info, stack, constant, factor, last);
}

void grd_INSIDE_GEO (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{	/* Suitable for geographic (lon, lat) data and polygons */
	int i, j, P, inside;
	size_t k;
	BOOLEAN greenwich;
	struct GMT_TABLE *polygon;

	greenwich = (info->header.x_min < 0.0 && info->header.x_max > 0.0);

	GMT_import_table ((void *)(info->ASCII_file), GMT_IS_FILE, &polygon, 0.0, greenwich, TRUE, TRUE);	/* Closed polygons */

	for (i = j = k = 0; k < info->nm; k++) {	/* Visit each node */
		for (P = inside = 0; !inside && P < polygon->n_segments; P++) {
			inside = GMT_inonout_sphpol (info->grd_x[i], info->grd_y[j], polygon->segment[P]);
		}
		stack[last][k] = (float)((inside) ? 1.0 : 0.0);
		i++;
		if (i == info->header.nx) {	/* Go to next row */
			i = 0, j++;
		}
	}

	/* Free memory used for pol */

	GMT_free_table (polygon);
}

void grd_INSIDE_XY (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{	/* Version for Cartesian data */
	int i, j, P, inside;
	size_t k;
	struct GMT_TABLE *polygon;

	GMT_import_table ((void *)(info->ASCII_file), GMT_IS_FILE, &polygon, 0.0, FALSE, TRUE, TRUE);	/* Closed polygons */

	for (i = j = k = 0; k < info->nm; k++) {	/* Visit each node */
		for (P = inside = 0; !inside && P < polygon->n_segments; P++) {
			if (info->grd_y[j] < polygon->segment[P]->min[GMT_Y] || info->grd_y[j] > polygon->segment[P]->max[GMT_Y]) continue;	/* Outside y range */
			if (info->grd_x[j] < polygon->segment[P]->min[GMT_X] || info->grd_x[j] > polygon->segment[P]->max[GMT_X]) continue;	/* Outside x range */
			inside = GMT_non_zero_winding (info->grd_x[i], info->grd_y[j], polygon->segment[P]->coord[GMT_X], polygon->segment[P]->coord[GMT_Y], polygon->segment[P]->n_rows);	/* Must test */
		}
		stack[last][k] = (float)((inside) ? 1.0 : 0.0);
		i++;
		if (i == info->header.nx) {	/* Go to next row */
			i = 0, j++;
		}
	}

	/* Free memory used for pol */

	GMT_free_table (polygon);
}

void grd_INV (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t 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 < info->nm; i++) {
		a = (constant[last]) ? factor[last] : 1.0 / stack[last][i];
		stack[last][i] = (float)a;
	}
}

void grd_ISNAN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = GMT_is_dnan (factor[last]);
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_is_fnan (stack[last][i]));
}

void grd_J0 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_J1 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_JN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev, order = 0;
	BOOLEAN simple = FALSE;
	double b = 0.0;

	prev = last - 1;
	if (constant[last]) {
		if (gmtdefs.verbose && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for JN!\n", GMT_program);
		if (gmtdefs.verbose && fabs (rint(factor[last]) - factor[last]) > GMT_SMALL) fprintf (stderr, "%s: Warning, order not an integer for JN!\n", GMT_program);
		order = irint (fabs (factor[last]));
		if (constant[prev]) {
			b = jn (order, fabs (factor[prev]));
			simple = TRUE;
		}
	}
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_K1 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_KEI (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_KER (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_KN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev, order = 0;
	BOOLEAN simple = FALSE;
	double b = 0.0;

	prev = last - 1;
	if (constant[last]) {
		if (gmtdefs.verbose && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for KN!\n", GMT_program);
		if (gmtdefs.verbose && fabs (rint(factor[last]) - factor[last]) > GMT_SMALL) fprintf (stderr, "%s: Warning, order not an integer for KN!\n", GMT_program);
		order = irint (fabs (factor[last]));
		if (constant[prev]) {
			b = GMT_kn (order, fabs (factor[prev]));
			simple = TRUE;
		}
	}
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j;
	size_t k;
	double d;
	BOOLEAN greenwich;
	PFI near_a_line;
	struct GMT_TABLE *line;

	if (GMT_io.in_col_type[0] & GMT_IS_GEO) {	/* Geographic data */
		GMT_distance_func = (PFD) ((GMT_IS_SPHERICAL) ? GMT_great_circle_dist_km : GMT_geodesic_dist_km);
		near_a_line = (PFI) GMT_near_a_line_spherical;
		greenwich = (info->header.x_min < 0.0 && info->header.x_max > 0.0);
	}
	else {
		GMT_distance_func = (PFD) GMT_cartesian_dist;
		near_a_line = (PFI) GMT_near_a_line_cartesian;
		greenwich = FALSE;
	}

	GMT_import_table ((void *)(info->ASCII_file), GMT_IS_FILE, &line, 0.0, greenwich, FALSE, TRUE);

	for (i = j = k = 0; k < info->nm; k++) {	/* Visit each node */
		(void) near_a_line (info->grd_x[i], info->grd_y[j], line, TRUE, &d, NULL, NULL);
		stack[last][k] = (float)d;
		i++;
		if (i == info->header.nx) {
			if (gmtdefs.verbose == 2) fprintf (stderr, "Line %d\n", j);
			i = 0, j++;
		}
	}

	/* Free memory used for line */

	GMT_free_table (line);
}

void grd_LE (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_f_NaN : (float)(a <= b);
	}
}

void grd_LOG (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_log (fabs ((double)stack[last][i])));
}

void grd_LOG10 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_log10 (fabs ((double)stack[last][i])));
}

void grd_LOG1P (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : d_log1p (fabs ((double)stack[last][i])));
}

void grd_LOG2 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)(((constant[last]) ? a : d_log (fabs ((double)stack[last][i]))) * M_LN2_INV);
}

void grd_LMSSCL (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int GMT_mode_selection = 0, GMT_n_multiples = 0;
	size_t i;
	double mode, lmsscl;
	float lmsscl_f;

	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last], 0, info->nm * sizeof (float));
		return;
	}

	/* Sort will put any NaNs to the end - we then count to find the real data */

	qsort ((void *)stack[last], info->nm, sizeof (float), GMT_comp_float_asc);
	for (i = info->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 < info->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);
}

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

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < info->nm; i++) stack[last][i] = (float)factor[last];
		return;
	}

	for (i = 0, low = FLT_MAX; i < info->nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		if (stack[last][i] < low) low = stack[last][i];
	}
	for (i = 0; i < info->nm; i++) if (!GMT_is_fnan (stack[last][i])) stack[last][i] = low;
}

void grd_LRAND (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a = 0.0, b = 0.0;

	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		b = (constant[last]) ? factor[last] : stack[last][i];
		stack[prev][i] = (GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_f_NaN : (float)(a < b);
	}
}

void grd_MAD (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double mad, med;
	float mad_f;

	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last], 0, info->nm * sizeof (float));
		return;
	}

	/* Sort will put any NaNs to the end - we then count to find the real data */

	qsort ((void *)stack[last], info->nm, sizeof (float), GMT_comp_float_asc);
	for (i = info->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 < info->nm; i++) stack[last][i] = mad_f;
}

void grd_MAX (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i, n_a = 0;
	double sum_a = 0.0;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < info->nm; i++) stack[last][i] = (float)factor[last];
		return;
	}

	for (i = 0; i < info->nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		sum_a += stack[last][i];
		n_a++;
	}
	sum_a = (n_a) ? sum_a / (double)n_a : 0.0;
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)sum_a;
}

void grd_MED (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	float med;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < info->nm; i++) stack[last][i] = (float)factor[last];
		return;
	}

	qsort ((void *)stack[last], info->nm, sizeof (float), GMT_comp_float_asc);
	for (i = info->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 < info->nm; i++) stack[last][i] = med;
}

void grd_MIN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int GMT_mode_selection = 0, GMT_n_multiples = 0;
	size_t i;
	double mode = 0.0;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < info->nm; i++) stack[last][i] = (float)factor[last];
		return;
	}

	qsort ((void *)stack[last], info->nm, sizeof (float), GMT_comp_float_asc);
	for (i = info->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 < info->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);
}

void grd_MUL (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a = 0.0, b = 0.0;

	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : -stack[last][i]);
}

void grd_NEQ (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a = 0.0, b = 0.0;

	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a, b;

	prev = last - 1;
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, dummy[2];
	size_t k;
	struct GMT_TABLE *T;

	if (GMT_io.in_col_type[0] & GMT_IS_GEO) {	/* Geographic data */
		GMT_distance_func = (PFD) ((GMT_IS_SPHERICAL) ? GMT_great_circle_dist_km : GMT_geodesic_dist_km);
	}
	else {
		GMT_distance_func = (PFD) GMT_cartesian_dist;
	}

	GMT_import_table ((void *)info->ASCII_file, GMT_IS_FILE, &T, 0.0, FALSE, FALSE, FALSE);

	for (i = j = k = 0; k < info->nm; k++) {	/* Visit each node */
		stack[last][k] = (float)GMT_dist_to_point (info->grd_x[i], info->grd_y[j], T, dummy);
		i++;
		if (i == info->header.nx) i = 0, j++;
	}

	/* Free memory used for points */

	GMT_free_table (T);
}

void grd_POP (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{

	/* Dummy routine that does nothing but consume the top element of stack */
}

void grd_PLM (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int prev, first, L, M;
	size_t i;
	double a = 0.0;
				/* 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 < info->nm; i++) stack[first][i] = (float)((constant[first]) ? a : GMT_plm (L, M, stack[first][i]));
}

void grd_POW (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = R2D * factor[last];
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : R2D * stack[last][i]);
}

void grd_RAND (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a = 0.0, b = 0.0;

	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_ROTX (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev, shift, *node;
	size_t k, nx;
	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] / info->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 += info->header.nx;	/* Same thing */
	nx = (size_t) info->header.nx;
	/* Set up permutation vector */

	node = (int *) GMT_memory (VNULL, nx, sizeof (int), GMT_program);
	z =  (float *) GMT_memory (VNULL, nx, sizeof (float), GMT_program);
	for (i = 0; i < info->header.nx; i++) node[i] = (i + shift) % info->header.nx;	/* Move by shift but rotate around */
	for (j = k = 0; j < info->header.ny; j++, k += nx) {	/* For each row */
		for (i = 0; i < info->header.nx; i++) z[node[i]] = stack[prev][k+i];
		memcpy ((void *)&stack[prev][k], (void *)z, nx * sizeof (float));
	}
	GMT_free ((void *)z);
	GMT_free ((void *)node);
}

void grd_ROTY (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev, shift, *node;
	size_t nx;
	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] / info->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 += info->header.ny;	/* Same thing */
	nx = (size_t) info->header.nx;
	/* Set up permutation vector */

	node = (int *) GMT_memory (VNULL, (size_t)info->header.ny, sizeof (int), GMT_program);
	z =  (float *) GMT_memory (VNULL, (size_t)info->header.ny, sizeof (float), GMT_program);
	for (j = 0; j < info->header.ny; j++) node[j] = (j + info->header.ny - shift) % info->header.ny;	/* Move by shift but rotate around */
	for (i = 0; i < info->header.nx; i++) {	/* For each column */
		for (j = 0; j < info->header.ny; j++) z[node[j]] = stack[prev][j*nx+i];
		for (j = 0; j < info->header.ny; j++) stack[prev][j*nx+i] = z[j];
	}
	GMT_free ((void *)z);
	GMT_free ((void *)node);
}

void grd_SDIST (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev;
	size_t k;
	double a, b;

	GMT_distance_func = (PFD) ((GMT_IS_SPHERICAL) ? GMT_geodesic_dist_degree : GMT_great_circle_dist);
	prev = last - 1;
	for (k = j = i = 0; k < info->nm; k++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][k];
		b = (constant[last]) ? factor[last] : stack[last][k];
		stack[prev][k] = (float)(GMT_distance_func) (a, b, (double)info->grd_x[i], (double)info->grd_y[j]);
		i++;
		if (i == info->header.nx) i = 0, j++;
	}
}

void grd_SAZ (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev;
	size_t k;
	double x0 = 0.0, y0 = 0.0, 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)info->header.ny, sizeof (double), GMT_program);
	sy = (double *) GMT_memory (VNULL, (size_t)info->header.ny, sizeof (double), GMT_program);
	for (j = 0; j < info->header.ny; j++) sincos (info->grd_y[j] * D2R, &sy[j], &cy[j]);
	for (k = j = i = 0; k < info->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 - info->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 == info->header.nx) i = 0, j++;
	}
	GMT_free ((void *)cy);
	GMT_free ((void *)sy);
}

void grd_SIGN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : copysign (1.0, (double)stack[last][i]));
}

void grd_SIN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_SINC (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_SIND (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_SINH (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_SQRT (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	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 < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : sqrt ((double)stack[last][i]));
}

void grd_STD (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i, n_a = 0;
	double sum_a = 0.0, sum_a2 = 0.0, std;

	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last], 0, info->nm * sizeof (float));
		return;
	}

	for (i = 0; i < info->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 < info->nm; i++) stack[last][i] = (float)std;
}

void grd_STEP (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = factor[last];
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a;

	for (i = 0; i < info->nm; i++) {
		a = info->grd_x[i%info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a;

	for (i = 0; i < info->nm; i++) {
		a = info->grd_y[i/info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_TAND (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_TANH (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_TN (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev, n;
	double a = 0.0;

	prev = last - 1;
	for (i = 0; i < info->nm; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][i];
		n = irint ((double)((constant[last]) ? factor[last] : stack[last][i]));
		stack[prev][i] = (float)GMT_chebyshev (a, n);
	}
}

void grd_TCRIT (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int 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 < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	float high;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < info->nm; i++) stack[last][i] = (float)factor[last];
		return;
	}

	for (i = 0, high = -FLT_MAX; i < info->nm; i++) {
		if (GMT_is_fnan (stack[last][i])) continue;
		if (stack[last][i] > high) high = stack[last][i];
	}
	for (i = 0; i < info->nm; i++) if (!GMT_is_fnan (stack[last][i])) stack[last][i] = high;
}

void grd_XOR (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev;
	double a = 0.0, b = 0.0;

	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_Y1 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

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

void grd_YLM (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, prev, L, M;
	size_t k;
	double x, z, normal, P, *C = VNULL, *S = VNULL;

	/* 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)info->header.nx, sizeof (double), GMT_program);
		S = (double *) GMT_memory (VNULL, (size_t)info->header.nx, sizeof (double), GMT_program);
		z = M * D2R;
		for (i = 0; i < info->header.nx; i++) {
			x = z * info->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 < info->header.ny; j++) {	/* For each latitude */

		x = cosd (90.0 - info->grd_y[j]);	/* plm takes cos(colatitude) */
		P = normal * GMT_plm (L, M, x);
		if (M > 0) {
			for (i = 0; i < info->header.nx; i++, k++) {
				stack[prev][k] = (float)(P * C[i]);
				stack[last][k] = (float)(P * S[i]);
			}
		}
		else {
			for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	int prev, order = 0;
	double b = 0.0;
	BOOLEAN simple = FALSE;

	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, argument = 0 for YN!\n", GMT_program);
	if (constant[last]) {
		if (gmtdefs.verbose && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for YN!\n", GMT_program);
		if (gmtdefs.verbose && (rint(factor[last]) != factor[last])) fprintf (stderr, "%s: Warning, order not an integer for YN!\n", GMT_program);
		order = irint (fabs (factor[last]));
		if (constant[prev]) {
			b = yn (order, fabs (factor[prev]));
			simple = TRUE;
		}
	}
	for (i = 0; i < info->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 (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = GMT_zcrit (factor[last]);
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_zcrit (a));
}

void grd_ZDIST (struct GRDMATH_INFO *info, float *stack[], BOOLEAN *constant, double *factor, int last)
{
	size_t i;
	double a = 0.0;

	if (constant[last]) a = GMT_zdist (factor[last]);
	for (i = 0; i < info->nm; i++) stack[last][i] = (float)((constant[last]) ? a : GMT_zdist (a));
}

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

int decode_argument (char *txt, struct GMT_HASH *H)
{
	int number, minus, plus, exp, dec, n_digits = 0, i;
	char *s, file[BUFSIZ];
	int get_operator (char *choice, struct GMT_HASH *H);

	if (!strcmp (txt, "=")) return GRDMATH_ARG_IS_SAVE;	/* Time to save stack; next arg is filename */

	/* Check if argument is operator */

	if ((i = get_operator (txt, H)) >= GRDMATH_ARG_IS_OPERATOR) return (i);

	/* Next look for symbols with special meaning */

	if (!(strcmp (txt, "PI") && strcmp (txt, "pi"))) return GRDMATH_ARG_IS_PI;
	if (!(strcmp (txt, "E") && strcmp (txt, "e"))) return GRDMATH_ARG_IS_E;
	if (!strcmp (txt, "X")) return GRDMATH_ARG_IS_X_MATRIX;
	if (!strcmp (txt, "Xn")) return GRDMATH_ARG_IS_x_MATRIX;
	if (!strcmp (txt, "Y")) return GRDMATH_ARG_IS_Y_MATRIX;
	if (!strcmp (txt, "Yn")) return GRDMATH_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);
	}
	if (number && n_digits > 0) return GRDMATH_ARG_IS_NUMBER;

	/* Determine if argument is file */

	strcpy (file, txt);
	for (i = 0; file[i]; i++) if (file[i] == '=' || file[i] == '?') file[i] = 0;
	if (!GMT_access (file, R_OK)) return GRDMATH_ARG_IS_FILE;

	fprintf (stderr, "%s: GMT SYNTAX ERROR: %s is not a number, operator or file name\n", GMT_program, txt);
	exit (EXIT_FAILURE);
}

int get_operator (char *choice, struct GMT_HASH *H) {
	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, H, GMT_HASH_SIZE, GMT_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 *New_Grdmath_Ctrl () {	/* Allocate and initialize a new control structure */
	struct GRDMATH_CTRL *C;
	
	C = (struct GRDMATH_CTRL *) GMT_memory (VNULL, 1, sizeof (struct GRDMATH_CTRL), "New_Grdmath_Ctrl");
	
	/* Initialize values whose defaults are not 0/FALSE/NULL */
			
	return ((void *)C);
}

void Free_Grdmath_Ctrl (struct GRDMATH_CTRL *C) {	/* Deallocate control structure */
	GMT_free ((void *)C);	
}

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