Implements interpolation functions including linear and Lagrange interpolation. More...

#include "mdb.h"

Functions
double	interp (double f, double x, long n, double xo, long warnings, long order, long *returnCode)
	Performs simple linear interpolation of data.

double	LagrangeInterp (double x, double f, long order1, double x0, long *returnCode)
	Performs Lagrange interpolation of data.

double	interpolate (double f, double x, int64_t n, double xo, OUTRANGE_CONTROL belowRange, OUTRANGE_CONTROL aboveRange, long order, unsigned long *returnCode, long M)
	Performs interpolation with range control options.

short	interp_short (short f, double x, int64_t n, double xo, long warnings, short order, unsigned long returnCode, long next_start_pos)
	Performs interpolation for short integer data types.

Detailed Description

Implements interpolation functions including linear and Lagrange interpolation.

Copyright

(c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
(c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License: This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author: M. Borland, C. Saunders, R. Soliday, H. Shang, N. Kuklev

Definition in file interp.c.

Function Documentation

◆ interp()

double interp	(	double *	f,
		double *	x,
		long	n,
		double	xo,
		long	warnings,
		long	order,
		long *	returnCode )

Performs simple linear interpolation of data.

This function interpolates the value at a given point using linear interpolation. It handles boundary conditions by returning the nearest data point value when the interpolation point is outside the data range.

Parameters

f	Pointer to the array of function values.
x	Pointer to the array of independent variable values.
n	Number of data points.
xo	The point at which to interpolate.
warnings	Flag to enable or disable warning messages.
order	The order of interpolation.
returnCode	Pointer to a variable to store the return code status.

Returns: The interpolated value at point xo.

Definition at line 34 of file interp.c.

                                                                                                    {
  long hi, lo, mid, offset;
 
  lo = 0;
  hi = n - 1;
  if (lo == hi) {
    if (warnings)
      printf("warning: only one point--returning value for that point\n");
    *returnCode = 0;
    return (f[0]);
  }
  if (x[lo] < x[hi]) { /* indep variable ordered from small to large */
    if (xo < x[lo = 0]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[0], x[n - 1]);
      *returnCode = 0;
      return (f[lo]);
    }
    if (xo > x[hi = n - 1]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[0], x[n - 1]);
      *returnCode = 0;
      return (f[hi]);
    }
 
    /* do binary search for closest point */
    while ((hi - lo) > 1) {
      mid = (lo + hi) / 2;
      if (xo < x[mid])
        hi = mid;
      else
        lo = mid;
    }
  } else { /* indep variable ordered from large to small */
    if (xo > x[lo = 0]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[n - 1], x[0]);
      *returnCode = 0;
      return (f[lo]);
    }
    if (xo < x[hi = n - 1]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[n - 1], x[0]);
      *returnCode = 0;
      return (f[hi]);
    }
 
    /* do binary search for closest point */
    while ((hi - lo) > 1) {
      mid = (lo + hi) / 2;
      if (xo > x[mid])
        hi = mid;
      else
        lo = mid;
    }
  }
  /* L.Emery's contribution */
  offset = lo - (order - 1) / 2; /* offset centers the argument in the set of points. */
  offset = MAX(offset, 0);
  offset = MIN(offset, n - order - 1);
  offset = MAX(offset, 0);
  return LagrangeInterp(x + offset, f + offset, order + 1, xo, returnCode);
}

◆ interp_short()

short interp_short	(	short *	f,
		double *	x,
		int64_t	n,
		double	xo,
		long	warnings,
		short	order,
		unsigned long *	returnCode,
		long *	next_start_pos )

Performs interpolation for short integer data types.

This function interpolates the value at a given point for short integer data, handling boundary conditions and allowing for different interpolation orders.

Parameters

f	Pointer to the array of short integer function values.
x	Pointer to the array of independent variable values.
n	Number of data points.
xo	The point at which to interpolate.
warnings	Flag to enable or disable warning messages.
order	The order of interpolation.
returnCode	Pointer to a variable to store the return code status.
next_start_pos	Pointer to store the next starting position for interpolation.

Returns: The interpolated short integer value at point xo.

