#ifndef lint
static char *RCSid() { return RCSid("$Id: interpol.c,v 1.29 2004/04/13 17:23:58 broeker Exp $"); }
#endif

/* GNUPLOT - interpol.c */

/*[
 * Copyright 1986 - 1993, 1998, 2004   Thomas Williams, Colin Kelley
 *
 * Permission to use, copy, and distribute this software and its
 * documentation for any purpose with or without fee is hereby granted,
 * provided that the above copyright notice appear in all copies and
 * that both that copyright notice and this permission notice appear
 * in supporting documentation.
 *
 * Permission to modify the software is granted, but not the right to
 * distribute the complete modified source code.  Modifications are to
 * be distributed as patches to the released version.  Permission to
 * distribute binaries produced by compiling modified sources is granted,
 * provided you
 *   1. distribute the corresponding source modifications from the
 *    released version in the form of a patch file along with the binaries,
 *   2. add special version identification to distinguish your version
 *    in addition to the base release version number,
 *   3. provide your name and address as the primary contact for the
 *    support of your modified version, and
 *   4. retain our contact information in regard to use of the base
 *    software.
 * Permission to distribute the released version of the source code along
 * with corresponding source modifications in the form of a patch file is
 * granted with same provisions 2 through 4 for binary distributions.
 *
 * This software is provided "as is" without express or implied warranty
 * to the extent permitted by applicable law.
]*/


/*
 * C-Source file identification Header
 *
 * This file belongs to a project which is:
 *
 * done 1993 by MGR-Software, Asgard  (Lars Hanke)
 * written by Lars Hanke
 *
 * Contact me via:
 *
 *  InterNet: mgr@asgard.bo.open.de
 *      FIDO: Lars Hanke @ 2:243/4802.22   (as long as they keep addresses)
 *
 **************************************************************************
 *
 *   Project: gnuplot
 *    Module:
 *      File: interpol.c
 *
 *   Revisor: Lars Hanke
 *   Revised: 26/09/93
 *  Revision: 1.0
 *
 **************************************************************************
 *
 * LEGAL
 *  This module is part of gnuplot and distributed under whatever terms
 *  gnuplot is or will be published, unless exclusive rights are claimed.
 *
 * DESCRIPTION
 *  Supplies 2-D data interpolation and approximation routines
 *
 * IMPORTS
 *  plot.h
 *    - cp_extend()
 *    - structs: curve_points, coordval, coordinate
 *
 *  setshow.h
 *    - samples, axis array[] variables 
 *    - plottypes
 *
 *  proto.h
 *    - solve_tri_diag()
 *    - typedef tri_diag
 *
 * EXPORTS
 *  gen_interp()
 *  sort_points()
 *  cp_implode()
 *
 * BUGS and TODO
 *  I would really have liked to use Gershon Elbers contouring code for
 *  all the stuff done here, but I failed. So I used my own code.
 *  If somebody is able to consolidate Gershon's code for this purpose
 *  a lot of gnuplot users would be very happy - due to memory problems.
 *
 **************************************************************************
 *
 * HISTORY
 * Changes:
 *  Nov 24, 1995  Markus Schuh (M.Schuh@meteo.uni-koeln.de):
 *      changed the algorithm for csplines
 *      added algorithm for approximation csplines
 *      copied point storage and range fix from plot2d.c
 *
 *  Dec 12, 1995 David Denholm
 *      oops - at the time this is called, stored co-ords are
 *      internal (ie maybe log of data) but min/max are in
 *      user co-ordinates. 
 *      Work with min and max of internal co-ords, and
 *      check at the end whether external min and max need to
 *      be increased. (since samples_1 is typically 100 ; we
 *      dont want to take more logs than necessary)
 *      Also, need to take into account which axes are active
 *
 *  Jun 30, 1996 Jens Emmerich
 *      implemented handling of UNDEFINED points
 */

#include "interpol.h"

#include "alloc.h"
#include "axis.h"
#include "contour.h"
#include "graphics.h"
#include "misc.h"
#include "plot2d.h"
/*  #include "setshow.h" */
#include "util.h"


/* in order to support multiple axes, and to simplify ranging in
 * parametric plots, we use arrays to store some things. For 2d plots,
 * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
 * names in plot.h
 */


