.

1 files changed, 0 insertions, 735 deletions

diff --git a/src/spline.c b/src/spline.c
deleted file mode 100644
index 0189dec8f..000000000
--- a/src/spline.c
+++ /dev/null
@@ -1,735 +0,0 @@
-/* spline.c -- Do spline interpolation. */
-
-/******************************************************************************
-	David L. Fox
-	2140 Braun Dr.
-	Golden, CO 80401
-
-This program has been placed in the public domain by its author.
-
-Version Date		Change
-1.1	17 Dec 1985	Modify getdata() to realloc() more memory if needed.
-1.0	14 May 1985
-
-spline [options] [file]
-
-Spline reads pairs of numbers from file (or the standard input, if file is missing).
-The pairs are interpreted as abscissas and ordinates of a function.  The output of
-spline consists of similar pairs generated from the input by interpolation with
-cubic splines. (See R. W. Hamming, Numerical Methods for Scientists and Engineers,
-2nd ed., pp. 349ff.)  There are sufficiently many points in the output to appear
-smooth when plotted (e.g., by graph(1)).  The output points are approximately evenly
-spaced and include the input points.
-
-There may be one or more options of the form: -o [argument [arg2] ].
-
-The available options are:
-
--a	Generate abscissa values automatically.  The input consists of a list of
-	ordinates.  The generated abscissas are evenly spaced with spacing given by
-	the argument, or 1 if the next argument is not a number.
-
--f	Set the first derivative of the spline function at the left and right end
-	points to the first and second arguments following -f.  If only one numerical
-	argument follows -f then that value is used for both left and right end points.
-
--k	The argument of the k option is the value of the constant used to calculate
-	the boundary values by y''[0] = ky''[1] and y''[n] = ky''[n-1].  The default
-	value is k = 0.
-
--m	The argument gives the maximum number of input data points.  The default is 1000.
-	If the input contains more than this many points additional memory will be
-	allocated.  This option may be used to slightly increase efficiency for large
-	data sets or reduce the amount of memory required for small data sets.
-
--n	The number of output points is given by the argument.  The default value is 100.
-
--p	The splines used for interpolation are forced to be periodic, i.e. y'[0] = y'[n].
-	The first and last ordinates should be equal.
-
--s	Set the second derivative of the spline function at the left and right end
-	points to the first and second arguments following -s.  If only one numerical
-	argument follows -s then that value is used for both left and right end points.
-
--x	The argument (and arg2, if present) are the lower (and upper) limits on the
-	abscissa values in the output.  If the x option is not specified these limits
-	are calculated from the data.  Automatic abscissas start at the lower limit
-	(default 0).
-
-The data need not be monotonic in x but all the x values must be distinct.
-Non-numeric data in the input is ignored.
-******************************************************************************/
-
-#include	<a:stdio.h>
-#include	<ctype.h>
-
-/* The constant DOUBLE is a compile time switch.
-	If #define DOUBLE appears here double pecision variables are
-	used to store all data and parameters.
-	Otherwise, float variables are used for the data.
-	For most smoothing and and interpolation applications single
-	precision is more than adequate.
-	Double precision is used to solve the system of linear equations
-	in either case.
-*/
-/* #define DOUBLE	*/
-
-#ifdef	DOUBLE
-# define double	real;		/* Type used for data storage. */
-# define IFMT	"%F"		/* Input format. */
-# define OFMT	"%18.14g %18.14g\n"	/* Output format. */
-#else
-# define float	real;		/* Type used for data storage. */
-# define IFMT	"%f"		/* Input format. */
-# define OFMT	"%8.5g %8.5g\n"	/* Output format. */
-#endif
-
-/* Numerical constants: These may be machine and/or precision dependent. */
-#define HUGE	(1.e38)		/* Largest floating point number. */
-#define EPS	1.e-5		/* Test for zero when solving linear equations. */
-
-/* Default parameters */
-#define NPTS	1000
-#define NINT	100
-#define DEFSTEP	1.	/* Default step size for automatic abcissas. */
-#define DEFBVK	0.	/* Default boundary value constant. */
-
-/* Boundary condition types. */
-#define EXTRAP	0	/* Extrapolate second derivative:
-			   y''(0) = k * y''(1), y''(n) = k * y''(n-1) */
-#define FDERIV	1	/* Fixed first derivatives y'(0) and y'(n). */
-#define SDERIV	2	/* Fixed second derivatives y''(0) and y''(n). */
-#define PERIOD	3	/* Periodic:  derivatives equal at end points. */
-
-/* Token types for command line processing. */
-#define OPTION	1
-#define NUMBER	2
-#define OTHER	3
-
-/* Define error numbers. */
-#define MEMERR	1
-#define NODATA	2
-#define BADOPT	4
-#define BADFILE	5
-#define XTRAARG	6
-#define XTRAPTS	7
-#define SINGERR	8
-#define DUPX	9
-#define BADBC	10
-#define RANGERR	11
-
-/* Constants and flags are global. */
-int aflag = FALSE;	/* Automatic abscissa flag. */
-real step = DEFSTEP;	/* Spacing for automatic abscissas. */
-int bdcode = EXTRAP;	/* Type of boundary conditions:
-			0 extrapolate f'' with coeficient k.
-			1 first derivatives given
-			2 second derivatives given
-			3 periodic */
-real leftd = 0.,	/* Boundary values of derivatives. */
-   rightd = 0.;
-real k = DEFBVK;	/* Boundary value constant. */
-int rflag = 0;		/* 0: take range from data, 1: min. given, 2: min. & max. given. */
-real xmin = HUGE, xmax;	/* Output range. */
-unsigned nint = NINT;	/* Number of intervals in output. */
-unsigned maxpts = NPTS;	/* Maximun number of points. */
-unsigned nknots = 0;	/* Number of input points. */
-int sflag = FALSE;	/* If TRUE data must be sorted. */
-char *datafile;		/* Name of data file */
-
-int xcompare();		/* Function to compare x values of two data points. */
-
-struct pair {
-	real x, y;
-	} *data;	/* Points to list of data points. */
-
-struct bandm {
-	double diag, upper, rhs;
-	} *eqns;	/* Points to elements of band matrix used to calculate
-			coefficients of the splines. */
-
-main(argc, argv)
-int argc;
-char **argv;
-{
-	setup(argc, argv);
-	if (nknots == 1) {	/* Cannot interpolate with one data point.  */
-		printf(OFMT,data->x, data->y);
-		exit(0);
-	}
-	if (sflag)	/* Sort data if needed. */
-		qsort(data, nknots, sizeof(struct pair), xcompare);
-	calcspline();
-	interpolate();
-}
-
-/* xcompare -- Compare abcissas of two data points (for qsort). */
-xcompare(arg1, arg2)
-struct pair *arg1, *arg2;
-{
-	if (arg1->x > arg2->x)
-		return(1);
-	else if (arg1->x < arg2->x)
-		return(-1);
-	else
-		return(0);
-}
-
-/* error -- Print error message and abort fatal errors. */
-error(errno, auxmsg)
-int errno;
-char *auxmsg;
-{	char *msg;
-	int fatal, usemsg;
-	static char *usage =
-	"usage: spline [options] [file]\noptions:\n",
-	*options = "-a spacing\n-k const\n-n intervals\n-m points\n-p\n-x xmin xmax\n";
-
-	fatal = TRUE;	/* Default is fatal error. */
-	usemsg = FALSE;	/* Default no usage message. */
-	fprintf(stderr, "spline: ");
-	switch(errno) {
-	case MEMERR:
-		msg = "not enough memory for %u data points\n";
-		break;
-	case NODATA:
-		msg = "data file %s empty\n";
-		break;
-	case BADOPT:
-		msg = "unknown option: %c\n";
-		usemsg = TRUE;
-		break;
-	case BADFILE:
-		msg = "cannot open file: %s\n";
-		break;
-	case XTRAARG:
-		msg = "extra argument ignored: %s\n";
-		fatal = FALSE;
-		usemsg = TRUE;
-		break;
-	case XTRAPTS:
-		fatal = FALSE;
-		msg = "%s";
-		break;
-	case DUPX:
-		msg = "duplicate abcissa value: %s\n";
-		break;
-	case RANGERR:
-		msg = "xmax < xmin not allowed %s\n";
-		break;
-	/* The following errors "can't happen." */
-	/* If they occur some sort of bug is indicated. */
-	case SINGERR:
-		msg = "singular matrix encountered %s\n";
-		break;
-	case BADBC:
-		msg = "internal error: bad boundary value code %d\n";
-		break;
-	default:
-		fprintf(stderr, "unknown error number: %d\n", errno);
-		exit(1);
-	}
-	fprintf(stderr, msg, auxmsg);
-	if (usemsg)
-		fprintf(stderr,"%s%s", usage, options);
-	if (fatal)
-		exit(1);
-}
-
-/* setup -- Initalize all constants and read data. */
-setup(argc, argv)
-int argc;
-char **argv;
-{	char *malloc();
-	FILE fdinp,		/* Source of input. */
-		doarg();
-
-	fdinp = doarg(argc, argv);	/* Process command line arguments. */
-
-	/* Allocate memory for data and band matrix of coefficients. */
-	if ((data = malloc(maxpts*sizeof(struct pair))) == NULL) {
-		error(MEMERR, (char *)maxpts);
-	}
-
-	getdata(fdinp);		/* Read data from fdinp. */
-	if (fdinp != stdin)
-		fclose(fdinp);	/* Close input data file. */
-	if (nknots == 0) {
-		error(NODATA, datafile);
-	}
-	/* Allocate memory for calculation of spline coefficients. */
-	if ((eqns = malloc((nknots+1)*sizeof(struct bandm))) == NULL) {
-		error(MEMERR, (char *)nknots);
-	}
-}
-
-/* doarg -- Process arguments. */
-FILE
-doarg(argc, argv)
-int argc;
-char **argv;
-{	int i, type;
-	double atof();
-	char *s, str[15];
-	FILE fdinp;
-
-	s = argv[i=1];
-	type = gettok(&i, &s, argv, str);
-	do {
-		if (type == OPTION) {
-			switch(*str) {
-			case 'a':	/* Automatic abscissa. */
-				aflag = TRUE;
-				rflag = rflag < 2 ? 1 : rflag;
-				if (xmin == HUGE)	/* Initialize xmin, if needed. */
-					xmin = 0.;
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					if ((step = atof(str)) <= 0.)
-						error(RANGERR,"");
-					type = gettok(&i, &s, argv, str);
-				}
-				break;
-			case 'f':	/* Fix first derivative at boundaries. */
-				bdcode = FDERIV;
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					leftd = atof(str);
-					if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-						rightd = atof(str);
-						type = gettok(&i, &s, argv, str);
-					}
-					else {
-						rightd = leftd;
-					}
-				}
-				break;
-			case 'k':	/* Set boundary value constant. */
-				bdcode = EXTRAP;
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					k = atof(str);
-					type = gettok(&i, &s, argv, str);
-				}
-				break;
-			case 'm':	/* Set number of intervals in output. */
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					maxpts = (unsigned)atoi(str);
-					type = gettok(&i, &s, argv, str);
-				}
-				break;
-			case 'n':	/* Set number of intervals in output. */
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					nint = (unsigned)atoi(str);
-					type = gettok(&i, &s, argv, str);
-				}
-				break;
-			case 'p':	/* Require periodic interpolation function. */
-				bdcode = PERIOD;
-				type = gettok(&i, &s, argv, str);
-				break;
-			case 's':	/* Fix second derivative at boundaries. */
-				bdcode = SDERIV;
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					leftd = atof(str);
-					if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-						rightd = atof(str);
-						type = gettok(&i, &s, argv, str);
-					}
-					else {
-						rightd = leftd;
-					}
-				}
-				break;
-			case 'x':	/* Set range of x values. */
-				rflag = 1;
-				if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-					xmin = atof(str);
-					if ((type = gettok(&i, &s, argv, str)) == NUMBER) {
-						xmax = atof(str);
-						rflag = 2;
-						type = gettok(&i, &s, argv, str);
-						if (xmax <= xmin)
-							error(RANGERR, "");
-					}
-				}
-				break;
-			default:
-				error(BADOPT, (char *)argv[i][1]);
-			}
-		}
-		else {
-			if (argc > i) {
-				datafile = argv[i];
-				if ((fdinp = fopen(argv[i++], "r")) == NULL) {
-					error(BADFILE, argv[i-1]);
-				}
-				if (argc > i)
-					error(XTRAARG, argv[i]);
-			}
-			else
-				fdinp = stdin;
-			break;
-		}
-	} while (i < argc);
-	return fdinp;
-}
-
-/* gettok -- Get one token from command line, return type. */
-gettok(indexp, locp, argv, str)
-int *indexp;	/* Pointer to index in argv array. */
-char **locp;	/* Pointer to current location in argv[*indexp]. */
-char **argv;
-char *str;
-{	char *s;
-	char *strcpy(), *strchr();
-	int type;
-
-	s = *locp;
-	while (isspace(*s) || *s == ',')
-		++s;
-	if (*s == '\0')		/* Look at next element in argv. */
-		s = argv[++*indexp];
-	if (*s == '-' && isalpha(s[1])) {
-		/* Found an option. */
-		*str = *++s;
-		str[1] = '\0';
-		++s;
-		type = OPTION;
-	}
-	else if (is_number(s)) {
-		while (is_number(s))
-			*str++ = *s++;
-		*str = '\0';
-		type = NUMBER;
-	}
-	else {
-		strcpy(str, s);
-		s = strchr(s, '\0');
-		type = OTHER;
-	}
-	*locp = s;
-	return(type);
-}
-
-/* is_number -- Return TRUE if argument is the ASCII representation of a number. */
-is_number(string)
-char *string;
-{
-	if (isdigit(*string) ||
-	   *string == '.'    ||
-	   *string == '+'   ||
-	   (*string == '-' && (isdigit(string[1]) || string[1] == '.')))
-		return(TRUE);
-	else
-		return(FALSE);
-}
-
-/* getdata -- Read data points from fdinp. */
-getdata(fdinp)
-FILE fdinp;
-{	int n, i;
-	real lastx, min, max;
-	struct pair *dp;
-	char msg[60], *realloc();
-
-	nknots = 0;
-	lastx = -HUGE;
-	dp = data;	/* Point to head of list of data. */
-	min = HUGE;
-	max = -HUGE;
-	do {
-		if (aflag) {	/* Generate abcissa.  */
-			dp->x = xmin + nknots*step;
-			n = 0;
-		}
-		else {		/* Read abcissa. */
-			while ((n = (fscanf(fdinp,IFMT,&dp->x))) == 0)
-				;	/* Skip non-numeric input. */
-		}
-		if (n == 1) {
-			if (min > dp->x)
-				min = dp->x;
-			if (max < dp->x)
-				max = dp->x;
-			if (lastx > dp->x) {		/* Check for monotonicity. */
-				sflag = TRUE;	/* Data must be sorted. */
-			}
-			lastx = dp->x;
-		}
-		/* Read ordinate. */
-		while ((n = (fscanf(fdinp, IFMT, &dp->y))) == 0)
-			;	/* Skip non-numeric input. */
-		if (++nknots >= maxpts) {		/* Too many points, allocate more memory. */
-			if ((data = realloc(data, (maxpts *= 2)*sizeof(struct pair))) == NULL) {
-				error(MEMERR, (char *)maxpts);
-			}
-			dp = data + nknots;
-			sprintf(msg, "more than %d input points, more memory allocated\n",
-			   maxpts/2);
-			error(XTRAPTS,msg);
-		}
-		else
-			++dp;
-	} while (n == 1);
-	--nknots;
-	if (aflag) {	/* Compute maximum ordinate. */
-		max = xmin + (nknots-1)*step;
-	}
-	if (rflag < 2) {
-		xmax = max;
-		if (rflag < 1)
-			xmin = min;
-	}
-}
-
-/* calcspline -- Calculate coeficients of spline functions. */
-calcspline()
-{
-	calccoef();
-	solvband();
-}
-
-/* calccoef -- Calculate coefficients of linear equations determining spline functions. */
-calccoef()
-{	int i, j;
-	struct bandm *ptr, tmp1, tmp2;
-	double h, h1;
-	char str[10];
-
-	ptr = eqns;
-	/* Initialize first row of matrix. */
-	if ((h1 = data[1].x - data[0].x) == 0.) {
-		sprintf(str, "%8.5g", data[0].x);
-		error(DUPX, str);
-	}
-	switch(bdcode) {	/* First equation depends on boundary conditions. */
-	case EXTRAP:
-		ptr->upper = -k;
-		ptr->diag = 1.;
-		ptr->rhs = 0.;
-		break;
-	case FDERIV:		/* Fixed first derivatives at ends. */
-		ptr->upper = 1.;
-		ptr->diag = 2.;
-		h = data[1].x - data[0].x;
-		ptr->rhs = 6.*((data[1].y - data[0].y)/(h*h) - leftd/h);
-		break;
-	case SDERIV:
-		ptr->upper = 0.;
-		ptr->diag = 1.;
-		ptr->rhs = leftd;
-		break;
-	case PERIOD:		/* Periodic splines. */
-		ptr->upper = 1.;
-		ptr->diag = 2.;
-		break;
-	default:
-		error(BADBC, (char *) bdcode);
-	}
-	++ptr;
-
-	/* Initialize rows 1 to n-1, sub-diagonal band is assumed to be 1. */
-	for (i=1; i	< nknots-1; ++i, ++ptr) {
-		h = h1;
-		if ((h1 = data[i+1].x - data[i].x) == 0.) {
-			sprintf(str, "%8.5g", data[i].x);
-			error(DUPX, str);
-		}
-		ptr->diag = 2.*(h + h1)/h;
-		ptr->upper = h1/h;
-		ptr->rhs = 6.*((data[i+1].y-data[i].y)/(h*h1) -
-		   (data[i].y - data[i-1].y)/(h*h));
-	}
-	/* Initialize last row. */
-	switch(bdcode) {
-	case EXTRAP:		/* Standard case. */
-		ptr->diag = 1.;
-		ptr->upper = -k;	/* Upper holds actual sub-diagonal value. */
-		ptr->rhs = 0.;
-		break;
-	case FDERIV:		/* Fixed first derivatives at ends. */
-		ptr->upper = 1.;
-		ptr->diag = 2.;
-		h = data[nknots-1].x - data[nknots-2].x;
-		ptr->rhs = 6.*((data[nknots-2].y - data[nknots-1].y)/(h*h) + rightd/h);
-		break;
-	case SDERIV:
-		ptr->upper = 0.;
-		ptr->diag = 1.;
-		ptr->rhs = rightd;
-		break;
-	case PERIOD:	/* Use periodic boundary conditions. */
-		/* First and last row are not in tridiagonal form. */
-		h = data[1].x - data[0].x;
-		h1 = data[nknots-1].x - data[nknots-2].x;
-		ptr->diag = 1.;
-		ptr->upper = 0.;
-		tmp1.diag = -1.;
-		tmp1.upper = 0;
-		tmp1.rhs = 0.;
-		tmp2.diag = 2.*h1/h;
-		tmp2.upper = h1/h;
-		tmp2.rhs = 6.*((data[1].y - data[0].y)/(h*h) -
-			   (data[nknots-1].y - data[nknots-2].y)/(h1*h));
-		/* Transform periodic boundary equations to tri-diagonal form. */
-		for (i = 1; i < nknots - 1; ++i) {
-			tmp1.upper /= tmp1.diag;
-			tmp1.rhs /= tmp1.diag;
-			ptr->diag /= tmp1.diag;
-			ptr->upper /= tmp1.diag;
-			tmp1.diag = tmp1.upper - eqns[i].diag;
-			tmp1.upper = -eqns[i].upper;
-			tmp1.rhs -= eqns[i].rhs;
-			tmp2.upper /= tmp2.diag;
-			tmp2.rhs /= tmp2.diag;
-			eqns->diag /= tmp2.diag;
-			eqns->upper /= tmp2.diag;
-			tmp2.diag = tmp2.upper - eqns[nknots-1-i].diag/eqns[nknots-1-i].upper;
-			tmp2.upper = -1./eqns[nknots-1-i].upper;
-			tmp2.rhs -= eqns[nknots-1-i].rhs/eqns[nknots-1-i].upper;
-		}
-		/* Add in remaining terms of boundary condition equation. */
-		ptr->upper += tmp1.diag;
-		ptr->diag += tmp1.upper;
-		ptr->rhs = tmp1.rhs;
-		eqns->diag += tmp2.upper;
-		eqns->upper += tmp2.diag;
-		eqns->rhs = tmp2.rhs;
-		break;
-	default:
-		error(BADBC, (char *) bdcode);
-	}
-}
-
-/* solvband -- Solve band matrix for spline functions. */
-solvband()
-{	int i, flag;
-	struct bandm *ptr;
-	double k1;
-	double fabs();
-	int fcompare();
-
-	ptr = eqns;
-	flag = FALSE;
-	/* Make a pass to triangularize matrix. */
-	for (i=1; i < nknots - 1; ++i, ++ptr) {
-		if (fabs(ptr->diag) < EPS) {
-		/* Near zero on diagonal, pivot. */
-			if (fabs(ptr->upper) < EPS)
-				error(SINGERR, "");
-			flag = TRUE;
-			ptr->diag = i;		/* Keep row index in diag.
-						Actual value of diag is always 1. */
-			if (i == nknots - 2) {
-				flag = 2;
-				/* Exchange next to last and last rows. */
-				k1 = ptr->rhs/ptr->upper;
-				if (fabs((ptr+1)->upper) < EPS)
-					error(SINGERR, "");
-				ptr->rhs = (ptr+1)->rhs/(ptr+1)->upper;
-				ptr->upper = (ptr+1)->diag/(ptr+1)->upper;
-				(ptr+1)->upper = 0.;
-				(ptr+1)->rhs = k1;
-			}
-			else {
-				ptr->rhs = (ptr+1)->rhs - (k1 = ptr->rhs/ptr->upper)*(ptr+1)->diag;
-				ptr->upper = (ptr+1)->upper;	/* This isn't super-diagonal element
-								but rather one to its right.
-								Must watch for this when
-								back substituting. */
-				(++ptr)->diag = i-1;
-				++i;
-				ptr->upper = 0;
-				ptr->rhs = k1;
-				(ptr+1)->rhs -= ptr->rhs;
-			}
-		}
-		else {
-			ptr->upper /= ptr->diag;
-			ptr->rhs /= ptr->diag;
-			ptr->diag = i-1;		/* Used to reorder solution if needed. */
-			(ptr+1)->diag -= ptr->upper;
-			(ptr+1)->rhs -= ptr->rhs;
-		}
-	}
-	/* Last row is a special case. */
-	if (flag != 2) {
-		/* If flag == 2 last row is already computed. */
-		ptr->upper /= ptr->diag;
-		ptr->rhs /= ptr->diag;
-		ptr->diag = ptr - eqns;
-		++ptr;
-		ptr->diag -= (ptr-1)->upper * ptr->upper;
-		ptr->rhs -= (ptr-1)->rhs * ptr->upper;
-		ptr->rhs /= ptr->diag;
-		ptr->diag = ptr - eqns;
-	}
-
-	/* Now make a pass back substituting for solution. */
-	--ptr;
-	for ( ; ptr >= eqns; --ptr) {
-		if ((int)ptr->diag != ptr - eqns) {
-			/* This row and one above have been exchanged in a pivot. */
-			--ptr;
-			ptr->rhs -= ptr->upper * (ptr+2)->rhs;
-		}
-		else
-			ptr->rhs -= ptr->upper * (ptr+1)->rhs;
-	}
-	if (flag) {	/* Undo reordering done by pivoting. */
-		qsort(eqns, nknots, sizeof(struct bandm), fcompare);
-	}
-}
-
-/* fcompare -- Compare two floating point numbers. */
-fcompare(arg1, arg2)
-real *arg1, *arg2;
-{
-	if (arg1 > arg2)
-		return(1);
-	else if (arg1 < arg2)
-		return(-1);
-	else
-		return(0);
-}
-
-/* interpolate -- Do spline interpolation. */
-interpolate()
-{	int i;
-	struct pair *dp;
-	struct bandm *ep;
-	real h, xi, yi, hi, xu, xl, limit;
-
-	h = (xmax - xmin)/nint;
-	ep = eqns;
-	dp = data + 1;
-	for (xi = xmin; xi < xmax + 0.25*h; xi += h) {
-		while (dp->x < xi && dp < data + nknots - 1) {
-			++dp;	/* Skip to correct interval. */
-			++ep;
-		}
-		if (dp < data + nknots - 1 && dp->x < xmax)
-			limit = dp->x;
-		else
-			limit = xmax;
-		for ( ; xi < limit - 0.25*h; xi += h) {
-			/* Do interpolation. */
-			hi = dp->x - (dp-1)->x;
-			xu = dp->x - xi;
-			xl = xi - (dp-1)->x;
-			yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
-			     (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
-			printf(OFMT, xi, yi);
-		}
-		if (limit != dp->x) {	/* Interpolate. */
-			hi = dp->x - (dp-1)->x;
-			xu = dp->x - xmax;
-			xl = xmax - (dp-1)->x;
-			yi = ((ep+1)->rhs*xl*xl/(6.*hi) + dp->y/hi - (ep+1)->rhs*hi/6.)*xl +
-			     (ep->rhs*xu*xu/(6.*hi) + (dp-1)->y/hi - ep->rhs*hi/6.)*xu;
-			printf(OFMT, xmax, yi);
-		}
-		else {		/* Print knot. */
-			printf(OFMT, dp->x, dp->y);
-			xi = dp->x;
-		}
-	}
-}