/*
******************************************************************************** * 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 Rings Node, April 1998 *******************************************************************************/ #include*/#include #include #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; } /*******************************************************************************