Definition at line 281 of file interp.c.

                                                                    {
  int64_t hi, lo, mid, i;
  short value;
 
  lo = 0;
  hi = n - 1;
 
  *returnCode = 0;
  /*if the interpolate point is less or equal to the start point, return the value of the start point */
  if (xo <= x[0]) {
    *next_start_pos = 0;
    return f[0];
  }
 
  /*if the value is in one of the indepent, return the corresponding f value */
  for (i = 0; i < n; i++)
    if (xo == x[i]) {
      *next_start_pos = i;
      return f[i];
    }
  if (lo == hi) {
    if (warnings)
      printf("warning: only one point--returning value for that point\n");
    *returnCode = 0;
    *next_start_pos = lo;
    return (f[0]);
  }
  if (x[lo] < x[hi]) { /* indep variable ordered from small to large */
    if (xo < x[lo = 0]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[0], x[n - 1]);
      *returnCode = 0;
      *next_start_pos = lo;
      return (f[lo]);
    }
    if (xo > x[hi = n - 1]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[0], x[n - 1]);
      *returnCode = 0;
      *next_start_pos = hi;
      return (f[hi]);
    }
 
    /* do binary search for closest point */
    while ((hi - lo) > 1) {
      mid = (lo + hi) / 2;
      if (xo < x[mid])
        hi = mid;
      else
        lo = mid;
    }
  } else { /* indep variable ordered from large to small */
    if (xo > x[lo = 0]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[n - 1], x[0]);
      *returnCode = 0;
      *next_start_pos = lo;
      return (f[lo]);
    }
    if (xo < x[hi = n - 1]) {
      if (warnings)
        printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
               xo, x[n - 1], x[0]);
      *returnCode = 0;
      *next_start_pos = hi;
      return (f[hi]);
    }
 
    /* do binary search for closest point */
    while ((hi - lo) > 1) {
      mid = (lo + hi) / 2;
      if (xo > x[mid])
        hi = mid;
      else
        lo = mid;
    }
  }
  /*remember this position so that one can start from her to speed up the interp of
    other points assume the interping points are sorted.*/
  *next_start_pos = lo;
  *returnCode = 0;
  if (order == -1) {
    /* inherit value from previous point*/
    value = (short)f[lo];
  } else if (order == -2) {
    /* inherit value from next point */
    value = (short)f[hi];
  } else {
    value = (f[hi] - f[lo]) / (xo - x[lo]);
  }
  return value;
}

◆ interpolate()

double interpolate	(	double *	f,
		double *	x,
		int64_t	n,
		double	xo,
		OUTRANGE_CONTROL *	belowRange,
		OUTRANGE_CONTROL *	aboveRange,
		long	order,
		unsigned long *	returnCode,
		long	M )

Performs interpolation with range control options.

This function interpolates the value at a given point with additional control over out-of-range conditions, such as skipping, aborting, warning, wrapping, saturating, or using a specified value.

Parameters

f	Pointer to the array of function values.
x	Pointer to the array of independent variable values.
n	Number of data points.
xo	The point at which to interpolate.
belowRange	Pointer to structure controlling behavior below the data range.
aboveRange	Pointer to structure controlling behavior above the data range.
order	The order of interpolation.
returnCode	Pointer to a variable to store the return code status.
M	Multiplier to adjust the interpolation condition based on the order.

Returns: The interpolated value at point xo.