/*
 * IMHO, code is getting too cluttered with repeated chunks of
 * code. Some macros to simplify, I hope.
 *
 * do { } while(0) is comp.lang.c recommendation for complex macros
 * also means that break can be specified as an action, and it will
 * 
 */


/* store VALUE or log(VALUE) in STORE, set TYPE as appropriate Do
 * OUT_ACTION or UNDEF_ACTION as appropriate. Adjust range provided
 * type is INRANGE (ie dont adjust y if x is outrange). VALUE must not
 * be same as STORE */
/* FIXME 20010610: UNDEF_ACTION is completely unused ??? Furthermore,
 * this is so similar to STORE_WITH_LOG_AND_UPDATE_RANGE() from axis.h
 * that the two should probably be merged.  */
#define STORE_AND_FIXUP_RANGE(store, value, type, min, max, auto,	\
			      out_action, undef_action)			\
do {									\
    store=value;							\
    if (type != INRANGE)						\
	break;  /* don't set y range if x is outrange, for example */	\
    if ((value) < (min)) {						\
       if ((auto) & AUTOSCALE_MIN)					\
	   (min) = (value);						\
       else {								\
	   (type) = OUTRANGE;						\
	   out_action;							\
	   break;							\
       }								\
    }									\
    if ((value) > (max)) {						\
       if ((auto) & AUTOSCALE_MAX)					\
	   (max) = (value);						\
       else {								\
	   (type) = OUTRANGE;						\
	   out_action;							\
       }								\
    }									\
} while(0)

#define UPDATE_RANGE(TEST,OLD,NEW,AXIS)		\
do {						\
    if (TEST)					\
	(OLD) = AXIS_DE_LOG_VALUE(AXIS,NEW);	\
} while(0)

#define spline_coeff_size 4
typedef double spline_coeff[spline_coeff_size];
typedef double five_diag[5];

static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
static int num_curves __PROTO((struct curve_points * plot));
static double *cp_binomial __PROTO((int points));
static void eval_bezier __PROTO((struct curve_points * cp, int first_point,
				 int num_points, double sr, coordval * px,
				 coordval *py, double *c));
static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
static void do_freq __PROTO((struct curve_points *plot,	int first_point, int num_points));		
int compare_points __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));


/*
 * position curve_start to index the next non-UNDEFINDED point,
 * start search at initial curve_start,
 * return number of non-UNDEFINDED points from there on,
 * if no more valid points are found, curve_start is set
 * to plot->p_count and 0 is returned
 */

static int
next_curve(plot, curve_start)
struct curve_points *plot;
int *curve_start;
{
    int curve_length;

    /* Skip undefined points */
    while (*curve_start < plot->p_count
	   && plot->points[*curve_start].type == UNDEFINED) {
	(*curve_start)++;
    };
    curve_length = 0;
    /* curve_length is first used as an offset, then the correct # points */
    while ((*curve_start) + curve_length < plot->p_count
	   && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
	curve_length++;
    };
    return (curve_length);
}


/* 
 * determine the number of curves in plot->points, separated by
 * UNDEFINED points
 */

static int
num_curves(plot)
struct curve_points *plot;
{
    int curves;
    int first_point;
    int num_points;

    first_point = 0;
    curves = 0;
    while ((num_points = next_curve(plot, &first_point)) > 0) {
	curves++;
	first_point += num_points;
    }
    return (curves);
}



/*
 * build up a cntr_struct list from curve_points
 * this funtion is only used for the alternate entry point to
 * Gershon's code and thus commented out
 ***deleted***
 */


/* HBB 990205: rewrote the 'bezier' interpolation routine,
 * to prevent numerical overflow and other undesirable things happening
 * for large data files (num_data about 1000 or so), where binomial
 * coefficients would explode, and powers of 'sr' (0 < sr < 1) become
 * extremely small. Method used: compute logarithms of these
 * extremely large and small numbers, and only go back to the
 * real numbers once they've cancelled out each other, leaving
 * a reasonable-sized one. */

/*
 * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
 *   and stores them into an array which is hooked to sdat.
 * (MGR 1992)
 */
