/* 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;
}

/*******************************************************************************
*/