diff options
Diffstat (limited to 'src/jtutil/interpolate.hh')
| -rw-r--r-- | src/jtutil/interpolate.hh | 96 |
1 files changed, 67 insertions, 29 deletions
diff --git a/src/jtutil/interpolate.hh b/src/jtutil/interpolate.hh index c636fa5..715ac65 100644 --- a/src/jtutil/interpolate.hh +++ b/src/jtutil/interpolate.hh @@ -68,7 +68,7 @@ // interpolation position lies (fuzzily) outside the data range. // '*' ==> Like 'o', but also allow the interpolation position to lie // (fuzzily) outside the data range. In other words, this case -// case accepts *any* interpolation position. +// accepts *any* interpolation position. // // For Lagrange interpolation all three possibilities are valid; // for spline interpolation the "off-centering as needed" is intrinsic @@ -89,15 +89,11 @@ // 4. Set up the x coordinates where interpolated data is desired. // 5. Do the actual interpolation. // -// The routines here all allow 2-5, 3-5, and/or 4-5 to be repeated as needed. // The idea is that earlier steps can preallocate storage arrays and/or // precompute coefficients to be used by later stages, so repeating the // later stages may have lower overheads than redoing the entire process. -// -// One could also imagine doing these steps in the order 12435, i.e. -// interpolating many different sets of date to the same position(s). -// The Cactus interpolator interface supports this nicely, but alas -// the interface here doesn't support this at all. :( :( +// The interface here requires that 1 be done first, then allows 2-4 +// to be done in any order before 5. // // @@ -109,15 +105,17 @@ //****************************************************************************** +namespace interpolation + { + +//************************************** + // // This class provides 1D Lagrange global polynomial interpolation // for uniformly spaced data. // // ... const qualifiers refer to results of interpolate() // - -namespace interpolation - { template <typename fpt> class Lagrange_uniform_1d { @@ -127,12 +125,21 @@ public: void setup_data_x(const linear_map<fpt>& x_map_in); // this function sets up the y data values + // n.b. indexing off this pointer is [i], i.e. it need not be 0-origin void setup_data_y(const fpt y_array_in[]); + // this function sets up + // * the x coordinate where interpolation is to be done + // * what sort of off-the-end-of-the-data handling is desired + void setup_interp_x(fpt x_interp_in, char end_action_in = 'o'); + // the actual interpolation - // ... needs setup_data_x() and setup_data_y() - // to have already been called - fpt interpolate(fpt x, char end_action = 'o') const; + // ... requires that + // setup_data_x() + // setup_data_y() + // setup_interp_x() + // have already been called + fpt interpolate(); // this function computes the backwards domain of dependence // of the interpolation, i.e. the interval of [i] such that @@ -141,20 +148,24 @@ public: // for interpolate() // // ... returns domain-of-dependence interval in [i_lo,i_hi] - // ... returns nominal center point (i_center in source code) - // as function value - // ... needs setup_data_x() to have already been called, - // but does *not* need setup_data_y() to have already been called - int bdod(fpt x, char end_action = 'o', int& i_lo, int& i_hi) const; - - // Jacobian d interpolate(x,end_action) / d y[i] - // ... needs setup_data_x() to have already been called, - // but does *not* need setup_data_y() to have already been called - fpt Jacobian(fpt x, char end_action = 'o', int i) const; + // ... returns nominal center point of underlying local interpolant + // (i_center in source code) as function value + // ... requires that + // setup_data_x() + // setup_interp_x() + // have already been called + int bdod(int& i_lo, int& i_hi) const; + + // Jacobian d (value returned by interpolate()) / d y[i] + // ... requires that + // setup_data_x() + // setup_interp_x() + // have already been called + fpt Jacobian(int i) const; // constructor, destructor + // ... compiler-generated no-op destructor is ok Lagrange_uniform_1d(int order_in); - ~Lagrange_uniform_1d() { } private: // we forbid copying and passing by value @@ -166,13 +177,40 @@ private: private: const int order; - const int N_interp; // number of points in local - // interpolator, i.e. - // linear = 2 points, - // quadratic = 3 points, etc + const int N_interp; // number of points in local + // interpolator, i.e. + // linear = 2 points, + // quadratic = 3 points, etc + + bool x_interp_valid; // valid flag for x_interp, end_action + fpt x_interp; + char end_action; // n.b. non-const pointers to const data! const linear_map<fpt> *x_map_ptr; // --> client-owned data const fpt *y_array; // --> client-owned data }; - }; + +//************************************** + +// +// ***** underlying local interpolants of varying orders ***** +// + +// ... x here is in units of the grid spacing, +// relative to the center point +// ... these make no checks for data being in range +template <typename fpt> + fpt interpolate_local_order1(const fpt x, const fpt *y_center); +template <typename fpt> + fpt interpolate_local_order2(const fpt x, const fpt *y_center); +template <typename fpt> + fpt interpolate_local_order3(const fpt x, const fpt *y_center); +template <typename fpt> + fpt interpolate_local_order4(const fpt x, const fpt *y_center); +fpt interpolate_local_order5(const fpt x, const fpt *y_center); +fpt interpolate_local_order6(const fpt x, const fpt *y_center); + +//************************************** + + }; // close namespace interpolation:: |