static double *
cp_binomial(points)
    int points;
{
    register double *coeff;
    register int n, k;
    int e;

    e = points;			/* well we're going from k=0 to k=p_count-1 */
    coeff = gp_alloc(e * sizeof(double), "bezier coefficients");

    n = points - 1;
    e = n / 2;
    /* HBB 990205: calculate these in 'logarithmic space',
     * as they become _very_ large, with growing n (4^n) */
    coeff[0] = 0.0;

    for (k = 0; k < e; k++) {
	coeff[k + 1] = coeff[k] + log(((double) (n - k)) / ((double) (k + 1)));
    }

    for (k = n; k >= e; k--)
	coeff[k] = coeff[n - k];

    return (coeff);
}


/* This is a subfunction of do_bezier() for BEZIER style computations.
 * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
 * the double values holding the next x and y coordinates.
 * (MGR 1992)
 */

static void
eval_bezier(cp, first_point, num_points, sr, px, py, c)
    struct curve_points *cp;
    int first_point;		/* where to start in plot->points (to find x-range) */
    int num_points;		/* to determine end in plot->points */
    double sr;			/* position inside curve, range [0:1] */
    coordval *px;		/* OUTPUT: x and y */
    coordval *py;
    double *c;			/* Bezier coefficient array */
{
    unsigned int n = num_points - 1;
    struct coordinate GPHUGE *this_points;

    this_points = (cp->points) + first_point;

    if (sr == 0.0) {
	*px = this_points[0].x;
	*py = this_points[0].y;
    } else if (sr == 1.0) {
	*px = this_points[n].x;
	*py = this_points[n].y;
    } else {
	/* HBB 990205: do calculation in 'logarithmic space',
	 * to avoid over/underflow errors, which would exactly cancel
	 * out each other, anyway, in an exact calculation
	 */
	unsigned int i;
	double lx = 0.0, ly = 0.0;
	double log_dsr_to_the_n = n * log(1 - sr);
	double log_sr_over_dsr = log(sr) - log(1 - sr);

	for (i = 0; i <= n; i++) {
	    double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr);

	    lx += this_points[i].x * u;
	    ly += this_points[i].y * u;
	}

	*px = lx;
	*py = ly;
    }
}

/*
 * generate a new set of coordinates representing the bezier curve and
 * set it to the plot
 */

static void
do_bezier(cp, bc, first_point, num_points, dest)
    struct curve_points *cp;
    double *bc;			/* Bezier coefficient array */
    int first_point;		/* where to start in plot->points */
    int num_points;		/* to determine end in plot->points */
    struct coordinate *dest;	/* where to put the interpolated data */
{
    int i;
    coordval x, y;

    /* min and max in internal (eg logged) co-ordinates. We update
     * these, then update the external extrema in user co-ordinates
     * at the end.
     */

    double ixmin, ixmax, iymin, iymax;
    double sxmin, sxmax, symin, symax;	/* starting values of above */

    x_axis = cp->x_axis;
    y_axis = cp->y_axis;

    ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
    ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
    iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
    iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);

    for (i = 0; i < samples_1; i++) {
	eval_bezier(cp, first_point, num_points,
		    (double) i / (double) (samples_1 - 1),
		    &x, &y, bc);

	/* now we have to store the points and adjust the ranges */

	dest[i].type = INRANGE;
	STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
	STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);

	dest[i].xlow = dest[i].xhigh = dest[i].x;
	dest[i].ylow = dest[i].yhigh = dest[i].y;

	dest[i].z = -1;
    }

    UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
    UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
    UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
    UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
}

/*
 * call contouring routines -- main entry
 */

/* 
 * it should be like this, but it doesn't run. If you find out why,
 * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
 *
 * Well, all this had originally been inside contour.c, so maybe links
 * to functions and of contour.c are broken.
 * ***deleted***
 * end of unused entry point to Gershon's code
 *
 */

/* 
 * Solve five diagonal linear system equation. The five diagonal matrix is
 * defined via matrix M, right side is r, and solution X i.e. M * X = R.
 * Size of system given in n. Return TRUE if solution exist.
 *  G. Engeln-Muellges/ F.Reutter: 
 *  "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
 *  ISBN 3-411-01677-9
 *
 * /  m02 m03 m04   0   0   0   0    .       .       . \   /  x0  \    / r0  \
 * I  m11 m12 m13 m14   0   0   0    .       .       . I   I  x1  I   I  r1  I
 * I  m20 m21 m22 m23 m24   0   0    .       .       . I * I  x2  I = I  r2  I
 * I    0 m30 m31 m32 m33 m34   0    .       .       . I   I  x3  I   I  r3  I
 *      .   .   .   .   .   .   .    .       .       .        .        .
 * \                           m(n-3)0 m(n-2)1 m(n-1)2 /   \x(n-1)/   \r(n-1)/
 * 
 */
