#include "interpolate.h"
/* returns piecewise-cubic approximation of f(x) at X */
double pwhermite
(
const double *x, /* array of strictly increasing points */
const double *fx, /* array of known values of f(x) */
const double *fpx, /* array of known values of f'(x) */
int n, /* number of values in x, fx, and fpx */
double X /* value of x at which to compute f(x) */
)
{
double a[4], m[3], dd, dx;
int k;
/*
** Locate the interval on which X lies. If this interval is
** not found by k = n - 3, [x[n - 2], x[n - 1]] is used.
*/
for (k = 0; k < n - 2; k++) {
if (X < x[k + 1])
break;
}
/*
** Compute values to be used in the computation of
** coefficients for the cubic interpolating polynomial.
*/
dx = x[k + 1] - x[k];
dd = (fx[k + 1] - fx[k]) / dx;
/*
** Compute the coefficients of the interpolating polynomial
** and determine multipliers based on center calculations
** involving X and the interpolation points.
*/
a[0] = fx[k];
a[1] = fpx[k];
a[2] = (dd - fpx[k]) / dx;
a[3] = (fpx[k + 1] + fpx[k] - (2.0 * dd)) / (dx * dx);
m[0] = X - x[k];
m[1] = m[0] * (X - x[k]);
m[2] = m[1] * (X - x[k + 1]);
/*
** Compute p3(x), the cubic interpolating polynomial,
** to obtain a value for f(x) at X.
*/
return a[0] + (a[1] * m[0]) + (a[2] * m[1]) + (a[3] * m[2]);
}