/*--------------------------------------------------------------------
 *	$Id: grdmath_main.c,v 1.5.4.4 2002/02/14 22:56:03 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * grdmath.c is a reverse polish calculator that operates on grd files 
 * (and constants) and perform basic matheatical 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-JUN-1995
 * Version:	3.0 based on old 2.1
 * Revised:	22-JUN-1998
 * Version:	3.1 based on old 3.0
 * Version:	3.3.2 based on old 3.1
 * Revised:	13-JUL-2000
 * Version:	3.3.6
 * 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
 * Version:	3.4
 *		3.4.1	EXTREMA did not consider edges
 *
 */
 
#include "gmt.h"

#define N_OPERATORS	99
#define STACK_SIZE	50

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

#include "grdmath_def.h"

PFV call_operator[N_OPERATORS];

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

struct GMT_HASH hashnode[HASH_SIZE];

float *grd_x, *grd_y;

struct GRD_HEADER header;

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

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

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

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

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

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

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

	GMT_hash_init (hashnode, operator, HASH_SIZE, N_OPERATORS);

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

	for (arg = 1; nm == 0 && arg < last_arg; arg++) {
	
		if (argv[arg][0] == '-') continue;	/* Command line option */
		if (decode_argument (argv[arg]) != ARG_IS_FILE) continue;
		
		strcpy (file, argv[arg]);
		for (j = 0; file[j]; j++) if (file[j] == '=') file[j] = 0;
		if (GMT_read_grd_info (argv[arg], &header)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file);
			exit (EXIT_FAILURE);
		}
		
		nm = header.nx * header.ny;
	}
	
	/* Scan command line for -R, -I, -F */

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

			switch (argv[arg][1]) {

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

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

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

	/* Get x and y vectors */

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

	grdmath_init (call_operator, consumed_operands, produced_operands);

	nstack = 0;
	
	for (arg = 1; !error && arg < last_arg; arg++) {
	
		/* First check if we should skip optional arguments */

		if (!(strncmp (argv[arg], "-R", 2) && strncmp (argv[arg], "-I", 2) && strncmp (argv[arg], "-F", 2) && strcmp (argv[arg], "-V"))) continue;

		if ((type = decode_argument (argv[arg])) != ARG_IS_OPERATOR) {	/* File name or factor */
		
			if (nstack == STACK_SIZE) {	/* Stack overflow */
				error = TRUE;
				continue;
			}

			if (type == ARG_IS_NUMBER) {
				constant[nstack] = TRUE;
				ok = sscanf (argv[arg], "%lf", &factor[nstack]);
				error = !ok;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (type == ARG_IS_PI) {
				constant[nstack] = TRUE;
				factor[nstack] = M_PI;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (type == ARG_IS_E) {
				constant[nstack] = TRUE;
				factor[nstack] = M_E;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}

			/* Here we need a matrix */

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

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

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

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

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

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

			/* Must make space for more */

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

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

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

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

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

	GMT_grd_init (&header, argc, argv, TRUE);
	
	if (new_stack < 0 && constant[0]) {	/* Only a constant provided, set grid accordingly */
		for (k = 0; k < nm; k++) stack[0][k] = (float)factor[0];
	}

	if (GMT_write_grd (outfile, &header, stack[0], 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);
	}
	
	for (i = 0; i < STACK_SIZE; i++) if (stack[i]) GMT_free ((void *)stack[i]);
	GMT_free ((void *)grd_x);
	GMT_free ((void *)grd_y);
	
	if (gmtdefs.verbose) fprintf (stderr, "\n");

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

	GMT_end (argc, argv);
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	/* Curvature (Laplacian) */

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

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to CURV is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	
	z = (float *) (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
 
	nx = header.nx;
	cx = 0.5 / (header.x_inc * header.x_inc);
	cy = 0.5 / (header.y_inc * header.y_inc);

	/* First left/right */

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

	/* Then top/bottom */

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

	/* Then inside */

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

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

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

	/* Central 2nd difference in x */

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

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

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

	/* Central 2nd difference in y */

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

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

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

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

	/* Central 1st difference in x */

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

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

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

	/* Central 1st difference in y */

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

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

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

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

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

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

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

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

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

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

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

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

void grd_EXTREMA (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, j, k, nx, nx2, ny1, nx1, x_status, y_status, xy_status[11];
	float *z, *z_left, *z_right, *z_top, *z_bot, left, right, top, bot;

	/* Find local extrema in grid */

	if (constant[last]) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand to EXTREMA is constant!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	if (header.nx < 2 || header.ny < 2) {	/* Cannot deal with profiles camouflaging as grids */
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning, operand EXTREMA requires both nx,ny to be >= 2!\n", GMT_program);
		memset ((void *)stack[last], 0, (size_t)(nm * sizeof (float)));
		return;
	}
	
	z = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
	z_left  = (float *) GMT_memory (VNULL, (size_t)header.ny, sizeof (float), GMT_program);
	z_right = (float *) GMT_memory (VNULL, (size_t)header.ny, sizeof (float), GMT_program);
	z_top   = (float *) GMT_memory (VNULL, (size_t)header.nx, sizeof (float), GMT_program);
	z_bot   = (float *) GMT_memory (VNULL, (size_t)header.nx, sizeof (float), GMT_program);
 
	memset ((void *)xy_status, 0, (size_t)(11 * sizeof (int)));
	xy_status[5] = -2;	xy_status[6] = +1;	xy_status[9] = -1;	xy_status[10] = +2;
	
	nx = header.nx;
	nx1 = nx - 1;
	ny1 = header.ny - 1;
	nx2 = nm - nx;
	k = nx2 - nx;
 	for (i = 0; i < header.nx; i++) {	/* Get extra rows using Natural Boundary Cconditions */
		z_top[i] = (float)(2.0 * stack[last][i] - stack[last][nx+i]);
		z_bot[i] = (float)(2.0 * stack[last][nx2+i] - stack[last][k+i]);
	}
 	for (j = 0; j < header.ny; j++) {	/* Get extra columns using Natural Boundary Cconditions */
		k = j * nx;
		z_left[j]  = (float)(2.0 * stack[last][k] - stack[last][k+1]);
		z_right[j] = (float)(2.0 * stack[last][k+nx-1] - stack[last][k+nx-2]);
	}
	
 	/* Search the entire grid, using boundary points when needed */

	for (j = k = 0; j < header.ny; j++) {
		/* k is middle (current) node */
		for (i = 0; i < header.nx; i++, k++) {
		
			if (i == 0) {		/* Along left edge */
				left = z_left[j];
				right = stack[last][k+1];
			}
			else if (i == nx1) {	/* Along right edge */
				right = z_right[j];
				left = stack[last][k-1];
			}
			else {			/* In the middle somewhere */
				left = stack[last][k-1];
				right = stack[last][k+1];
			}
			if (j == 0) {		/* Along top edge */
				top = z_top[i];
				bot = stack[last][k+nx];
			}
			else if (j == ny1) {	/* Along bottom edge */
				top = stack[last][k-nx];
				bot = z_bot[i];
			}
			else {			/* In the middle somewhere */
				top = stack[last][k-nx];
				bot = stack[last][k+nx];
			}
			/* Examine x-slice for local maxima */
			if (stack[last][k] >= left && stack[last][k] > right)
				x_status = 2;
			else if (stack[last][k] > left && stack[last][k] >= right)
				x_status = 2;
			else if (stack[last][k] <= left && stack[last][k] < right)
				x_status = 1;
			else if (stack[last][k] < left && stack[last][k] <= right)
				x_status = 1;
			else
				x_status = 0;
			/* Examine y-slice for local maxima */
			if (stack[last][k] >= top && stack[last][k] > bot)
				y_status = 8;
			else if (stack[last][k] > top && stack[last][k] >= bot)
				y_status = 8;
			else if (stack[last][k] <= top && stack[last][k] < bot)
				y_status = 4;
			else if (stack[last][k] < top && stack[last][k] <= bot)
				y_status = 4;
			else
				y_status = 0;
				
			/* Combine x- and y-status to yield nature of extremum if found:
			 * the x_status | y_status binary OR will yield values from 0 to 10 but
			 * only 5, 6, 9, and 10 are useful here.  These values are mapped via the
			 * xy_status look-up array to yield +-2, +-1, or 0, so that
			 * +2 means a local maximum, -2 a local minimum
			 * +1 means a saddle with local x-maximum, -1 saddle with y-minimum, and
			 * 0 means no extrema found */
			
			z[k] = (float)xy_status[x_status|y_status]; 
		}
	}

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

void grd_FDIST (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i, c, d, prev1, prev2, prev3;
	double a, b, prob;
	
	prev1 = last - 1;
	prev2 = last - 2;
	prev3 = last - 3;
	if (gmtdefs.verbose && constant[prev3] && factor[prev3] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev2] && factor[prev2] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev1] && factor[prev1] == 0.0) fprintf (stderr, "%s: Warning, operand three == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand four == 0 for FDIST!\n", GMT_program);
	for (i = 0; i < nm; i++) {
		a = (constant[prev3]) ? factor[prev3] : stack[prev3][i];
		b = (constant[prev2]) ? factor[prev2] : stack[prev2][i];
		c = irint ((double)((constant[prev1]) ? factor[prev1] : stack[prev1][i]));
		d = irint ((double)((constant[last]) ? factor[last] : stack[last][i]));
		(void) GMT_f_q (a * c, c, b * d, d, &prob);
		stack[prev3][i] = (float)prob;
	}
}
void grd_FLOOR (float *stack[], BOOLEAN *constant, double *factor, int last)
{
	int i;
	double a;
	
	for (i = 0; i < nm; i++) {
		a = (constant[last]) ? factor[last] : stack[last][i];
		stack[last][i] = (float)floor (a);
	}
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

	/* Next look for symbols with special meaning */

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

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

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

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

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

	op = GMT_hash_lookup (choice, hashnode, HASH_SIZE);

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

	return (op);
}

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

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

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