static int
solve_five_diag(m, r, x, n)
    five_diag m[];
    double r[], x[];
    int n;
{
    int i;
    five_diag *hv;

    hv = gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");

    hv[0][0] = m[0][2];
    if (hv[0][0] == 0) {
	free(hv);
	return FALSE;
    }
    hv[0][1] = m[0][3] / hv[0][0];
    hv[0][2] = m[0][4] / hv[0][0];

    hv[1][3] = m[1][1];
    hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
    if (hv[1][0] == 0) {
	free(hv);
	return FALSE;
    }
    hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
    hv[1][2] = m[1][4] / hv[1][0];

    for (i = 2; i < n; i++) {
	hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
	hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
	if (hv[i][0] == 0) {
	    free(hv);
	    return FALSE;
	}
	hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
	hv[i][2] = m[i][4] / hv[i][0];
    }

    hv[0][4] = 0;
    hv[1][4] = r[0] / hv[0][0];
    for (i = 1; i < n; i++) {
	hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
    }

    x[n - 1] = hv[n][4];
    x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
    for (i = n - 3; i >= 0; i--)
	x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];

    free(hv);
    return TRUE;
}


/*
 * Calculation of approximation cubic splines
 * Input:  x[i], y[i], weights z[i]
 *         
 * Returns matrix of spline coefficients
 */
static spline_coeff *
cp_approx_spline(plot, first_point, num_points)
    struct curve_points *plot;
    int first_point;		/* where to start in plot->points */
    int num_points;		/* to determine end in plot->points */
{
    spline_coeff *sc;
    five_diag *m;
    double *r, *x, *h, *xp, *yp;
    struct coordinate GPHUGE *this_points;
    int i;

    x_axis = plot->x_axis;
    y_axis = plot->y_axis;

    sc = gp_alloc((num_points) * sizeof(spline_coeff),
				   "spline matrix");

    if (num_points < 4)
	int_error(plot->token, "Can't calculate approximation splines, need at least 4 points");

    this_points = (plot->points) + first_point;

    for (i = 0; i < num_points; i++)
	if (this_points[i].z <= 0)
	    int_error(plot->token, "Can't calculate approximation splines, all weights have to be > 0");

    m = gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");

    r = gp_alloc((num_points - 2) * sizeof(double), "spline right side");
    x = gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
    h = gp_alloc((num_points - 1) * sizeof(double), "spline help vector");

    xp = gp_alloc((num_points) * sizeof(double), "x pos");
    yp = gp_alloc((num_points) * sizeof(double), "y pos");

    /* KB 981107: With logarithmic axis first convert back to linear scale */

    xp[0] = AXIS_DE_LOG_VALUE(x_axis, this_points[0].x);
    yp[0] = AXIS_DE_LOG_VALUE(y_axis, this_points[0].y);
    for (i = 1; i < num_points; i++) {
	xp[i] = AXIS_DE_LOG_VALUE(x_axis, this_points[i].x);
	yp[i] = AXIS_DE_LOG_VALUE(y_axis, this_points[i].y);
	h[i - 1] = xp[i] - xp[i - 1];
    }

    /* set up the matrix and the vector */

    for (i = 0; i <= num_points - 3; i++) {
	r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
		    - (yp[i + 1] - yp[i]) / h[i]);

	if (i < 2)
	    m[i][0] = 0;
	else
	    m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];

	if (i < 1)
	    m[i][1] = 0;
	else
	    m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
		- 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);

	m[i][2] = 2 * (h[i] + h[i + 1])
	    + 6 / this_points[i].z / h[i] / h[i]
	    + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
	    + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];

	if (i > num_points - 4)
	    m[i][3] = 0;
	else
	    m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
		- 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);

	if (i > num_points - 5)
	    m[i][4] = 0;
	else
	    m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
    }

    /* solve the matrix */
    if (!solve_five_diag(m, r, x, num_points - 2)) {
	free(h);
	free(x);
	free(r);
	free(m);
	free(xp);
	free(yp);
	int_error(plot->token, "Can't calculate approximation splines");
    }
    sc[0][2] = 0;
    for (i = 1; i <= num_points - 2; i++)
	sc[i][2] = x[i - 1];
    sc[num_points - 1][2] = 0;

    sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
    for (i = 1; i <= num_points - 2; i++)
	sc[i][0] = yp[i] - 2 / this_points[i].z *
	    (sc[i - 1][2] / h[i - 1]
	     - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
	     + sc[i + 1][2] / h[i]);
    sc[num_points - 1][0] = yp[num_points - 1]
	- 2 / this_points[num_points - 1].z / h[num_points - 2]
	* (sc[num_points - 2][2] - sc[num_points - 1][2]);

    for (i = 0; i <= num_points - 2; i++) {
	sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
	    - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
	sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
    }

    free(h);
    free(x);
    free(r);
    free(m);
    free(xp);
    free(yp);

    return (sc);
}


