Canoe
Comprehensive Atmosphere N' Ocean Engine
interpn.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 
3 #include "interpolation.h"
4 
12 void interpn(double *val, double const *coor, double const *data,
13  double const *axis, size_t const *len, int ndim, int nval) {
14  int i1, i2;
15  i1 = locate(axis, *coor, *len);
16 
17  // if the interpolation value is out of bound
18  // use the closest value
19  if (i1 == -1) {
20  i1 = 0;
21  i2 = 0;
22  } else if (i1 == *len - 1) {
23  i1 = *len - 1;
24  i2 = *len - 1;
25  } else
26  i2 = i1 + 1;
27 
28  double *v1 = (double *)malloc(nval * sizeof(double));
29  double *v2 = (double *)malloc(nval * sizeof(double));
30 
31  double x1 = axis[i1];
32  double x2 = axis[i2];
33 
34  if (ndim == 1) {
35  for (int j = 0; j < nval; ++j) {
36  v1[j] = data[i1 * nval + j];
37  v2[j] = data[i2 * nval + j];
38  }
39  } else {
40  int s = nval;
41  for (int j = 1; j < ndim; ++j) s *= len[j];
42  interpn(v1, coor + 1, data + i1 * s, axis + *len, len + 1, ndim - 1, nval);
43  interpn(v2, coor + 1, data + i2 * s, axis + *len, len + 1, ndim - 1, nval);
44  }
45 
46  if (x2 != x1)
47  for (int j = 0; j < nval; ++j)
48  val[j] = ((*coor - x1) * v2[j] + (x2 - *coor) * v1[j]) / (x2 - x1);
49  else
50  for (int j = 0; j < nval; ++j) val[j] = (v1[j] + v2[j]) / 2.;
51 
52  free(v1);
53  free(v2);
54 }
55 
61 double interp1(double x, double const *data, double const *axis, size_t len) {
62  double value;
63  interpn(&value, &x, data, axis, &len, 1, 1);
64  return value;
65 }
double interp1(double x, double const *data, double const *axis, size_t len)
Definition: interpn.c:61
void interpn(double *val, double const *coor, double const *data, double const *axis, size_t const *len, int ndim, int nval)
Definition: interpn.c:12
int locate(double const *xx, double x, int n)
Definition: locate.c:7