Hex  1.0 Hydrogen-electron collision solver
interpolate.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2  * *
3  * / / / / __ \ \ / / *
4  * / /__ / / / _ \ \ \/ / *
5  * / ___ / | |/_/ / /\ \ *
6  * / / / / \_\ / / \ \ *
7  * *
8  * Jakub Benda (c) 2014 *
9  * Charles University in Prague *
10  * *
11 \* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12
13 #ifndef HEX_INTERPOLATE
14 #define HEX_INTERPOLATE
15
16 #include <algorithm>
17
18 #include <gsl/gsl_interp.h>
19
20 #include "arrays.h"
21
30 template <typename T> NumberArray<T> interpolate (rArray const & x0, NumberArray<T> const & y0, rArray const & x)
31 {
32 // if (x0.size() == 0)
33 // throw exception ("Nothing to interpolate.\n");
34
35  if (x0.size() < 2)
36  return y0;
37
38  // output array
39  NumberArray<T> y(x.size());
40
41  for (size_t i = 0; i < x.size(); i++)
42  {
43  // right guardian
44  rArray::const_iterator right = std::upper_bound(x0.begin(), x0.end(), x[i]);
45
46  if (right == x0.end() or right == x0.begin())
47  {
48  y[i] = T(0);
49  continue;
50  }
51
52  // neighbours
53  double x0_left = *(right - 1);
54  double x0_right = *right;
55  T y0_left = y0[right - 1 - x0.begin()];
56  T y0_right = y0[right - x0.begin()];
57
58  // compute linear average
59  y[i] = (y0_left * (x0_right - x[i]) + y0_right * (x[i] - x0_left)) / (x0_right - x0_left);
60  }
61
62  // FIXME implement better schemes
63
64  return y;
65 }
66
97 inline rArray interpolate_real (rArray const & x0, rArray const & y0, rArray const & x, const gsl_interp_type * interpolation)
98 {
99  // return 0 if the array to interpolate is empty
100  if (x0.size() == 0)
101  return rArray(x.size(), 0.);
102
103  // return the source element if the array to interpolate contains just the one element
104  if (x0.size() == 1)
105  return rArray(x.size(), y0[0]);
106
107  // setup the interpolator
108  gsl_interp *itp = gsl_interp_alloc (interpolation, x0.size());
109  gsl_interp_init (itp, x0.data(), y0.data(), x0.size());
110  gsl_interp_accel *accel = gsl_interp_accel_alloc ();
111
112  // interpolate
113  rArray y(x.size());
114  for (size_t i = 0; i < x.size(); i++)
115  {
116  // check that we are not extrapolating
117  if (x0.front() <= x[i] and x[i] <= x0.back())
118  y[i] = gsl_interp_eval (itp, x0.data(), y0.data(), x[i], accel);
119  else
120  y[i] = 0.;
121  }
122
123  // release memory
124  gsl_interp_accel_free (accel);
125  gsl_interp_free (itp);
126
127  return y;
128 }
129
130 #endif
virtual T * data()
Data pointer.
Definition: arrays.h:757
T const * const_iterator
Definition: arrays.h:611
T & front(int i=0)
Definition: arrays.h:777
A comfortable number array class.
Definition: arrays.h:171
T & back(int i=0)
Definition: arrays.h:779
iterator end()
Definition: arrays.h:774
iterator begin()
Definition: arrays.h:771
NumberArray< double > rArray
Definition: arrays.h:1609
rArray interpolate_real(rArray const &x0, rArray const &y0, rArray const &x, const gsl_interp_type *interpolation)
Return values interpolated by O₂scl.
Definition: interpolate.h:97
NumberArray< T > interpolate(rArray const &x0, NumberArray< T > const &y0, rArray const &x)
Return linearly interpolated values.
Definition: interpolate.h:30
size_t size() const
Item count.
Definition: arrays.h:673