/*
 * Calculation of cubic splines
 *
 * This can be treated as a special case of approximation cubic splines, with
 * all weights -> infinity.
 *         
 * Returns matrix of spline coefficients
 */
static spline_coeff *
cp_tridiag(plot, first_point, num_points)
    struct curve_points *plot;
    int first_point, num_points;
{
    spline_coeff *sc;
    tri_diag *m;
    double *r, *x, *h, *xp, *yp;
    struct coordinate GPHUGE *this_points;
    int i;

    x_axis = plot->x_axis;
    y_axis = plot->y_axis;
    if (num_points < 3)
	int_error(plot->token, "Can't calculate splines, need at least 3 points");

    this_points = (plot->points) + first_point;

    sc = gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
    m = gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");

    r = gp_alloc((num_points - 2) * sizeof(double), "spline right side");
    x = gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
    h = gp_alloc((num_points - 1) * sizeof(double), "spline help vector");

    xp = gp_alloc((num_points) * sizeof(double), "x pos");
    yp = gp_alloc((num_points) * sizeof(double), "y pos");

    /* KB 981107: With logarithmic axis first convert back to linear scale */

    xp[0] = AXIS_DE_LOG_VALUE(x_axis,this_points[0].x);
    yp[0] = AXIS_DE_LOG_VALUE(y_axis,this_points[0].y);
    for (i = 1; i < num_points; i++) {
	xp[i] = AXIS_DE_LOG_VALUE(x_axis,this_points[i].x);
	yp[i] = AXIS_DE_LOG_VALUE(y_axis,this_points[i].y);
	h[i - 1] = xp[i] - xp[i - 1];
    }

    /* set up the matrix and the vector */

    for (i = 0; i <= num_points - 3; i++) {
	r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
		    - (yp[i + 1] - yp[i]) / h[i]);

	if (i < 1)
	    m[i][0] = 0;
	else
	    m[i][0] = h[i];

	m[i][1] = 2 * (h[i] + h[i + 1]);

	if (i > num_points - 4)
	    m[i][2] = 0;
	else
	    m[i][2] = h[i + 1];
    }

    /* solve the matrix */
    if (!solve_tri_diag(m, r, x, num_points - 2)) {
	free(h);
	free(x);
	free(r);
	free(m);
	free(xp);
	free(yp);
	int_error(plot->token, "Can't calculate cubic splines");
    }
    sc[0][2] = 0;
    for (i = 1; i <= num_points - 2; i++)
	sc[i][2] = x[i - 1];
    sc[num_points - 1][2] = 0;

    for (i = 0; i <= num_points - 1; i++)
	sc[i][0] = yp[i];

    for (i = 0; i <= num_points - 2; i++) {
	sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
	    - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
	sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
    }

    free(h);
    free(x);
    free(r);
    free(m);
    free(xp);
    free(yp);

    return (sc);
}

static void
do_cubic(plot, sc, first_point, num_points, dest)
    struct curve_points *plot;	/* still containes old plot->points */
    spline_coeff *sc;		/* generated by cp_tridiag */
    int first_point;		/* where to start in plot->points */
    int num_points;		/* to determine end in plot->points */
    struct coordinate *dest;	/* where to put the interpolated data */
{
    double xdiff, temp, x, y;
    double xstart, xend;	/* Endpoints of the sampled x range */
    int i, l;
    struct coordinate GPHUGE *this_points;