Definition at line 160 of file interp.c.

                                                                                                {
  long code;
  int64_t hi, lo, mid, offset;
  double result, below, above;
 
  *returnCode = 0;
 
  lo = 0;
  hi = n - 1;
  if (M > 0) {
    above = f[n - 1];
    below = f[0];
  } else {
    above = f[0];
    below = f[n - 1];
  }
  if ((xo * M > x[hi] * M && M > 0) || (xo * M < x[lo] * M && M < 0)) {
    if (aboveRange->flags & OUTRANGE_SKIP) {
      *returnCode = OUTRANGE_SKIP;
      return above;
    } else if (aboveRange->flags & OUTRANGE_ABORT) {
      *returnCode = OUTRANGE_ABORT;
      return above;
    } else if (aboveRange->flags & OUTRANGE_WARN)
      *returnCode = OUTRANGE_WARN;
    if (aboveRange->flags & OUTRANGE_VALUE) {
      *returnCode |= OUTRANGE_VALUE;
      return aboveRange->value;
    }
    if (aboveRange->flags & OUTRANGE_WRAP) {
      double delta;
      *returnCode |= OUTRANGE_WRAP;
      if ((delta = x[hi] - x[lo]) == 0)
        return f[0];
      while (xo * M > x[hi] * M)
        xo -= delta;
    } else if (aboveRange->flags & OUTRANGE_SATURATE || !(aboveRange->flags & OUTRANGE_EXTRAPOLATE)) {
      *returnCode |= OUTRANGE_SATURATE;
      return above;
    }
  }
  if ((xo * M < x[lo] * M && M > 0) || (xo * M > x[hi] * M && M < 0)) {
    if (belowRange->flags & OUTRANGE_SKIP) {
      *returnCode = OUTRANGE_SKIP;
      return below;
    } else if (belowRange->flags & OUTRANGE_ABORT) {
      *returnCode = OUTRANGE_ABORT;
      return below;
    } else if (belowRange->flags & OUTRANGE_WARN)
      *returnCode = OUTRANGE_WARN;
    if (belowRange->flags & OUTRANGE_VALUE) {
      *returnCode |= OUTRANGE_VALUE;
      return belowRange->value;
    }
    if (belowRange->flags & OUTRANGE_WRAP) {
      double delta;
      *returnCode |= OUTRANGE_WRAP;
      if ((delta = x[hi] - x[lo]) == 0)
        return below;
      while (xo * M < x[lo] * M)
        xo += delta;
    } else if (belowRange->flags & OUTRANGE_SATURATE || !(belowRange->flags & OUTRANGE_EXTRAPOLATE)) {
      *returnCode |= OUTRANGE_SATURATE;
      return below;
    }
  }
 
  if (lo == hi) {
    if (xo == x[lo]) {
      if (aboveRange->flags & OUTRANGE_WARN || belowRange->flags & OUTRANGE_WARN)
        *returnCode = OUTRANGE_WARN;
    }
    return f[0];
  }
 
  lo = 0;
  hi = n - 1;
  if (xo * M < x[0] * M)
    hi = 1;
  else if (xo * M > x[n - 1] * M)
    lo = hi - 1;
  else {
    /* do binary search for closest point */
    while ((hi - lo) > 1) {
      mid = (lo + hi) / 2;
      if (xo * M < x[mid] * M)
        hi = mid;
      else
        lo = mid;
    }
  }
 
  /* L.Emery's contribution */
  if (order > n - 1)
    order = n - 1;
  offset = lo - (order - 1) / 2; /* offset centers the argument in the set of points. */
  offset = MAX(offset, 0);
  offset = MIN(offset, n - order - 1);
  result = LagrangeInterp(x + offset, f + offset, order + 1, xo, &code);
  if (!code)
    bomb("zero denominator in LagrangeInterp", NULL);
  return result;
}

◆ LagrangeInterp()

double LagrangeInterp	(	double *	x,
		double *	f,
		long	order1,
		double	x0,
		long *	returnCode )

Performs Lagrange interpolation of data.

This function computes the interpolated value at a given point using Lagrange polynomial interpolation of a specified order.

Parameters

x	Pointer to the array of independent variable values.
f	Pointer to the array of function values.
order1	The order of the Lagrange polynomial.
x0	The point at which to interpolate.
returnCode	Pointer to a variable to store the return code status.

Returns: The interpolated value at point x0.

Definition at line 115 of file interp.c.

                                                                                      {
  long i, j;
  double denom, numer, sum;
 
  for (i = sum = 0; i < order1; i++) {
    denom = 1;
    numer = 1;
    for (j = 0; j < order1; j++) {
      if (i != j) {
        denom *= (x[i] - x[j]);
        numer *= (x0 - x[j]);
        if (numer == 0) {
          *returnCode = 1;
          return f[j];
        }
      }
    }
    if (denom == 0) {
      *returnCode = 0;
      return 0.0;
    }
    sum += f[i] * numer / denom;
  }
  *returnCode = 1;
  return sum;
}

Functions

Detailed Description

Function Documentation

◆ interp()

◆ interp_short()

◆ interpolate()

◆ LagrangeInterp()