SDDS ToolKit Programs and Libraries for C and Python
Loading...
Searching...
No Matches
Functions
interp.c File Reference

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.

#include "mdb.h"

Go to the source code of this file.

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.
 

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
fPointer to the array of function values.
xPointer to the array of independent variable values.
nNumber of data points.
xoThe point at which to interpolate.
warningsFlag to enable or disable warning messages.
orderThe order of interpolation.
returnCodePointer to a variable to store the return code status.
Returns
The interpolated value at point xo.

Definition at line 34 of file interp.c.

34 {
35 long hi, lo, mid, offset;
36
37 lo = 0;
38 hi = n - 1;
39 if (lo == hi) {
40 if (warnings)
41 printf("warning: only one point--returning value for that point\n");
42 *returnCode = 0;
43 return (f[0]);
44 }
45 if (x[lo] < x[hi]) { /* indep variable ordered from small to large */
46 if (xo < x[lo = 0]) {
47 if (warnings)
48 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
49 xo, x[0], x[n - 1]);
50 *returnCode = 0;
51 return (f[lo]);
52 }
53 if (xo > x[hi = n - 1]) {
54 if (warnings)
55 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
56 xo, x[0], x[n - 1]);
57 *returnCode = 0;
58 return (f[hi]);
59 }
60
61 /* do binary search for closest point */
62 while ((hi - lo) > 1) {
63 mid = (lo + hi) / 2;
64 if (xo < x[mid])
65 hi = mid;
66 else
67 lo = mid;
68 }
69 } else { /* indep variable ordered from large to small */
70 if (xo > x[lo = 0]) {
71 if (warnings)
72 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
73 xo, x[n - 1], x[0]);
74 *returnCode = 0;
75 return (f[lo]);
76 }
77 if (xo < x[hi = n - 1]) {
78 if (warnings)
79 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
80 xo, x[n - 1], x[0]);
81 *returnCode = 0;
82 return (f[hi]);
83 }
84
85 /* do binary search for closest point */
86 while ((hi - lo) > 1) {
87 mid = (lo + hi) / 2;
88 if (xo > x[mid])
89 hi = mid;
90 else
91 lo = mid;
92 }
93 }
94 /* L.Emery's contribution */
95 offset = lo - (order - 1) / 2; /* offset centers the argument in the set of points. */
96 offset = MAX(offset, 0);
97 offset = MIN(offset, n - order - 1);
98 offset = MAX(offset, 0);
99 return LagrangeInterp(x + offset, f + offset, order + 1, xo, returnCode);
100}
LagrangeInterp
double LagrangeInterp(double *x, double *f, long order1, double x0, long *returnCode)
Performs Lagrange interpolation of data.
Definition interp.c:115

◆ 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
fPointer to the array of short integer function values.
xPointer to the array of independent variable values.
nNumber of data points.
xoThe point at which to interpolate.
warningsFlag to enable or disable warning messages.
orderThe order of interpolation.
returnCodePointer to a variable to store the return code status.
next_start_posPointer 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.

282 {
283 int64_t hi, lo, mid, i;
284 short value;
285
286 lo = 0;
287 hi = n - 1;
288
289 *returnCode = 0;
290 /*if the interpolate point is less or equal to the start point, return the value of the start point */
291 if (xo <= x[0]) {
292 *next_start_pos = 0;
293 return f[0];
294 }
295
296 /*if the value is in one of the indepent, return the corresponding f value */
297 for (i = 0; i < n; i++)
298 if (xo == x[i]) {
299 *next_start_pos = i;
300 return f[i];
301 }
302 if (lo == hi) {
303 if (warnings)
304 printf("warning: only one point--returning value for that point\n");
305 *returnCode = 0;
306 *next_start_pos = lo;
307 return (f[0]);
308 }
309 if (x[lo] < x[hi]) { /* indep variable ordered from small to large */
310 if (xo < x[lo = 0]) {
311 if (warnings)
312 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
313 xo, x[0], x[n - 1]);
314 *returnCode = 0;
315 *next_start_pos = lo;
316 return (f[lo]);
317 }
318 if (xo > x[hi = n - 1]) {
319 if (warnings)
320 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
321 xo, x[0], x[n - 1]);
322 *returnCode = 0;
323 *next_start_pos = hi;
324 return (f[hi]);
325 }
326
327 /* do binary search for closest point */
328 while ((hi - lo) > 1) {
329 mid = (lo + hi) / 2;
330 if (xo < x[mid])
331 hi = mid;
332 else
333 lo = mid;
334 }
335 } else { /* indep variable ordered from large to small */
336 if (xo > x[lo = 0]) {
337 if (warnings)
338 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
339 xo, x[n - 1], x[0]);
340 *returnCode = 0;
341 *next_start_pos = lo;
342 return (f[lo]);
343 }
344 if (xo < x[hi = n - 1]) {
345 if (warnings)
346 printf("warning: %22.15e outside [%22.15e,%22.15e] (interp)\n",
347 xo, x[n - 1], x[0]);
348 *returnCode = 0;
349 *next_start_pos = hi;
350 return (f[hi]);
351 }
352
353 /* do binary search for closest point */
354 while ((hi - lo) > 1) {
355 mid = (lo + hi) / 2;
356 if (xo > x[mid])
357 hi = mid;
358 else
359 lo = mid;
360 }
361 }
362 /*remember this position so that one can start from her to speed up the interp of
363 other points assume the interping points are sorted.*/
364 *next_start_pos = lo;
365 *returnCode = 0;
366 if (order == -1) {
367 /* inherit value from previous point*/
368 value = (short)f[lo];
369 } else if (order == -2) {
370 /* inherit value from next point */
371 value = (short)f[hi];
372 } else {
373 value = (f[hi] - f[lo]) / (xo - x[lo]);
374 }
375 return value;
376}