    /* min and max in internal (eg logged) co-ordinates. We update
     * these, then update the external extrema in user co-ordinates
     * at the end.
     */
    double ixmin, ixmax, iymin, iymax;
    double sxmin, sxmax, symin, symax;	/* starting values of above */

    x_axis = plot->x_axis;
    y_axis = plot->y_axis;

    ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
    ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
    iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
    iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);

    this_points = (plot->points) + first_point;

    l = 0;

    /* HBB 20010727: Sample only across the actual x range, not the
     * full range of input data */
#if SAMPLE_CSPLINES_TO_FULL_RANGE
    xstart = this_points[0].x;
    xend = this_points[num_points - 1].x;
#else
    xstart = GPMAX(this_points[0].x, sxmin);
    xend = GPMIN(this_points[num_points - 1].x, sxmax);

    if (xstart >= xend)
	int_error(plot->token,
		  "Cannot smooth: no data within fixed xrange!");
#endif
    xdiff = (xend - xstart) / (samples_1 - 1);

    for (i = 0; i < samples_1; i++) {
	x = xstart + i * xdiff;

	/* Move forward to the spline interval this point is in */
	while ((x >= this_points[l + 1].x) && (l < num_points - 2))
	    l++;

	/* KB 981107: With logarithmic x axis the values were
         * converted back to linear scale before calculating the
         * coefficients. Use exponential for log x values. */
	temp = AXIS_DE_LOG_VALUE(x_axis, x)
	    - AXIS_DE_LOG_VALUE(x_axis, this_points[l].x);

	/* Evaluate cubic spline polynomial */
	y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];

	/* With logarithmic y axis, we need to convert from linear to
         * log scale now. */
	if (Y_AXIS.log) {
	    if (y > 0.)
		y = AXIS_DO_LOG(y_axis, y);
	    else
		y = symin - (symax - symin);
	}

	dest[i].type = INRANGE;
	STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
	STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);

	dest[i].xlow = dest[i].xhigh = dest[i].x;
	dest[i].ylow = dest[i].yhigh = dest[i].y;

	dest[i].z = -1;

    }

    UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
    UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
    UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
    UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);

}


/*
 * do_freq() is like the other smoothers only in that it
 * needs to adjust the plot ranges. We don't have to copy
 * approximated curves or anything like that.
 */

