/* lspline.c
********************************************************************************
* lspline.c -- routines for linear spline functions
*
* User routines:
* Pro_LSplineFunc(series) creates a function that returns the
* value interpolated between the nearest
* entries in a series.
*
* Mark Showalter & Neil Heather, 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_LSPLINE_STRUCT {
XPRO_CLASS class;
PRO_OBJECT *series;
RL_INT4 kmax;
} ZPRO_LSPLINE;
/********************
* Static variables *
********************/
static XPRO_CLASS lspline_class = {XPRO_LSPLINE_CLASS, "linear-spline", NULL};
/********************************
* Internal function prototypes *
********************************/
static RL_FLT8 ZPro_EvalLSpline RL_PROTO((RL_VOID *pointer, RL_FLT8 x));
static void ZPro_FreeLSpline RL_PROTO((RL_VOID *pointer));
static void ZPro_PrintLSpline RL_PROTO((RL_VOID *pointer));
/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
* Pro_LSplineFunc (lspline.c)
*$ Abstract:
* This routine generates a new linear spline function object, which
* performs linear interpolation between the nearest valid series values.
*$ Keywords:
* PROFILE, FUNCTION
* C, PUBLIC, SUBROUTINE
*$ Declarations:
* PRO_OBJECT *Pro_LSplineFunc(series)
* PRO_OBJECT *series;
*$ Inputs:
* series pointer to the series object
*$ Outputs:
* none
*$ Returns:
* a pointer to a new linear spline object, or NULL on non-fatal error.
*$ Detailed_description:
* This routine generates a new linear spline function object, which
* performs linear interpolation between the nearest valid series values.
*$ External_references:
* Profile toolkit
*$ Side_effects:
* Memory is allocated. The new object retains a link to the original
* series.
*$ Examples:
* This snippet of code takes a series that contains only three valid
* samples, and uses it to construct an linear spline function. The
* function returns f(x)=x for 1 <= x <= 5.
*
* PRO_OBJECT *series, *func;
* RL_FLT8 table[] = {1., 2., 98., 99., 5.};
*
* series = Pro_ArraySeries(table, 1, 5, 1., 1.,
* 98., 99., TRUE, FALSE);
* func = Pro_LSplineFunc(series);
*$ Error_handling:
* Profile library error handling is in effect.
*
* Conditions raised:
* PRO_CLASS_ERROR if series object is NULL or not a series.
* PRO_SETUP_FAILURE if series contains less than two valid values.
* 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_LSplineFunc(series)
PRO_OBJECT *series;
{
ZPRO_LSPLINE *lspline;
PRO_OBJECT *new;
RL_INT4 k, k1, k2, flag;
RL_FLT8 x1, x2, temp;
/* Skip over invalid entries at each end of series */
(void) Pro_SeriesIndices(series, &k1, &k2);
for (k = k1; k <= k2; k++) {
(void) Pro_SeriesValue(series, k, &flag);
if (flag == 0) break;
}
k1 = k;
for (k = k2; k >= k1; k--) {
(void) Pro_SeriesValue(series, k, &flag);
if (flag == 0) break;
}
k2 = k;
if (k1 >= k2) {
(void) sprintf(xpro_message, "function setup failure\n\
insufficient samples to set up linear spline function\n\
series \"%s\" contains fewer than two valid samples",
Pro_ObjectName(series,0));
RL_RaiseError("PRO_SETUP_FAILURE", xpro_message);
}
/* Determine domain of function */
x1 = Pro_SeriesXValue(series, (RL_FLT8) k1);
x2 = Pro_SeriesXValue(series, (RL_FLT8) k2);
if (x1 > x2) {temp = x1; x1 = x2; x2 = temp;}
/* Allocate and initialize structure */
lspline = (ZPRO_LSPLINE *) XRL_Malloc(sizeof(ZPRO_LSPLINE));
if (lspline == NULL) return NULL;
lspline->class = lspline_class;
lspline->series = series;
lspline->kmax = k2;
/* Create function object */
new = XPro_MakeFunc(x1, x2,
ZPro_EvalLSpline, ZPro_FreeLSpline, ZPro_PrintLSpline,
lspline);
/* Enslave the series and check for error */
if (new != NULL && XPro_EnslaveObject(new, series)) {
Pro_FreeObject(new);
}
/* Transfer coordinate names to new function */
Pro_RenameObject(new, 1, Pro_ObjectName(series,1));
Pro_RenameObject(new, 2, Pro_ObjectName(series,2));
return new;
}
/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_EvalLSpline(pointer, x)
*
* This internal routine calculates the value of a linear spline function.
*
* Inputs:
* pointer pointer to a ZPRO_LSPLINE data structure.
* x location at which to evaluate linear spline function.
*
* Return: value of function at x.
*
* Errors:
* PRO_INTERPOLATION_ERROR on interpolation failure due to lack of valid
* data above or below x-value. In this situation, the
* nearest valid y-value is returned if possible; otherwise
* zero. In general, this should not happen because the
* series is checked for valid end points when it is
* created.
*******************************************************************************/
static RL_FLT8 ZPro_EvalLSpline(pointer, x)
RL_VOID *pointer;
RL_FLT8 x;
{
ZPRO_LSPLINE *lspline;
RL_INT4 k1, k2, kmin, kmax, k, lower_flag, upper_flag;
RL_FLT8 findex, x1, x2, y1, y2;
/* Get linear spline pointer */
lspline = (ZPRO_LSPLINE *) pointer;
/* Retrieve the domain of indexes for the function */
Pro_SeriesIndices(lspline->series, &kmin, &kmax);
/* Find index range that straddles x */
findex = Pro_SeriesIndex(lspline->series, x);
k1 = (RL_INT4) floor(findex);
k2 = k1 + 1;
/* Handle special case of evaluation at upper endpoint */
if (k2 > lspline->kmax) {k1--; k2--;}
/* Find first unflagged series value below x */
for (k = k1; k >= kmin; k--) {
y1 = Pro_SeriesValue(lspline->series, k, &lower_flag);
if (lower_flag == 0) break;
}
x1 = Pro_SeriesXValue(lspline->series, (RL_FLT8) k);
/* Find first unflagged series value above x */
for (k = k2; k <= kmax; k++) {
y2 = Pro_SeriesValue(lspline->series, k, &upper_flag);
if (upper_flag == 0) break;
}
x2 = Pro_SeriesXValue(lspline->series, (RL_FLT8) k);
/* Handle case of successful interpolation */
if (lower_flag == 0 && upper_flag == 0)
return y1 + (y2-y1) * (x-x1)/(x2-x1);
/* This situation should never occur since the end points of the domain are
* checked for validity when the object is created. However, somebody could
* conceivably muck around with the series values after the linear spline
* function has been created.
*/
(void) sprintf(xpro_message, "function evaluation failure\n\
missing adjacent value for linear spline interpolation\n\
values near index %#g (x = %#g) have been changed in series \"%s\"",
findex, x, Pro_ObjectName(lspline->series,0));
RL_RaiseError("PRO_EVALUATION_FAILURE", xpro_message);
/* Return best possible value */
if (lower_flag == 0) return y1;
if (upper_flag == 0) return y2;
return 0.;
}
/*
********************************************************************************
* ZPro_FreeLSpline(pointer)
*
* This internal routine frees the memory used by a ZPRO_LSPLINE structure.
*
* Inputs:
* pointer pointer to the ZPRO_LSPLINE data structure.
*******************************************************************************/
static void ZPro_FreeLSpline(pointer)
RL_VOID *pointer;
{
XRL_Free(pointer);
}
/*
********************************************************************************
* ZPro_PrintLSpline(pointer)
*
* This internal routine prints information about the ZPRO_LSPLINE structure.
*
* Input:
* pointer pointer to soft function object.
****************************************************************************/
static void ZPro_PrintLSpline(pointer)
RL_VOID *pointer;
{
ZPRO_LSPLINE *lspline;
lspline = (ZPRO_LSPLINE *) pointer;
/* Make sure object is not NULL */
if (lspline == NULL) {
printf("PRINT ERROR: Linear spline function pointer is NULL\n");
return;
}
/* Make sure object is a linear spline function */
if (lspline->class.id != XPRO_LSPLINE_CLASS) {
printf("PRINT ERROR: Object is not a linear spline function\n");
return;
}
/* Print object info... */
printf("\nLinear spline function parameters...\n");
printf("series = "); XPro_PrintInfo(lspline->series);
return;
}
/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
* FPro_LSplineFunc (lspline.c)
*$ Abstract:
* This routine generates a new linear spline function object, which
* performs linear interpolation between the nearest valid series values.
*$ Keywords:
* PROFILE, FUNCTION
* FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
* integer*4 function FPro_LSplineFunc(series)
* integer*4 series
*$ Inputs:
* series FORTRAN pointer to the series object.
*$ Outputs:
* none
*$ Returns:
* FORTRAN pointer to a new linear spline object, or 0 on non-fatal error.
*$ Detailed_description:
* This routine generates a new linear spline function object, which
* performs linear interpolation between the nearest valid series values.
*$ External_references:
* Profile toolkit
*$ Side_effects:
* Memory is allocated. The new object retains a link to the original
* series.
*$ Examples:
* This snippet of code takes a series that contains only three valid
* samples, and uses it to construct an linear spline function. The
* function returns f(x)=x for 1 <= x <= 5.
*
* integer*4 series, func
* real*8 table(5)/ 1.d0, 2.d0, 98.d0, 99.d0, 5.d0 /
*
* series = FPro_ArraySeries(table, 1, 5, 1.d0, 1.d0,
* & 98.d0, 99.d0, .TRUE., .FALSE.)
* func = FPro_LSplineFunc(series)
*$ Error_handling:
* Profile library error handling is in effect.
*
* Conditions raised:
* PRO_CLASS_ERROR if series object is NULL or not a series.
* PRO_SETUP_FAILURE if series contains less than two valid values.
* RL_MEMORY_ERROR on memory allocation error.
* FORTRAN_POINTER_ERROR if series is not a valid FORTRAN pointer.
*$ 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_lsplinefunc) (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_LSplineFunc((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;
}
/*******************************************************************************
*/