/*--------------------------------------------------------------------
 *	$Id: grdvolume.c,v 1.4.4.4 2002/02/27 17:41:10 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
 *--------------------------------------------------------------------*/
/*
 * grdvolume reads a 2d binary gridded grdfile, and calculates the volume
 * under the surface using exact integration of the bilinear interpolating
 * surface.  As an option, the user may supply a contour value; then the
 * volume is only integrated inside the chosen contour.
 *
 * Author:	Paul Wessel
 * Date:	23-SEP-1997
 * Revised:	02-JUN-1999
 * Version:	3.4.1
 */

#include "gmt.h"

void SW_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da);
void NE_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da);
void SE_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da);
void NW_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da);
void NS_trapezoid (float f[], int ij, int nx, BOOLEAN right, double *dv, double *da);
void EW_trapezoid (float f[], int ij, int nx, BOOLEAN top, double *dv, double *da);
int ors_find_kink (double y[], int n, int mode);
double vol_prism_frac_x (float *z, int ij, int nx, double x0, double x1, double a, double b, double c, double d);
double vol_prism_frac_y (float *z, int ij, int nx, double y0, double y1, double a, double b, double c, double d);
double median3 (double x[]);

main (int argc, char **argv)
{
	int i, j, n = 0, c, ij, nx, ny, ij_inc[4], one_or_zero, k, pos, neg, nc, n_contours, mode = 0, nz = 0, nm; 
	
	BOOLEAN error = FALSE, mapping = FALSE, contour = FALSE, full = FALSE, bad, cut[4];
	BOOLEAN set_base = FALSE;
	
	double take_out, west, east, south, north, dv, da, cval = 0.0, x_size, y_size, m_pr_deg, x_fact, shift = 0.0, fact = 1.0;
	double *area, *vol, *height, clow, chigh, deltac, base, this_base, small;
	
	float *f;

	char *grdfile = CNULL, format[BUFSIZ], map_units = '\0';
	
	struct GRD_HEADER grd;
	
	west = east = south = north = 0.0;
	chigh = deltac = clow = 0.0;
	
	argc = GMT_begin (argc, argv);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
			
				/* Common parameters */
			
				case 'V':
					if (argv[i][2] == 'L' || argv[i][2] == 'l') full = TRUE;
				case 'R':
				case ':':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;


				/* Supplemental parameters */
				
				case 'C':
					contour = TRUE;
					n = sscanf (&argv[i][2], "%lf/%lf/%lf", &clow, &chigh, &deltac);
					if (n == 1) chigh = clow, deltac = 1.0;	/* So calculation of ncontours will yield 1 */
					break;
				case 'L':
					base = (argv[i][2]) ? atof (&argv[i][2]) : GMT_d_NaN;
					set_base = TRUE;
					break;
				case 'S':
					mapping = TRUE;
                                        if (argv[i][2]) map_units = argv[i][2];
					break;
				case 'T':
 					mode = 1;
					break;
				case 'Z':
					nz = (argv[i][2]) ? sscanf (&argv[i][2], "%lf/%lf", &fact, &shift) : -1;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			grdfile = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdvolume %s - Calculating volume under a surface within a contour\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdvolume <grdfile> [-C<cval> or -C<low/high/delta>] [-L<base>] [-S[k]] [-T] [-Rw/s/e/n[r]] [-V] [-Z<fact>[/<shift>]]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<grdfile> is the name of the 2-D binary data set\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C find area and volume inside the <cval> contour\n");
		fprintf (stderr, "\t   OR search using all contours from low to high\n");
		fprintf (stderr, "\t   [Default returns entire area and volume of grid]\n");
		fprintf (stderr, "\t-L Add volume from <base> up to contour [Default is from contour and up only]\n");
		fprintf (stderr, "\t-S Convert degrees to m, append k for km [Default is Cartesian]\n");
		fprintf (stderr, "\t-T Use curvature rather than maximum to find best contour value\n");
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		fprintf (stderr, "\t   Append l for listing of all results (when contour search is selected)\n");
		fprintf (stderr, "\t-Z Subract <shift> and then multiply data by <fact> before processing [1/0]\n");
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!grdfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input grd file\n", GMT_program);
		error++;
	}
	if (contour && !(n == 1 || n == 3)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -C: Must specify 1 or 3 arguments\n", GMT_program);
		error++;
	}
	if (nz < 0 || nz > 2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -Z: Must specify <fact> and optionally <shift>\n", GMT_program);
		error++;
	}
	if (mapping && !(map_units == '\0' || map_units == 'k')) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -S: May append k only\n", GMT_program);
		error++;
	}
	if (set_base && GMT_is_dnan (base)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -L: Must specify base\n", GMT_program);
		error++;
	}

	if (error) exit (EXIT_FAILURE);

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

	if (GMT_read_grd_info (grdfile, &grd)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}
	
	if (!project_info.region_supplied) {	/* No subset asked for */
		west = grd.x_min;
		east = grd.x_max;
		south = grd.y_min;
		north = grd.y_max;
	}
	else if (!project_info.region)	/* Got w/s/e/n, make into w/e/s/n */
		d_swap (south, east);

	one_or_zero = (grd.node_offset) ? 0 : 1;
	nx = irint ( (east - west) / grd.x_inc) + one_or_zero;
	ny = irint ( (north - south) / grd.y_inc) + one_or_zero;
	nm = nx * ny;

	f = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);

	if (GMT_read_grd (grdfile, &grd, f, west, east, south, north, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}

	ij_inc[0] = 0;	ij_inc[1] = 1;	ij_inc[2] = 1 - nx;	ij_inc[3] = -nx;
	x_fact = grd.x_inc;
	if (mapping) {
		m_pr_deg = TWO_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0;
		if (map_units == 'k') m_pr_deg *= 0.001;	/* Use km instead */
		x_size = m_pr_deg * grd.x_inc;
		y_size = m_pr_deg * grd.y_inc;
	}
	else {
		x_size = grd.x_inc;
		y_size = grd.y_inc;
	}

	n_contours = (contour) ? irint ((chigh - clow) / deltac) + 1 : 1;

	height = (double *) GMT_memory (VNULL, (size_t)n_contours, sizeof (double), GMT_program);
	vol = (double *) GMT_memory (VNULL, (size_t)n_contours, sizeof (double), GMT_program);
	area = (double *) GMT_memory (VNULL, (size_t)n_contours, sizeof (double), GMT_program);

	if (!(fact == 1.0 && shift == 0.0)) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Subtracting %lg and multiplying by %lg\n", GMT_program, shift, fact);
		for (k = 0; k < nm; k++) f[k] = (float)((f[k] - shift) * fact);
		grd.z_min = (grd.z_min - shift) * fact;
		grd.z_max = (grd.z_max - shift) * fact;
		if (fact < 0.0) d_swap (grd.z_min, grd.z_max);
	}

	this_base = (set_base) ? base : 0.0;
	small = deltac * 1.0e-6;

	for (c = 0; contour && c < n_contours; c++) {	/* Trace contour, only count volumes inside contours */

		cval = clow + c * deltac;
		take_out = (c == 0) ? cval : deltac;	/* Take out start contour the first time and just the increment subsequent times */

		for (k = 0; k < nm; k++) {
			f[k] -= (float)take_out;		/* Take out the zero value */
			if (f[k] == 0.0) f[k] = (float)small;	/* But we dont want exactly zero, just + or - */
		}
		if (set_base) this_base -= take_out;

		if (set_base && this_base >= 0.0) {
			fprintf (stderr, "%s: Base is > than contour - exiting\n", GMT_program);
			exit (EXIT_FAILURE);
		}

		for (j = 1, ij = grd.nx; j < grd.ny; j++) {

			dv = da = 0.0;	/* Reset these for each row */

			for (i = 0; i < grd.nx-1; i++, ij++) {

				/* Find if a contour goes through this bin */

				for (k = neg = pos = 0, bad = FALSE; !bad && k < 4; k++) {
					(f[ij+ij_inc[k]] <= (float)small) ? neg++ : pos++;
					if (GMT_is_fnan (f[ij+ij_inc[k]])) bad = TRUE;
				}

                                if (bad || neg == 4) continue;	/* Contour not crossing, go to next bin */

				if (pos == 4) {	/* Need entire prism */
					dv += 0.25 * (f[ij] + f[ij+1] + f[ij-nx] + f[ij-nx+1]);
					da += 1.0;
				}
				else {	/* Need partial prisms */

					for (k = nc = 0; k < 4; k++) cut[k] = FALSE;
					if ((f[ij+1] * f[ij]) < 0.0)       nc++, cut[0] = TRUE;	/* Crossing the S border */
					if ((f[ij+1] * f[ij+1-nx]) < 0.0)  nc++, cut[1] = TRUE;	/* Crossing the E border */
					if ((f[ij-nx] * f[ij+1-nx]) < 0.0) nc++, cut[2] = TRUE;	/* Crossing the N border */
					if ((f[ij-nx] * f[ij]) < 0.0)      nc++, cut[3] = TRUE;	/* Crossing the W border */

					if (nc < 2) continue;	/* Can happen if some nodes were 0 and then reset to smal, thus passing the test */

					if (nc == 4) {	/* Saddle scenario */
						if (f[ij] > 0) {	/* Need both SW and NE triangles */
							SW_triangle (f, ij, nx, TRUE, &dv, &da);
							NE_triangle (f, ij, nx, TRUE, &dv, &da);
						}
						else {			/* Need both SE and NW corners */
							SE_triangle (f, ij, nx, TRUE, &dv, &da);
							NW_triangle (f, ij, nx, TRUE, &dv, &da);
						}
						
					}
					else if (cut[0]) {	/* Contour enters at S border ... */
						if (cut[1])	/* and exits at E border */
							SE_triangle (f, ij, nx, (f[ij+1] > 0.0), &dv, &da);
						else if (cut[2])	/* or exits at N border */
							NS_trapezoid (f, ij, nx, f[ij] < 0.0, &dv, &da);
						else			/* or exits at W border */
							SW_triangle (f, ij, nx, (f[ij] > 0.0), &dv, &da);
					}
					else if (cut[1]) {	/* Contour enters at E border */
						if (cut[2])	/* exits at N border */
							NE_triangle (f, ij, nx, (f[ij+1-nx] > 0.0), &dv, &da);
						else			/* or exits at W border */
							EW_trapezoid (f, ij, nx, f[ij] < 0.0, &dv, &da);
					}
					else			/* Contours enters at N border and exits at W */
						NW_triangle (f, ij, nx, (f[ij-nx] > 0.0), &dv, &da);
				}
			}
			ij++;

			if (set_base) dv -= this_base;	/* Adjust for lower starting base */

			/* Allow for shrinking of longitudes with latitude */
			if (mapping) x_fact = x_size * cosd (grd.y_max - (j+0.5) * grd.y_inc);

			dv *= x_fact;
			da *= x_fact;
			vol[c]  += dv;
			area[c] += da;
		}
	}
	if (!contour) {	/* Columns with bilinear tops */
		if (set_base && base >= grd.z_min) {
			fprintf (stderr, "%s: Base is > than minimum z - exiting\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		for (j = 1, ij = grd.nx; j < grd.ny; j++) {
			dv = da = 0.0;
			for (i = 0; i < grd.nx-1; i++, ij++) {
				for (k = 0, bad = FALSE; !bad && k < 4; k++) if (GMT_is_fnan (f[ij+ij_inc[k]])) bad = TRUE;
                                if (bad) continue;
				
				dv += 0.25 * (f[ij] + f[ij+1] + f[ij-nx] + f[ij-nx+1]);
				da += 1.0;
				if (set_base) dv -= this_base;	/* Adjust for lower starting base */
			}
			ij++;

			/* Allow for shrinking of longitudes with latitude */
			if (mapping) x_fact = x_size * cosd (grd.y_max - (j+0.5) * grd.y_inc);

			dv *= x_fact;
			da *= x_fact;
			vol[0]  += dv;
			area[0] += da;
		}
	}
	
	/* Adjust values for actual y dimension */

	for (c = 0; c < n_contours; c++) {

		vol[c]  *= y_size;
		area[c] *= y_size;
		height[c] = (area[c] > 0.0) ? vol[c] / area[c] : GMT_d_NaN;
	}


	/* Find the best contour that gives largest height */

	c = (contour) ? ors_find_kink (height, n_contours, mode) : 0;

	/* Print out final estimates */

        sprintf (format, "%s\t%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);

	if (full) {
		sprintf (format, "%s\t%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		for (c = 0; c < n_contours; c++) fprintf (GMT_stdout, format, clow + c * deltac, area[c], vol[c], height[c]);
	}
	else if (contour) {
        	sprintf (format, "%s\t%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (GMT_stdout, format, clow + c * deltac, area[c], vol[c], height[c]);
	}
	else {
        	sprintf (format, "%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (GMT_stdout, format, area[c], vol[c], height[c]);
	}

	GMT_free ((void *)f);
	GMT_free ((void *)area);
	GMT_free ((void *)vol);
	GMT_free ((void *)height);
	
	GMT_end (argc, argv);
}

void SW_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da)
{	/* Calculates area of a SW-corner triangle */
	/* triangle = TRUE gets triangle, FALSE gives the complementary area */
	double x1, y0;

	x1 = f[ij] / (f[ij] - f[ij+1]);
	y0 = f[ij] / (f[ij] - f[ij-nx]);
	if (triangle) {
		*dv += vol_prism_frac_x (f, ij, nx, 0.0, x1, 0.0, 0.0, -y0 / x1, y0);
		*da += 0.5 * x1 * y0;
	}
	else {
		*dv += 0.25 * (f[ij] + f[ij+1] + f[ij-nx] + f[ij-nx+1]) - vol_prism_frac_x (f, ij, nx, 0.0, x1, 0.0, 0.0, -y0 / x1, y0);
		*da += 1.0 - 0.5 * x1 * y0;
	}
}

void NE_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da)
{	/* Calculates area of a NE-corner triangle */
	/* triangle = TRUE gets triangle, FALSE gives the complementary area */
	double x0, y1, a, x0_1, y1_1;

	x0 = f[ij-nx] / (f[ij-nx] - f[ij+1-nx]);
	y1 = f[ij+1] / (f[ij+1] - f[ij+1-nx]);
	x0_1 = 1.0 - x0;
	y1_1 = y1 - 1.0;
	a = y1_1 / x0_1;
	if (triangle) {
		*dv += vol_prism_frac_x (f, ij, nx, x0, 1.0, a, 1.0 - a * x0, 0.0, 0.0);
		*da -= 0.5 * x0_1 * y1_1;	/* -ve because we need 1 - y1 */
	}
	else {
		*dv += 0.25 * (f[ij] + f[ij+1] + f[ij-nx] + f[ij-nx+1]) - vol_prism_frac_x (f, ij, nx, x0, 1.0, a, 1.0 - a * x0, 0.0, 0.0);
		*da += 1.0 + 0.5 * x0_1 * y1_1;	/* +ve because we need 1 - y1 */
	}
}

void SE_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da)
{	/* Calculates area of a SE-corner triangle */
	/* triangle = TRUE gets triangle, FALSE gives the complementary area */
	double x0, y1, c, x0_1;

	x0 = f[ij] / (f[ij] - f[ij+1]);
	y1 = f[ij+1] / (f[ij+1] - f[ij+1-nx]);
	x0_1 = 1.0 - x0;
	c = y1 / x0_1;
	if (triangle) {
		*dv += vol_prism_frac_x (f, ij, nx, x0, 1.0, 0.0, 0.0, c, -c * x0);
		*da += 0.5 * x0_1 * y1;
	}
	else {
		*dv += 0.25 * (f[ij] + f[ij+1] + f[ij-nx] + f[ij-nx+1]) - vol_prism_frac_x (f, ij, nx, x0, 1.0, 0.0, 0.0, c, -c * x0);
		*da += 1.0 - 0.5 * x0_1 * y1;
	}
}

void NW_triangle (float f[], int ij, int nx, BOOLEAN triangle, double *dv, double *da)
{	/* Calculates area of a NW-corner triangle */
	/* triangle = TRUE gets triangle, FALSE gives the complementary area */
	double x1, y0, y0_1;

	x1 = f[ij-nx] / (f[ij-nx] - f[ij+1-nx]);
	y0 = f[ij] / (f[ij] - f[ij-nx]);
	y0_1 = 1.0 - y0;
	if (triangle) {
		*dv += vol_prism_frac_x (f, ij, nx, 0.0, x1, y0_1 / x1, y0, 0.0, 1.0);
		*da += 0.5 * x1 * y0_1;
	}
	else {
		*dv += 0.25 * (f[ij] + f[ij+1] + f[ij-nx] + f[ij-nx+1]) - vol_prism_frac_x (f, ij, nx, 0.0, x1, y0_1 / x1, y0, 0.0, 1.0);
		*da += 1.0 - 0.5 * x1 * y0_1;
	}
}

void NS_trapezoid (float f[], int ij, int nx, BOOLEAN right, double *dv, double *da)
{	/* Calculates area of a NS trapezoid */
	/* right = TRUE gets the right trapezoid, FALSE gets the left */
	double x0, x1;

	x0 = f[ij] / (f[ij] - f[ij+1]);
	x1 = f[ij-nx] / (f[ij-nx] - f[ij+1-nx]);
	if (right) {	/* Need right piece */
		*dv += vol_prism_frac_y (f, ij, nx, 0.0, 1.0, x1 - x0, x0, 0.0, 1.0);
		*da += 0.5 * (2.0 - x0 - x1);
	}
	else {
		*dv += vol_prism_frac_y (f, ij, nx, 0.0, 1.0, 0.0, 0.0, x1 - x0, x0);
		*da += 0.5 * (x0 + x1);
	}
}

void EW_trapezoid (float f[], int ij, int nx, BOOLEAN top, double *dv, double *da)
{	/* Calculates area of a EW trapezoid */
	/* top = TRUE gets the top trapezoid, FALSE gets the bottom */
	double y0, y1;

	y0 = f[ij] / (f[ij] - f[ij-nx]);
	y1 = f[ij+1] / (f[ij+1] - f[ij+1-nx]);
	if (top) {	/* Need top piece */
		*dv += vol_prism_frac_x (f, ij, nx, 0.0, 1.0, y1 - y0, y0, 0.0, 1.0);
		*da += 0.5 * (2.0 - y0 - y1);
	}
	else {
		*dv += vol_prism_frac_x (f, ij, nx, 0.0, 1.0, 0.0, 0.0, y1 - y0, y0);
		*da += 0.5 * (y0 + y1);
	}
}

/* This function returns the volume bounded by a trapezoid based on two vertical
 * lines x0 and x1 and two horizontal lines y0 = ax +b and y1 = cx + d
 */

double vol_prism_frac_x (float *z, int ij, int nx, double x0, double x1, double a, double b, double c, double d)
{
	double dzdx, dzdy, dzdxy, ca, db, c2a2, d2b2, cdab, v, x02, x12, x03, x04, x13, x14;

	dzdx = (z[ij+1] - z[ij]);
	dzdy = (z[ij-nx] - z[ij]);
	dzdxy = (z[ij-nx+1] + z[ij] - z[ij+1] - z[ij-nx]);

	ca = c - a;
	db = d - b;
	c2a2 = c * c - a * a;
	d2b2 = d * d - b * b;
	cdab = c * d - a * b;
	x02 = x0 * x0;	x03 = x02 * x0;	x04 = x02 * x02;
	x12 = x1 * x1;	x13 = x12 * x1;	x14 = x12 * x12;

	v = (3.0 * dzdxy * c2a2 * (x14 - x04) +
	     4.0 * (2.0 * dzdx * ca + dzdy * c2a2 + 2.0 * dzdxy * cdab) * (x13 - x03) +
	     6.0 * (2.0 * z[ij] * ca + 2.0 * dzdx * db + 2.0 * dzdy * cdab + dzdxy * d2b2) * (x12 - x02) +
	     12.0 * (2.0 * z[ij] * db + dzdy * d2b2) * (x1 - x0)) / 24.0;

	return (v);
}

/* This function returns the volume bounded by a trapezoid based on two horizontal
 * lines y0 and y1 and two vertical lines x0 = ay +b and x1 = cy + d
 */

double vol_prism_frac_y (float *z, int ij, int nx, double y0, double y1, double a, double b, double c, double d)
{
	double dzdx, dzdy, dzdxy, ca, db, c2a2, d2b2, cdab, v, y02, y03, y04, y12, y13, y14;

	dzdx = (z[ij+1] - z[ij]);
	dzdy = (z[ij-nx] - z[ij]);
	dzdxy = (z[ij-nx+1] + z[ij] - z[ij+1] - z[ij-nx]);

	ca = c - a;
	db = d - b;
	c2a2 = c * c - a * a;
	d2b2 = d * d - b * b;
	cdab = c * d - a * b;
	y02 = y0 * y0;	y03 = y02 * y0;	y04 = y02 * y02;
	y12 = y1 * y1;	y13 = y12 * y1;	y14 = y12 * y12;

	v = (3.0 * dzdxy * c2a2 * (y14 - y04) +
	     4.0 * (2.0 * dzdy * ca + dzdx * c2a2 + 2.0 * dzdxy * cdab) * (y13 - y03) +
	     6.0 * (2.0 * z[ij] * ca + 2.0 * dzdy * db + 2.0 * dzdx * cdab + dzdxy * d2b2) * (y12 - y02) +
	     12.0 * (2.0 * z[ij] * db + dzdx * d2b2) * (y1 - y0)) / 24.0;

	return (v);
}

int ors_find_kink (double y[], int n, int mode)
{	/* mode: 0 = find maximum, 1 = find curvature kink */
	int i, ic, im;
	double *c, *f;

	if (mode == 0) {	/* Find maximum value */	

		for (i = im = 0; i < n; i++) if (y[i] > y[im]) im = i;
		return (im);
	}

	/* Calculate curvatures */

	c = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);

	for (i = 1; i < (n-1); i++) c[i] = y[i+1] - 2.0 * y[i] + y[i-1];
	c[0] = c[1];
	if (n > 1) c[n-1] = c[n-2];

	/* Apply 3-point median filter to curvatures  */

	f = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
	for (i = 1; i < (n-1); i++) f[i] = median3 (&c[i-1]);

	/* Find maximum negative filtered curvature */

	for (i = ic = 1; i < (n-1); i++) if (f[i] < f[ic]) ic = i;

	GMT_free ((void *)c);
	GMT_free ((void *)f);

	return (ic);
}

double median3 (double x[])
{

	if (x[0] < x[1]) {
		if (x[2] > x[1]) return (x[1]);
		if (x[2] > x[0]) return (x[2]);
		return (x[0]);
	}
	else {
		if (x[2] > x[0]) return (x[0]);
		if (x[2] < x[1]) return (x[1]);
		return (x[2]);
	}
}