static void
do_freq(plot, first_point, num_points)
    struct curve_points *plot;	/* still contains old plot->points */
    int first_point;		/* where to start in plot->points */
    int num_points;		/* to determine end in plot->points */
{
    double x, y;
    int i;
    int x_axis = plot->x_axis;
    int y_axis = plot->y_axis;
    struct coordinate GPHUGE *this;

    /* min and max in internal (eg logged) co-ordinates. We update
     * these, then update the external extrema in user co-ordinates
     * at the end.
     */

    double ixmin, ixmax, iymin, iymax;
    double sxmin, sxmax, symin, symax;	/* starting values of above */

    ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
    ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
    iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
    iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);

    this = (plot->points) + first_point;

    for (i=0; i<num_points; i++) {

	x = this[i].x;
	y = this[i].y;

	this[i].type = INRANGE;

	STORE_AND_FIXUP_RANGE(this[i].x, x, this[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
	STORE_AND_FIXUP_RANGE(this[i].y, y, this[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);

	this[i].xlow = this[i].xhigh = this[i].x;
	this[i].ylow = this[i].yhigh = this[i].y;
	this[i].z = -1;
    }

    UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
    UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
    UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
    UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
}


/*
 * Frequency plots have don't need new points allocated; we just need
 * to adjust the plot ranges. Wedging this into gen_interp() would
 * make that code even harder to read.
 */
 
void
gen_interp_frequency(plot)
    struct curve_points *plot;
{

    int i, curves;
    int first_point, num_points;

    curves = num_curves(plot);

    first_point = 0;
    for (i = 0; i < curves; i++) {
	num_points = next_curve(plot, &first_point);
	do_freq(plot, first_point, num_points);
	first_point += num_points + 1;
    }
    return;
}

/*
 * This is the main entry point used for everything except frequencies.
 * As stated in the header, it is fine, but I'm not too happy with it.
 */

void
gen_interp(plot)
    struct curve_points *plot;
{

    spline_coeff *sc;
    double *bc;
    struct coordinate *new_points;
    int i, curves;
    int first_point, num_points;

    curves = num_curves(plot);
    new_points = gp_alloc((samples_1 + 1) * curves * sizeof(struct coordinate),
			  "interpolation table");

    first_point = 0;
    for (i = 0; i < curves; i++) {
	num_points = next_curve(plot, &first_point);
	switch (plot->plot_smooth) {
	case SMOOTH_CSPLINES:
	    sc = cp_tridiag(plot, first_point, num_points);
	    do_cubic(plot, sc, first_point, num_points,
		     new_points + i * (samples_1 + 1));
	    free(sc);
	    break;
	case SMOOTH_ACSPLINES:
	    sc = cp_approx_spline(plot, first_point, num_points);
	    do_cubic(plot, sc, first_point, num_points,
		     new_points + i * (samples_1 + 1));
	    free(sc);
	    break;

	case SMOOTH_BEZIER:
	case SMOOTH_SBEZIER:
	    bc = cp_binomial(num_points);
	    do_bezier(plot, bc, first_point, num_points,
		      new_points + i * (samples_1 + 1));
	    free((char *) bc);
	    break;
	default:		/* keep gcc -Wall quiet */
	    ;
	}
	new_points[(i + 1) * (samples_1 + 1) - 1].type = UNDEFINED;
	first_point += num_points;
    }

    free(plot->points);
    plot->points = new_points;
    plot->p_max = curves * (samples_1 + 1);
    plot->p_count = plot->p_max - 1;

    return;
}

/* 
 * sort_points
 *
 * sort data succession for further evaluation by plot_splines, etc.
 * This routine is mainly introduced for compilers *NOT* supporting the
 * UNIX qsort() routine. You can then easily replace it by more convenient
 * stuff for your compiler.
 * (MGR 1992)
 */

/* HBB 20010720: To avoid undefined behaviour that would be caused by
 * casting functions pointers around, changed arguments to what
 * qsort() *really* wants */
/* HBB 20010720: removed 'static' to avoid HP-sUX gcc bug */
int
compare_points(arg1, arg2)
    SORTFUNC_ARGS arg1;
    SORTFUNC_ARGS arg2;
{
    struct coordinate const *p1 = arg1;
    struct coordinate const *p2 = arg2;

    if (p1->x > p2->x)
	return (1);
    if (p1->x < p2->x)
	return (-1);
    return (0);
}

void
sort_points(plot)
    struct curve_points *plot;
{
    int first_point, num_points;

    first_point = 0;
    while ((num_points = next_curve(plot, &first_point)) > 0) {
	/* Sort this set of points, does qsort handle 1 point correctly? */
	/* HBB 20010720: removed casts -- they don't help a thing, but
	 * may hide problems */
	qsort(plot->points + first_point, num_points,
	      sizeof(struct coordinate), compare_points);
	first_point += num_points;
    }
    return;
}


/*
 * cp_implode() if averaging is selected this function computes the new
 *              entries and shortens the whole thing to the necessary
 *              size
 * MGR Addendum
 */

void
cp_implode(cp)
    struct curve_points *cp;
{
    int first_point, num_points;
    int i, j, k;
    double x = 0., y = 0., sux = 0., slx = 0., suy = 0., sly = 0.;
    /* int x_axis = cp->x_axis; */
    /* int y_axis = cp->y_axis; */
    TBOOLEAN all_inrange = FALSE;

    x_axis = cp->x_axis;                                                   
    y_axis = cp->y_axis;                                                   
    j = 0;
    first_point = 0;
    while ((num_points = next_curve(cp, &first_point)) > 0) {
	k = 0;
	for (i = first_point; i < first_point + num_points; i++) {
	    /* HBB 20020801: don't try to use undefined datapoints */
	    if (cp->points[i].type == UNDEFINED)
	        continue;
	    if (!k) {
		x = cp->points[i].x;
		y = cp->points[i].y;
		sux = cp->points[i].xhigh;
		slx = cp->points[i].xlow;
		suy = cp->points[i].yhigh;
		sly = cp->points[i].ylow;
		all_inrange = (cp->points[i].type == INRANGE);
		k = 1;
	    } else if (cp->points[i].x == x) {
		y += cp->points[i].y;
		sux += cp->points[i].xhigh;
		slx += cp->points[i].xlow;
		suy += cp->points[i].yhigh;
		sly += cp->points[i].ylow;
		if (cp->points[i].type != INRANGE)
		    all_inrange = FALSE;
		k++;
	    } else {
		cp->points[j].x = x;
 		if ( cp->plot_smooth == SMOOTH_FREQUENCY )
		    k = 1;
		cp->points[j].y = y /= (double) k;
		cp->points[j].xhigh = sux / (double) k;
		cp->points[j].xlow = slx / (double) k;
		cp->points[j].yhigh = suy / (double) k;
		cp->points[j].ylow = sly / (double) k;
		/* HBB 20000405: I wanted to use STORE_AND_FIXUP_RANGE
		 * here, but won't: it assumes we want to modify the
		 * range, and that the range is given in 'input'
		 * coordinates.  For logarithmic axes, the overhead
		 * would be larger than the possible gain, so write it
		 * out explicitly, instead:
		 * */
		cp->points[j].type = INRANGE;
		if (! all_inrange) {
		    if (X_AXIS.log) {
			if (x <= -VERYLARGE) {
			    cp->points[j].type = OUTRANGE;
			    goto is_outrange;
			}
			x = AXIS_UNDO_LOG(x_axis, x);
		    }
		    if (((x < X_AXIS.min) && !(X_AXIS.autoscale & AUTOSCALE_MIN))
			|| ((x > X_AXIS.max) && !(X_AXIS.autoscale & AUTOSCALE_MAX))) {
			cp->points[j].type = OUTRANGE;
			goto is_outrange;
		    } 
		    if (Y_AXIS.log) {
			if (y <= -VERYLARGE) {
			    cp->points[j].type = OUTRANGE;
			    goto is_outrange;
			}
			y = AXIS_UNDO_LOG(y_axis, y);
		    }
		    if (((y < Y_AXIS.min) && !(Y_AXIS.autoscale & AUTOSCALE_MIN))
			|| ((y > Y_AXIS.max) && !(Y_AXIS.autoscale & AUTOSCALE_MAX)))
			cp->points[j].type = OUTRANGE;
		is_outrange:
		    ;
		} /* if(! all inrange) */

		j++;		/* next valid entry */
		k = 0;		/* to read */
		i--;		/* from this (-> last after for(;;)) entry */
	    } /* else (same x position) */
	} /* for(points in curve) */

	if (k) {
	    cp->points[j].x = x;
	    if ( cp->plot_smooth == SMOOTH_FREQUENCY )
		k = 1;
	    cp->points[j].y = y /= (double) k;
	    cp->points[j].xhigh = sux / (double) k;
	    cp->points[j].xlow = slx / (double) k;
	    cp->points[j].yhigh = suy / (double) k;
	    cp->points[j].ylow = sly / (double) k;
	    cp->points[j].type = INRANGE;
	    if (! all_inrange) {
		    if (X_AXIS.log) {
			if (x <= -VERYLARGE) {
			    cp->points[j].type = OUTRANGE;
			    goto is_outrange2;
			}
			x = AXIS_UNDO_LOG(x_axis, x);
		    }
		    if (((x < X_AXIS.min) && !(X_AXIS.autoscale & AUTOSCALE_MIN))
			|| ((x > X_AXIS.max) && !(X_AXIS.autoscale & AUTOSCALE_MAX))) {
			cp->points[j].type = OUTRANGE;
			goto is_outrange2;
		    } 
		    if (Y_AXIS.log) {
			if (y <= -VERYLARGE) {
			    cp->points[j].type = OUTRANGE;
			    goto is_outrange2;
			}
			y = AXIS_UNDO_LOG(y_axis, y);
		    }
		    if (((y < Y_AXIS.min) && !(Y_AXIS.autoscale & AUTOSCALE_MIN))
			|| ((y > Y_AXIS.max) && !(Y_AXIS.autoscale & AUTOSCALE_MAX)))
			cp->points[j].type = OUTRANGE;
		is_outrange2:
		    ;
	    }
	    j++;		/* next valid entry */
	}
	/* insert invalid point to separate curves */
	if (j < cp->p_count) {
	    cp->points[j].type = UNDEFINED;
	    j++;
	}
	first_point += num_points;
    }				/* end while */
    cp->p_count = j;
    cp_extend(cp, j);
}