80 lines
1.8 KiB
C++
80 lines
1.8 KiB
C++
#pragma once
|
|
|
|
#include "./numerics/interpolation1d/interpolation1d_base.h"
|
|
|
|
#include "./utils/vector.h"
|
|
#include "./numerics/abs.h"
|
|
|
|
namespace numerics{
|
|
|
|
template <typename T>
|
|
struct interp_rational : Base_interp<T> {
|
|
using Base = Base_interp<T>;
|
|
|
|
// bring base data members into scope (or use this->xx / this->yy below)
|
|
using Base::xx;
|
|
using Base::yy;
|
|
using Base::mm;
|
|
|
|
T dy;
|
|
|
|
|
|
|
|
interp_rational(const utils::Vector<T> &xv, const utils::Vector<T> &yv, uint64_t m)
|
|
: Base_interp<T>(xv, &yv[0], m), dy(T{0}){}
|
|
|
|
T rawinterp(int64_t jl, T x) override{
|
|
|
|
const T TINY = T{1.0e-99};
|
|
int64_t ns=0;
|
|
T y, w, t, hh, h, dd;
|
|
const T *xa = &xx[jl], *ya = &yy[jl];
|
|
utils::Vector<T> c(mm, T{0}), d(mm, T{0});
|
|
|
|
hh = numerics::abs(x - xa[0]);
|
|
|
|
for (int64_t i = 0; i < mm; ++i){
|
|
h = numerics::abs(x-xa[i]);
|
|
if (h == T{0}){
|
|
dy = T{0};
|
|
return ya[i];
|
|
}else if (h < hh){
|
|
ns = i;
|
|
hh = h;
|
|
}
|
|
c[i] = ya[i];
|
|
d[i] = ya[i] + TINY; // The TINY part is needed to prevent a rare zero-over-zero condition.
|
|
}
|
|
y = ya[ns];
|
|
ns -= 1;
|
|
for (int64_t m = 1; m < mm; ++m){
|
|
for (int64_t i = 0; i < mm-m; ++i){
|
|
w = c[i+1] - d[i];
|
|
h = xa[i+m] - x; // h will never be zero, since this was tested in the initializing loop.
|
|
t = (xa[i] - x)*d[i]/h;
|
|
dd = t - c[i+1];
|
|
if (dd == T{0}){ // This error condition indicates that the interpolating function has a pole at the requested value of x.
|
|
throw std::invalid_argument("Error in routine interp_rational"); //
|
|
}
|
|
dd = w/dd;
|
|
d[i] = c[i+1]*dd;
|
|
c[i] = t*dd;
|
|
}
|
|
const bool take_left = (2 * (ns + 1) < (mm - m));
|
|
|
|
if (take_left) {
|
|
dy = c[ns + 1];
|
|
} else {
|
|
dy = d[ns];
|
|
ns -= 1;
|
|
}
|
|
y += dy;
|
|
}
|
|
return y;
|
|
}
|
|
|
|
};
|
|
|
|
} // namespace numerics
|
|
|