/* spline.c
********************************************************************************
* spline.c -- Routines for cubic spline curve objects
*
* User routines:
* Pro_SplineCurve(series) creates a new curve object by fitting
* cubic splines to the valid points of a
* given series.
*
* Mark Showalter, PDS Ring-Moon Systems Node, April 1998
*******************************************************************************/
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "profile.h"
#include "fortran.h"
/********************
* Type definitions *
********************/
typedef struct ZPRO_SPLINE_STRUCT {
XPRO_CLASS class;
RL_INT4 last;
RL_FLT8 *table0, *table2;
RL_FLT8 x1, x2, dx, k1;
} ZPRO_SPLINE;
/********************
* Static variables *
********************/
static XPRO_CLASS spline_class = {XPRO_SPLINE_CLASS, "cubic-spline", NULL};
/********************************
* Internal function prototypes *
********************************/
static RL_FLT8 ZPro_SplineValue RL_PROTO((RL_VOID *pointer, RL_FLT8 x,
RL_FLT8 *slope));
static RL_FLT8 ZPro_SplineInverse RL_PROTO((RL_VOID *pointer, RL_FLT8 y,
RL_FLT8 x1, RL_FLT8 x2));
static void ZPro_FreeSpline RL_PROTO((RL_VOID *pointer));
static void ZPro_PrintSpline RL_PROTO((RL_VOID *pointer));
/* Spline Toolkit prototypes
* Note that the entire source code for the Spline Toolkit has been appended
* below. The only change has been to "comment out" all of the #include
* directives, since they are handled above.
*/
RL_BOOL Spline_Setup RL_PROTO((RL_INT4 last, RL_BOOL ignore[],
RL_FLT8 y[], RL_FLT8 y2[]));
RL_FLT8 Spline_Value RL_PROTO((RL_FLT8 x, RL_INT4 last, RL_FLT8 y[],
RL_FLT8 y2[], RL_FLT8 *slope));
RL_FLT8 Spline_Invert RL_PROTO((RL_FLT8 value, RL_FLT8 xlo, RL_FLT8 xhi,
RL_INT4 last, RL_FLT8 y[], RL_FLT8 y2[],
RL_BOOL *error_flag));
RL_INT4 Spline_Extrema RL_PROTO((RL_FLT8 x1, RL_FLT8 x2,
RL_INT4 last, RL_FLT8 y[], RL_FLT8 y2[],
RL_INT4 maxcount, RL_FLT8 xlist[],
RL_INT4 types[]));
/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
* Pro_SplineCurve (spline.c)
*$ Abstract:
* This routine generates a curve by fitting cubic splines to the samples
* of a series object.
*$ Keywords:
* PROFILE
* C, PUBLIC, SUBROUTINE
*$ Declarations:
* PRO_OBJECT *Pro_SplineCurve(series)
* PRO_OBJECT *series;
*$ Inputs:
* series pointer to a series object.
*$ Outputs:
* none
*$ Returns:
* pointer to a new spline curve object, or NULL on non-fatal error.
*$ Detailed_description:
* This routine generates a curve by fitting cubic splines to the samples
* of a series object. The domain and coordinate names of the new curve
* match those of the original series.
*
* A cubic spline curve has the property that it is everywhere smooth;
* its value and both its first and second derivatives are continuous.
* This is the optimal method for interpolating between samples of a
* smoothly varying function.
*$ External_references:
* Profile toolkit
*$ Side_effects:
* Memory is allocated.
*$ Examples:
* This snippet of code takes a series that contains only five valid
* samples, and uses it to construct a cubic spline function. Since the
* valid samples are all equal to x squared, the resultant spline function
* will return f(x)=x*x everywhere within the domain [-1,5].
*
* PRO_OBJECT *series, *curve;
* RL_FLT8 table[] = {-1., 0., 1., 4., 98., 99., 25.};
*
* series = Pro_ArraySeries(table, -1, 7, -1., 1.,
* 98., 99., TRUE, FALSE);
* curve = Pro_SplineCurve(series);
*$ Error_handling:
* Profile library error handling is in effect.
*
* Conditions raised:
* PRO_CLASS_ERROR if given object is NULL or is not a series.
* PRO_SETUP_FAILURE if insufficient valid values are found in the
* series object.
* RL_MEMORY_ERROR on memory allocation error.
*$ Limitations:
* none
*$ Author_and_institution:
* Mark R. Showalter
* NASA/Ames Research Center
*$ Version_and_date:
* 1.0: March 1998
*$ Change_history:
* none
*******************************************************************************/
PRO_OBJECT *Pro_SplineCurve(series)
PRO_OBJECT *series;
{
ZPRO_SPLINE *spline;
PRO_OBJECT *new;
RL_INT4 k, k1, k2, nsamples, flag, ix, xcount;
RL_BOOL status;
RL_BOOL *ignore = NULL;
RL_FLT8 *xlist = NULL;
/* Determine sample range in series */
Pro_SeriesIndices(series, &k1, &k2);
nsamples = k2 - k1 + 1;
/* Allocate data structure */
spline = (ZPRO_SPLINE *) XRL_Malloc(sizeof(ZPRO_SPLINE));
if (spline == NULL) return NULL;
spline->table0 = (RL_FLT8 *) XRL_Malloc(nsamples * sizeof(RL_FLT8));
spline->table2 = (RL_FLT8 *) XRL_Malloc(nsamples * sizeof(RL_FLT8));
if (spline->table2 == NULL) goto ERROR;
/* Initialize data structure */
spline->class = spline_class;
spline->last = nsamples - 1;
spline->k1 = (RL_FLT8) k1;
Pro_SeriesSampling(series, &(spline->x1), &(spline->x2), &(spline->dx));
/* Allocate temporary array for series flags */
ignore = (RL_BOOL *) XRL_Malloc(nsamples * sizeof(RL_BOOL));
if (ignore == NULL) goto ERROR;
/* Fill arrays for Spline_Setup() */
for (k = k1; k <= k2; k++) {
spline->table0[k - k1] = Pro_SeriesValue(series, k, &flag);
ignore[k - k1] = (flag != 0);
}
/* Set up spline tables (and override SPLINE_INSUFFICIENT_ARRAY error) */
(void) RL_ClearError1("SPLINE_INSUFFICIENT_ARRAY");
(void) RL_SetErrorType("SPLINE_INSUFFICIENT_ARRAY", RL_RECORD);
status = Spline_Setup(spline->last, ignore, spline->table0, spline->table2);
if (!status) {
if (RL_ClearError1("SPLINE_INSUFFICIENT_ARRAY")) {
(void) sprintf(xpro_message, "curve setup failure\n\
insufficient samples to set up cubic spline curve\n\
series \"%s\" contains fewer than four valid samples",
Pro_ObjectName(series,0));
RL_RaiseError("PRO_SETUP_FAILURE", xpro_message);
}
goto ERROR;
}
XRL_Free((RL_VOID *) ignore);
ignore = NULL;
/* Count and tabulate extrema */
xcount = Spline_Extrema(0., (RL_FLT8) spline->last, spline->last,
spline->table0, spline->table2, 0, NULL, NULL);
xlist = (RL_FLT8 *) XRL_Malloc(xcount * sizeof(RL_FLT8));
if (xcount > 0 && xlist == NULL) goto ERROR;
(void) Spline_Extrema(0., (RL_FLT8) spline->last, spline->last,
spline->table0, spline->table2, xcount, xlist, NULL);
for (ix = 0; ix < xcount; ix++) {
xlist[ix] = Pro_SeriesXValue(series, xlist[ix] + spline->k1);
}
/* Create curve object */
new = XPro_MakeCurve(spline->x1, spline->x2,
xcount, xlist, ZPro_SplineValue, ZPro_SplineInverse,
ZPro_FreeSpline, ZPro_PrintSpline,
(RL_VOID *) spline);
/* Free temporary arrays */
XRL_Free((RL_VOID *) ignore);
XRL_Free((RL_VOID *) xlist);
/* Enslave the series and check for error */
/* if (new != NULL && XPro_EnslaveObject(new, series)) {
Pro_FreeObject(new);
}
*/
/* Transfer coodinate names to new function */
Pro_RenameObject(new, 1, Pro_ObjectName(series,1));
Pro_RenameObject(new, 2, Pro_ObjectName(series,2));
return new;
ERROR:
XRL_Free((RL_VOID *) ignore);
XRL_Free((RL_VOID *) xlist);
ZPro_FreeSpline((RL_VOID *) spline);
return NULL;
}
/*
*******************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_SplineValue(pointer, x, slope)
*
* This internal routine calculates the value of a spline curve.
*
* Inputs:
* pointer pointer to the ZPRO_SPLINE data structure.
* x location at which to evaluate spline curve.
* *slope slope of spline curve at x.
*
* Return: value of spline curve at x.
*******************************************************************************/
static RL_FLT8 ZPro_SplineValue(pointer, x, slope)
RL_VOID *pointer;
RL_FLT8 x, *slope;
{
ZPRO_SPLINE *spline;
RL_FLT8 value;
spline = (ZPRO_SPLINE *) pointer;
value = Spline_Value((x - spline->x1) / spline->dx, spline->last,
spline->table0, spline->table2, slope);
if (slope != NULL) *slope /= spline->dx;
return value;
}
/*
*******************************************************************************
* ZPro_SplineInverse(pointer, y, x1, x2)
*
* This internal routine calculates the inverse of a spline curve
*
* Inputs:
* pointer pointer to the ZPRO_SPLINE data structure.
* y location at which to evaluate spline curve's inverse.
* x1, x2 limits on the search for the inverse.
*
* Return: location at which the value of the spline curve is y.
*
* Errors:
* PRO_EVAL_FAILURE if curve could not be inverted. In this case a
* guess at the nearest x within the specified
* bounds is returned.
*******************************************************************************/
static RL_FLT8 ZPro_SplineInverse(pointer, y, x1, x2)
RL_VOID *pointer;
RL_FLT8 y, x1, x2;
{
ZPRO_SPLINE *spline;
RL_FLT8 value;
RL_BOOL error_flag;
spline = (ZPRO_SPLINE *) pointer;
value = Spline_Invert(y, (x1 - spline->x1) / spline->dx,
(x2 - spline->x1) / spline->dx,
spline->last, spline->table0, spline->table2,
&error_flag);
if (error_flag) {
(void) sprintf(xpro_message, "function evaluation failure\n\
unable to invert cubic spline curve\n\
seeking f(x) = %#g where x is in [%#g,%#g]", y, x1, x2);
RL_RaiseError("PRO_EVALUATION_FAILURE", xpro_message);
}
return spline->dx * value + spline->x1;
}
/*
*******************************************************************************
* ZPro_FreeSpline(pointer)
*
* This internal routine frees the memory used by a ZPRO_SPLINE structure.
*
* Inputs:
* pointer pointer to the ZPRO_SPLINE data structure.
*******************************************************************************/
static void ZPro_FreeSpline(pointer)
RL_VOID *pointer;
{
ZPRO_SPLINE *spline;
spline = (ZPRO_SPLINE *) pointer;
XRL_Free((RL_VOID *) spline->table0);
XRL_Free((RL_VOID *) spline->table2);
XRL_Free(pointer);
}
/*
*******************************************************************************
* ZPro_PrintSpline(pointer)
*
* This internal routine prints out information about a cubic spline curve. It
* is called by Pro_PrintCurve() and is used mainly for debugging.
*
* Input:
* pointer pointer to the ZPRO_SPLINE structure.
*******************************************************************************/
static void ZPro_PrintSpline(pointer)
RL_VOID *pointer;
{
ZPRO_SPLINE *spline;
RL_INT4 i;
spline = (ZPRO_SPLINE *) pointer;
/* Make sure object is not NULL */
if (spline == NULL) {
printf("PRINT ERROR: Cubic spline curve pointer is NULL\n");
return;
}
/* Make sure object is a cubic spline curve */
if (spline->class.id != XPRO_SPLINE_CLASS) {
printf("PRINT ERROR: Object is not a cubic spline curve\n");
return;
}
/* Print object info... */
printf("\nCubic spline curve parameters...\n");
printf("elements = %1d\n", spline->last + 1);
printf("sampling = % #g, % #g, % #g\n", spline->x1, spline->x2,
spline->dx);
for (i = 0; i <= 5 && i <= spline->last; i++)
printf("series[%1d] = % #g, % #g\n", i,
spline->table0[i], spline->table2[i]);
return;
}
/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
* FPro_SplineCurve (spline.c)
*$ Abstract:
* This routine generates a curve by fitting cubic splines to the samples
* of a series object.
*$ Keywords:
* PROFILE
* FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
* integer*4 function FPro_SplineCurve(series)
* integer*4 series
*$ Inputs:
* series FORTRAN pointer to a series object.
*$ Outputs:
* none
*$ Returns:
* pointer to a new spline curve object, or 0 on non-fatal error.
*$ Detailed_description:
* This routine generates a curve by fitting cubic splines to the samples
* of a series object. The domain and coordinate names of the new curve
* match those of the original series.
*
* A cubic spline curve has the property that it is everywhere smooth;
* its value and both its first and second derivatives are continuous.
* This is the optimal method for interpolating between samples of a
* smoothly varying function.
*$ External_references:
* Profile toolkit
*$ Side_effects:
* Memory is allocated.
*$ Examples:
* This snippet of code takes a series that contains only five valid
* samples, and uses it to construct a cubic spline function. Since the
* valid samples are all equal to x squared, the resultant spline function
* will return f(x)=x*x everywhere within the domain [-1,5].
*
* integer*4 series, curve
* real*8 table(7)/ -1.d0, 0.d0, 1.d0, 4.d0, 98.d0, 99.d0, 25.d0/
*
* series = Pro_ArraySeries(table, -1, 7, -1.d0, 1.d0,
* & 98.d0, 99.d0, .TRUE., .FALSE.)
* curve = FPro_SplineCurve(series)
*$ Error_handling:
* Profile library error handling is in effect.
*
* Conditions raised:
* PRO_CLASS_ERROR if given object is NULL or is not a series.
* PRO_SETUP_FAILURE if insufficient valid values are found in the
* series object.
* RL_MEMORY_ERROR on memory allocation error.
*$ Limitations:
* none
*$ Author_and_institution:
* Mark R. Showalter
* NASA/Ames Research Center
*$ Version_and_date:
* 1.0: March 1998
*$ Change_history:
* none
*******************************************************************************/
RL_INT4 FORTRAN_NAME(fpro_splinecurve) (series)
RL_INT4 *series;
{
RL_VOID *ptr1, *ptr2;
RL_INT4 index;
/* Look up series pointer */
ptr1 = FORT_GetPointer(*series);
if (ptr1 == NULL) return 0;
/* Call function */
ptr2 = (RL_VOID *) Pro_SplineCurve((PRO_OBJECT *) ptr1);
if (ptr2 == NULL) return 0;
/* Save new pointer */
index = FORT_AddPointer(ptr2);
if (index == 0) Pro_FreeObject((PRO_OBJECT *) ptr2);
return index;
}
/*
********************************************************************************
* CUBIC SPLINE TOOLKIT
********************************************************************************
*$ Component_name:
* Spline_Setup (spline.c)
*$ Abstract:
* Fills an array with second derivative estimates for use in cubic spline
* interpolation.
*$ Keywords:
* PROFILE, NUMERICAL_ANALYSIS, SPLINE
* C, PUBLIC, SUBROUTINE
*$ Declarations:
* RL_BOOL Spline_Setup(last, ignore, y, y2)
* RL_INT4 last;
* RL_BOOL ignore;
* RL_FLT8 y[], y2[];
*$ Inputs:
* last index for the last tabulated value. Note that this is
* equal to the number of intervals between tabulated
* values, and is one less than the number of tabulated
* values.
* ignore[] TRUE to indicate that the corresponding y-value is
* invalid and should be ignored. You may pass a NULL
* pointer to indicate all y-values are valid.
* y[] tabulated y-values (which must be uniformly spaced).
*$ Outputs:
* y[] any y-values flagged with ignore[] == TRUE are filled in
* with an interpolated value.
* y2[] table of second derivatives, scaled by a factor to
* facilitate interpolation.
*$ Returns:
* TRUE if spline setup was successful; FALSE if an error occurred.
*$ Detailed_description:
* This function fills an array with second derivative estimates for use in
* determining cubic splines. These values may then be used as an input to
* function Spline_Value() for determining interpolated values. This
* interpolation scheme insures that a function and its first and second
* derivatives are everywhere continuous.
*
* In this implementation, the tabulated values must be uniformly spaced in
* the independent variable x. We treat this spacing as unity, which
* results in no loss of generality. The implementation does permit
* embedded invalid values (indicated by TRUE values in the ignore[]
* array); these values are then filled in by cubic spline interpolation so
* that the returned array contains only valid values.
*
* The procedure employed is documented within the code. It is based on
* the algorithm found in:
* Hornbeck, Robert W. Numerical Methods. Quantum Publishers,
* Inc., New York, 1975, pp. 47--50.
* The procedure is quite similar to that found in Numerical Recipes
* subroutine SPLINE, for the special case where x-values are uniformly
* spaced. See:
* Press, William H., et al. Numerical Recipes. Cambridge
* University Press, Cambridge, 1986, pp. 86--89.
*$ External_references:
* none
*$ Side_effects:
* A temporary data array is allocated and then deallocated.
*$ Examples:
* none
*$ Error_handling:
* This routine returns FALSE if an error occurred. It also uses standard
* RingLib error handling. Conditions raised:
*
* RL_MEMORY_ERROR on memory allocation failure.
* SPLINE_INSUFFICIENT_ARRAY if the y-array contains fewer than 4 valid
* values.
*$ Limitations:
* As noted above, the tabulated values must be uniformly spaced in the
* independent variable x. A minimum of four tabulated values must be
* given.
*$ Author_and_institution:
* Mark R. Showalter
* PDS Ring-Moon Systems Node, NASA/Ames Research Center
*$ Version_and_date:
* 1.0: October 1996
*$ Change_history:
* none
*******************************************************************************/
/*
#include <stdlib.h>
#include "spline.h"
#include "rl_memory.h"
#include "rl_error.h"
*/
RL_BOOL Spline_Setup(last, ignore, y, y2)
RL_INT4 last;
RL_BOOL ignore[];
RL_FLT8 y[], y2[];
{
RL_INT4 *i, j, ii, i1, i2, n;
RL_BOOL hasgaps;
RL_FLT8 *coeff, *rhval, dxj0, dxj1, dyj0, dyj1, a, b, r,
x, ytemp[2], y2temp[2];
/* Allocate necessary arrays */
i = (RL_INT4 *) XRL_Malloc((last+1) * sizeof(RL_INT4));
coeff = (RL_FLT8 *) XRL_Malloc((last+1) * sizeof(RL_FLT8));
rhval = (RL_FLT8 *) XRL_Malloc((last+1) * sizeof(RL_FLT8));
if (rhval == NULL) {
XRL_Free((RL_VOID *) i);
XRL_Free((RL_VOID *) coeff);
return FALSE;
}
/*******************************************************************************
* We assume that our sequence of values y[i] is uniformly spaced in x but that
* some points may be missing or invalid, as indicated by a non-zero value for
* ignore[i]. As a first step, we derive an index array i[j] to convert to a
* non-uniform sequence in which all values are valid.
*******************************************************************************/
/* Set up index table that skips over invalid entries */
hasgaps = FALSE;
j = 0;
for (ii=0; ii<=last; ii++) {
if (ignore != NULL && ignore[ii]) {
hasgaps = TRUE;
} else {
i[j] = ii;
j++;
}
}
n = j;
/* Check for valid number of entries */
if (n < 4) {
RL_RaiseError("SPLINE_INSUFFICIENT_ARRAY",
"spline array has insufficient valid values");
XRL_Free((RL_VOID *) i);
XRL_Free((RL_VOID *) coeff);
XRL_Free((RL_VOID *) rhval);
return FALSE;
}
/*******************************************************************************
* Cubic splines require tabulations of the function's second derivative y''
* in addition to the function's value y. The requirement that first derivatives
* be continuous leads to a sequence of coupled equations:
* dx[j] y''[j-1] + 2 (dx[j+1]+dx[j]) y''[j] + dx[j+1] y''[j+1]
* = 6 (dy[j+1]/dx[j+1] - dy[j]/dx[j])
* where
* dx[j] = x[j] - x[j-1];
* dy[j] = y[j] - y[j-1]
* (Press et al., eq. (3.3.7), p. 87). In addition, we use the following
* boundary conditions at the end points:
* y''[0] = y''[1]
* y''[N-2] = y''[N-1]
* These guarantee that the interpolation will match a polynomial up to second
* order exactly.
*
* Instead of returning the actual second derivatives yy'', we return the values
* divided by 6; this makes the interpolation function Spline_Value() faster.
*******************************************************************************/
/* Set up coefficients for first constraint based on the boundary condition:
* y''[0] + coeff[0] y''[1] = rhval[0]
*/
coeff[0] = -1.;
rhval[0] = 0.;
/* Set up coefficients for each new constraint in sequence... */
dxj1 = i[1] - i[0];
dyj1 = y[i[1]] - y[i[0]];
for (j=1; j < n-1; j++) {
/* Calculate differences */
dxj0 = dxj1;
dyj0 = dyj1;
dxj1 = i[j+1] - i[j];
dyj1 = y[i[j+1]] - y[i[j]];
/* Set up coefficients for next constraint:
* dx[j] y''[j-1] + 2 (dx[j+1]+dx[j]) y''[j] + dx[j+1] y''[j+1]
* = (dy[j+1]/dx[j+1] - dy[j]/dx[j])
* rewritten as...
* y''[j-1] + a y''[j] + b y''[j+1] = r[j]
*/
a = 2 * (dxj1 + dxj0) / dxj0;
b = dxj1 / dxj0;
r = (dyj1/dxj1 - dyj0/dxj0) / dxj0;
/* Now subtract off previous constraint and then rescale
* y''[j-1] + a y''[j] + b y''[j+1] = r
* - y''[j-1] + coeff[j-1] y''[j] = rhval[j-1]
* ------------------------------------------------------------------
* (a - coeff[j-1]) y''[j] + b y''[j+1] = (r - rhval[j-1])
*/
coeff[j] = b / (a - coeff[j-1]);
rhval[j] = (r - rhval[j-1]) / (a - coeff[j-1]);
}
/* At this point we have N-1 equations in N unknowns, of the form:
* y''[j] + coeff[j] y''[j+1] = rhval[j]
* for j = 0 to N-2.
*
* Use the last boundary condition to solve for y''[N-1]
* y''[N-2] + coeff[N-2] y''[N-1] = rhval[N-2]
* y''[N-2] + (-1) y''[N-1] = 0
*/
y2[i[n-1]] = rhval[n-2] / (coeff[n-2] + 1.);
/* Now back-substitute to solve for remaining values */
for (j = n-2; j >= 0; j--) {
y2[i[j]] = rhval[j] - coeff[j] * y2[i[j+1]];
}
/*******************************************************************************
* Now we have to fill in the missing points so the returned array has no
* invalid values.
*******************************************************************************/
if (hasgaps) {
j = 1;
for (ii = 0; ii <= last; ii++) {
if (ignore[ii]) {
while (i[j] < ii && j < n-1) j++;
i1 = i[j-1];
i2 = i[j];
ytemp[0] = y[i1];
ytemp[1] = y[i2];
y2temp[0] = y2[i1];
y2temp[1] = y2[i2];
x = (ii - i1) / ((RL_FLT8) (i2 - i1));
y[ii] = Spline_Value(x, 1, ytemp, y2temp, NULL);
y2[ii] = y2temp[0] + x * (y2temp[1] - y2temp[0]);
}
}
}
/* Free the temporary arrays */
XRL_Free((RL_VOID *) i);
XRL_Free((RL_VOID *) coeff);
XRL_Free((RL_VOID *) rhval);
return TRUE;
}
/*
********************************************************************************
*$ Component_name:
* Spline_Value (spline.c)
*$ Abstract:
* Returns a function value and (optionally) the local derivative based on
* cubic spline interpolation.
*$ Keywords:
* PROFILE, NUMERICAL_ANALYSIS, SPLINE
* C, PUBLIC, SUBROUTINE
*$ Declarations:
* RL_FLT8 Spline_Value(x, last, y, y2, slope)
* RL_FLT8 x, y[], y2[], *slope;
* RL_INT4 last;
*$ Inputs:
* x location at which to interpolate, assuming that the
* tabulated values start at x=0 and are spaced by dx=1.
* last last index into the array, equal to one less than the
* dimension of y0[] and y2[].
* y0[0...last] table of y-values.
* y2[0...last] table of second derivatives returned by Spline_Setup().
* slope a value of NULL will suspend the calculation of the
* tabulated function's derivative.
*$ Outputs:
* *slope first derivative of the tabulated function at the point
* selected. (Not calculated if slope==NULL). Note that
* this value is based upon a spacing of unity between
* tabulated values; it should be divided by dx if (dx!=1).
*$ Returns:
* interpolated value of the function.
*$ Detailed_description:
* This function returns a function value and (optionally) the local
* derivative based on cubic spline interpolation. It requires a set of
* tabulate y-values, along with the second derivatives of y as returned by
* Spline_Setup(). The returned value is an extrapolation from the
* endpoints when x is outside the table limits.
*
* The procedure employed is based on the algorithm found in:
* Hornbeck, Robert W. Numerical Methods. Quantum Publishers,
* Inc., New York, 1975, pp. 47--50.
* It is quite similar to that found in Numerical Recipes subroutine SPLINT
* for the special case where x-values are uniformly spaced. See:
* Press, William H., et al. Numerical Recipes. Cambridge
* University Press, Cambridge, 1986, pp. 86--89.
*$ External_references:
* none
*$ Side_effects:
* none
*$ Examples:
* none
*$ Error_handling:
* none
*$ Limitations:
* none
*$ Author_and_institution:
* Mark R. Showalter
* PDS Ring-Moon Systems Node, NASA/Ames Research Center
*$ Version_and_date:
* 1.0: October 1996
*$ Change_history:
* none
*******************************************************************************/
/*
#include <math.h>
#include "spline.h"
*/
RL_FLT8 Spline_Value(x, last, y, y2, slope)
RL_FLT8 x, y[], y2[], *slope;
RL_INT4 last;
{
RL_INT4 index;
RL_FLT8 xlo, xhi, xlo2, xhi2, value;
index = (RL_INT4) floor(x);
if (index < 0) index = 0;
if (index >= last) index = last - 1;
y += index;
y2 += index;
xlo = x - index;
xhi = 1. - xlo;
/*******************************************************************************
* The equations below simplify for the case y2[0] == y2[1]. This situation
* arises frequently because the cubic spline solution forces this condition at
* the end points. A simplified form of the solution avoids some loss of
* precision for x values far outside the range (0...last).
*******************************************************************************/
if (y2[0] != y2[1]) {
xlo2 = xlo * xlo;
xhi2 = xhi * xhi;
value = xhi*y[0] + xlo*y[1] +
xhi*(xhi2 - 1.)*y2[0] + xlo*(xlo2 - 1.)*y2[1];
if (slope != NULL)
*slope = y[1] - y[0] +
(3.*xlo2 - 1.)*y2[1] - (3.*xhi2 - 1.)*y2[0];
}
else {
value = xhi*y[0] + xlo*(y[1] + 3.*(xlo - 1.)*y2[0]);
if (slope != NULL)
*slope = y[1] - y[0] + (6.*xlo - 3.)*y2[0];
}
return value;
}
/*
********************************************************************************
*$ Component_name:
* Spline_Invert (spline.c)
*$ Abstract:
* Finds the index that maps to a particular tabulated spline value.
*$ Keywords:
* PROFILE, NUMERICAL_ANALYSIS, SPLINE
* C, PUBLIC, SUBROUTINE
*$ Declarations:
* RL_FLT8 Spline_Invert(value, xlo, xhi, last, y, y2, error_flag)
* RL_FLT8 value, xlo, xhi, y[], y2[];
* RL_INT4 last;
* RL_BOOL *error_flag;
*$ Inputs:
* value spline y-value for which the corresponding index is
* sought.
* xlo lower bound on the desired solution, or -HUGE_VAL for an
* unbounded search.
* xhi upper bound on the desired solution, or HUGE_VAL for an
* unbounded search.
* y[0...last] table of y-values as passed to Spline_Setup().
* y2[0...last] table of second derivatives as returned by
* Spline_Setup().
*$ Outputs:
* *error_flag TRUE if spline could not be inverted; FALSE otherwise.
*$ Returns:
* the index value that maps nearest to the given spline y-value.
*$ Detailed_description:
* This function returns the index that maps nearest to a particular
* tabulated spline value, within specified limits. It is the inverse of
* function Spline_Value(). The solution is found by Newton's Method.
*$ External_references:
* Spline_Value().
*$ Side_effects:
* none
*$ Examples:
* none
*$ Error_handling:
* *error_flag is set to TRUE if a solution was not found.
*$ Limitations:
* The nearest endpoint will be returned if the y-value sought does not
* fall within the selected limits. Convergence is not guaranteed if the
* spline is not monotonic between the given limits. Also, one has no
* control over which occurrence of the y-value is returned when more than
* one solution is present.
*$ Author_and_institution:
* Mark R. Showalter
* PDS Ring-Moon Systems Node, NASA/Ames Research Center
*$ Version_and_date:
* 1.0: October 1996
*$ Change_history:
* none
*******************************************************************************/
/*
#include <math.h>
#include "spline.h"
*/
#define SMALL_STEP 1.e-8 /* When steps drop below this cutoff, we will
consider an increase in the magnitude of the
step to be considered evidence for
convergence. This guarantees termination of
the loop. */
RL_FLT8 Spline_Invert(value, xlo, xhi, last, y, y2, error_flag)
RL_FLT8 value, xlo, xhi, y[], y2[];
RL_INT4 last;
RL_BOOL *error_flag;
{
RL_FLT8 ylo, yhi, xtest, ytest, dydx, resid, dx, absdx, prev_absdx;
RL_INT4 was_xlo, was_xhi;
/* Anticipate success */
*error_flag = FALSE;
/* Test for valid boundaries */
if (xlo >= xhi) {
*error_flag = TRUE;
return xlo;
}
/* Take an initial guess at x by linear interpolation between endpoints.
* However, for semi-bounded searches, just start near the finite limit; for
* completely unbounded searches, use linear interpolation between tabulated
* endpoints. */
if (xlo != -HUGE_VAL) {
if (xhi != HUGE_VAL) {
ylo = Spline_Value(xlo, last, y, y2, NULL);
yhi = Spline_Value(xhi, last, y, y2, NULL);
xtest = xlo + (value - ylo) * (xhi-xlo)/(yhi-ylo);
}
else {
xtest = xlo + last;
}
}
else {
if (xhi != HUGE_VAL) {
xtest = xhi - last;
}
else {
xtest = last * (value - y[0]) / (y[last] - y[0]);
}
}
/* Initialize flags to indicate when we hit bounds */
was_xlo = FALSE;
was_xhi = FALSE;
/* Avoid termination on first pass */
absdx = SMALL_STEP;
/*******************************************************************************
* Iterate using Newton's Method
*******************************************************************************/
while (1) {
/* Make sure x falls within bounds */
if (xtest < xlo) {
if (was_xlo) {
*error_flag = TRUE;
return xlo;
}
was_xlo = TRUE;
xtest = xlo;
}
if (xtest > xhi) {
if (was_xhi) {
*error_flag = TRUE;
return xhi;
}
was_xhi = TRUE;
xtest = xhi;
}
/* Evaluate spline function and derivative */
ytest = Spline_Value(xtest, last, y, y2, &dydx);
/* Simplest test for convergence */
resid = ytest - value;
if (resid == 0.) return xtest;
/* At local extremum, jump beyond more distant endpoint and try again */
if (dydx == 0.) {
if (xtest - xlo < xhi - xtest) xtest = xhi + 1.;
else xtest = xlo - 1.;
continue;
}
/* Generate next guess */
dx = resid / dydx;
xtest -= dx;
/* If changes stop shrinking, we have convergence */
prev_absdx = absdx;
absdx = fabs(dx);
if (absdx < SMALL_STEP && absdx >= prev_absdx) return xtest;
}
}
/*
********************************************************************************
*$ Component_name:
* Spline_Extrema (spline.c)
*$ Abstract:
* Locates all of the extrema within a cubic spline function.
*$ Keywords:
* PROFILE, NUMERICAL_ANALYSIS, SPLINE
* C, PUBLIC, SUBROUTINE
*$ Declarations:
* RL_INT4 Spline_Extrema(x1, x2, last, y, y2,
* maxcount, xlist, types)
* RL_INT4 last, maxcount, types[];
* RL_FLT8 x1, x2, y[], y2[], xlist[]
*$ Inputs:
* x1, x2 range of x-values within which to locate all extrema.
* last last index of y and y2 tables (one less than number of
* elements).
* y[0...last] table of y-values.
* y2[0...last] table of second derivatives, returned by Spline_Setup().
* maxcount dimensioned number of available elements in xlist[] and
* types[].
* xlist buffer where extremum locations should be placed, or
* NULL to suspend return of locations.
* types buffer where extremum types should be placed, or NULL
* to suspend return of types.
*$ Outputs:
* xlist[0...] locations of all extrema, up to first maxcount.
* types[0...] -1 for minima; 0 for inflection points; +1 for maxima.
*$ Returns:
* The number of local extrema within the interval of interest (which may
* be > maxcount).
*$ Detailed_description:
* This function locates the extrema within a single interval of cubic
* spline interpolation. Note that the interval includes the lower
* endpoint but not the upper endpoint. Inflection points are ignored.
*$ External_references:
* none
*$ Side_effects:
* none
*$ Examples:
* To locate and save all the extrema in a spline function...
*
* total = Spline_Extrema(-HUGE_VAL, HUGE_VAL, y, y2, 0, NULL, NULL);
* xlist = (RL_FLT8) XRL_Malloc(total * sizeof(RL_FLT8));
* types = (RL_INT4) XRL_Malloc(total * sizeof(RL_INT4));
* Spline_Extrema(-HUGE_VAL, HUGE_VAL, y, y2, total, xlist, types);
*$ Error_handling:
* none
*$ Limitations:
* none
*$ Author_and_institution:
* Mark R. Showalter
* PDS Ring-Moon Systems Node, NASA/Ames Research Center
*$ Version_and_date:
* 1.0: January 1997
*$ Change_history:
* none
*******************************************************************************/
/*
#include <math.h>
#include "spline.h"
*/
#define ISIGN(x) (x < 0. ? -1 : (x == 0. ? 0 : 1))
RL_INT4 Spline_Extrema(x1, x2, last, y, y2, maxcount, xlist, types)
RL_INT4 last, maxcount, types[];
RL_FLT8 x1, x2, y[], y2[], xlist[];
{
RL_INT4 count, i1, i2, imin, imax, nx, j, type;
RL_FLT8 a, b, c, discr, q, xs[2], temp, xtest, xtrue, xprev, second;
/*******************************************************************************
* Loop through all relevant intervals of spline
*******************************************************************************/
xprev = -HUGE_VAL;
count = 0;
imin = (x1 < 1.) ? 0 : (RL_INT4) floor(x1);
imax = (x2 > last) ? last : (RL_INT4) ceil(x2);
for (i1 = imin; i1 < imax; i1++) {
i2 = i1 + 1;
/*******************************************************************************
* Locate the extrema
*
* We wish to solve for the x-values at which the slope of the cubic spline is
* zero. We have
* slope = y[i2] - y[i1] +
* (3*x*x - 1)*y2[i2] - (3*(1-x)^2 - 1)*y2[i1]
* This can be solved via the quadratic formula.
*******************************************************************************/
/* Generate quadratic coefficients */
a = 3. * (y2[i2] - y2[i1]);
b = 6. * y2[i1];
c = y[i2] - y[i1] - y2[i2] - 2.*y2[i1];
/* Solve quadratic formula (See Numerical Methods, p. 145) */
nx = 0;
discr = b*b - 4.*a*c;
if (discr >= 0.) {
if (b > 0.) q = -0.5 * (b + sqrt(discr));
else q = -0.5 * (b - sqrt(discr));
if (q != 0.) {xs[nx] = c/q; nx++;}
if (a != 0.) {xs[nx] = q/a; nx++;}
}
/* Put pair of solutions in increasing order and check for duplicates */
if (nx == 2) {
if (xs[0] == xs[1]) {
nx = 1;
}
else if (xs[0] > xs[1]) {
temp = xs[0];
xs[0] = xs[1];
xs[1] = temp;
}
}
/*******************************************************************************
* Test each extremum
*******************************************************************************/
for (j=0; j<nx; j++) {
/* Exclude extrema outside the interval */
xtest = xs[j];
if (xtest < 0. && i1 != 0) continue;
if (xtest > 1. && i2 != last) continue;
xtrue = xtest + i1;
if (xtrue < x1) continue;
if (xtrue > x2) continue;
if (xtrue <= xprev) continue;
/* Test for minimum, maximum or inflection */
second = (a+a)*xtest + b;
type = -ISIGN(second);
/* When an inflection point falls on an interval boundary, we can look
at neighboring intervals for more information */
if (type == 0) {
if (xtrue == i1 && i1 > 0 &&
ISIGN(y2[i1-1]) * ISIGN(y2[i2]) == 1)
type = -ISIGN(y2[i1-1]);
if (xtrue == i2 && i2 < last &&
ISIGN(y2[i1]) * ISIGN(y2[i2+1]) == 1)
type = -ISIGN(y2[i2+1]);
}
/* Save this extremum in list, if possible */
if (count < maxcount) {
if (xlist != NULL) xlist[count] = xtrue;
if (types != NULL) types[count] = type;
}
count++;
xprev = xtrue;
}
}
return count;
}
/*******************************************************************************
*/