◆ 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
fPointer to the array of function values.
xPointer to the array of independent variable values.
nNumber of data points.
xoThe point at which to interpolate.
belowRangePointer to structure controlling behavior below the data range.
aboveRangePointer to structure controlling behavior above the data range.
orderThe order of interpolation.
returnCodePointer to a variable to store the return code status.
MMultiplier to adjust the interpolation condition based on the order.
Returns
The interpolated value at point xo.

Definition at line 160 of file interp.c.

161 {
162 long code;
163 int64_t hi, lo, mid, offset;
164 double result, below, above;
165
166 *returnCode = 0;
167
168 lo = 0;
169 hi = n - 1;
170 if (M > 0) {
171 above = f[n - 1];
172 below = f[0];
173 } else {
174 above = f[0];
175 below = f[n - 1];
176 }
177 if ((xo * M > x[hi] * M && M > 0) || (xo * M < x[lo] * M && M < 0)) {
178 if (aboveRange->flags & OUTRANGE_SKIP) {
179 *returnCode = OUTRANGE_SKIP;
180 return above;
181 } else if (aboveRange->flags & OUTRANGE_ABORT) {
182 *returnCode = OUTRANGE_ABORT;
183 return above;
184 } else if (aboveRange->flags & OUTRANGE_WARN)
185 *returnCode = OUTRANGE_WARN;
186 if (aboveRange->flags & OUTRANGE_VALUE) {
187 *returnCode |= OUTRANGE_VALUE;
188 return aboveRange->value;
189 }
190 if (aboveRange->flags & OUTRANGE_WRAP) {
191 double delta;
192 *returnCode |= OUTRANGE_WRAP;
193 if ((delta = x[hi] - x[lo]) == 0)
194 return f[0];
195 while (xo * M > x[hi] * M)
196 xo -= delta;
197 } else if (aboveRange->flags & OUTRANGE_SATURATE || !(aboveRange->flags & OUTRANGE_EXTRAPOLATE)) {
198 *returnCode |= OUTRANGE_SATURATE;
199 return above;
200 }
201 }
202 if ((xo * M < x[lo] * M && M > 0) || (xo * M > x[hi] * M && M < 0)) {
203 if (belowRange->flags & OUTRANGE_SKIP) {
204 *returnCode = OUTRANGE_SKIP;
205 return below;
206 } else if (belowRange->flags & OUTRANGE_ABORT) {
207 *returnCode = OUTRANGE_ABORT;
208 return below;
209 } else if (belowRange->flags & OUTRANGE_WARN)
210 *returnCode = OUTRANGE_WARN;
211 if (belowRange->flags & OUTRANGE_VALUE) {
212 *returnCode |= OUTRANGE_VALUE;
213 return belowRange->value;
214 }
215 if (belowRange->flags & OUTRANGE_WRAP) {
216 double delta;
217 *returnCode |= OUTRANGE_WRAP;
218 if ((delta = x[hi] - x[lo]) == 0)
219 return below;
220 while (xo * M < x[lo] * M)
221 xo += delta;
222 } else if (belowRange->flags & OUTRANGE_SATURATE || !(belowRange->flags & OUTRANGE_EXTRAPOLATE)) {
223 *returnCode |= OUTRANGE_SATURATE;
224 return below;
225 }
226 }
227
228 if (lo == hi) {
229 if (xo == x[lo]) {
230 if (aboveRange->flags & OUTRANGE_WARN || belowRange->flags & OUTRANGE_WARN)
231 *returnCode = OUTRANGE_WARN;
232 }
233 return f[0];
234 }
235
236 lo = 0;
237 hi = n - 1;
238 if (xo * M < x[0] * M)
239 hi = 1;
240 else if (xo * M > x[n - 1] * M)
241 lo = hi - 1;
242 else {
243 /* do binary search for closest point */
244 while ((hi - lo) > 1) {
245 mid = (lo + hi) / 2;
246 if (xo * M < x[mid] * M)
247 hi = mid;
248 else
249 lo = mid;
250 }
251 }
252
253 /* L.Emery's contribution */
254 if (order > n - 1)
255 order = n - 1;
256 offset = lo - (order - 1) / 2; /* offset centers the argument in the set of points. */
257 offset = MAX(offset, 0);
258 offset = MIN(offset, n - order - 1);
259 result = LagrangeInterp(x + offset, f + offset, order + 1, xo, &code);
260 if (!code)
261 bomb("zero denominator in LagrangeInterp", NULL);
262 return result;
263}
bomb
void bomb(char *error, char *usage)
Reports error messages to the terminal and aborts the program.
Definition bomb.c:26

◆ 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
xPointer to the array of independent variable values.
fPointer to the array of function values.
order1The order of the Lagrange polynomial.
x0The point at which to interpolate.
returnCodePointer to a variable to store the return code status.
Returns
The interpolated value at point x0.

Definition at line 115 of file interp.c.

115 {
116 long i, j;
117 double denom, numer, sum;
118
119 for (i = sum = 0; i < order1; i++) {
120 denom = 1;
121 numer = 1;
122 for (j = 0; j < order1; j++) {
123 if (i != j) {
124 denom *= (x[i] - x[j]);
125 numer *= (x0 - x[j]);
126 if (numer == 0) {
127 *returnCode = 1;
128 return f[j];
129 }
130 }
131 }
132 if (denom == 0) {
133 *returnCode = 0;
134 return 0.0;
135 }
136 sum += f[i] * numer / denom;
137 }
138 *returnCode = 1;
139 return sum;
140}
Generated by doxygen 1.12.0