87 lines
2.6 KiB
C++
87 lines
2.6 KiB
C++
#pragma once
|
||
|
||
#include "./numerics/interpolation1d/interpolation1d_base.h"
|
||
|
||
//#include "./numerics/abs.h"
|
||
#include "./utils/vector.h"
|
||
|
||
|
||
namespace numerics{
|
||
|
||
template <typename T>
|
||
struct interp_cubic_spline : 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;
|
||
|
||
utils::Vector<T> y2;
|
||
|
||
|
||
interp_cubic_spline(utils::Vector<T> &xv, utils::Vector<T> &yv, T yp1=T{1.e99}, T ypn=T{1.e99})
|
||
: Base_interp<T>(xv, &yv[0], 2), y2(xv.size(),T{0}) {
|
||
sety2(&xv[0], &yv[0], yp1, ypn);
|
||
}
|
||
|
||
interp_cubic_spline(utils::Vector<T> &xv, const T *yv, T yp1=T{1.e99}, T ypn=T{1.e99})
|
||
: Base_interp<T>(xv, yv, 2), y2(xv.size(),T{0}) {
|
||
sety2(&xv[0], yv, yp1, ypn);
|
||
}
|
||
|
||
|
||
void sety2(const T *xv, const T *yv, T yp1, T ypn){
|
||
|
||
T p, qn, sig, un;
|
||
uint64_t n = y2.size();
|
||
utils::Vector<T> u(n-1, T{0});
|
||
|
||
if (yp1 > static_cast<T>(0.99e99)){ // The lower boundary condition is set either to be “natural”
|
||
y2[0] = u[0] = T{0};
|
||
}else{ // or else to have a specified first derivative.
|
||
y2[0] = T{-0.5};
|
||
u[0] = (3.0/(xv[1]-xv[0]))*(((yv[1]-yv[0])/(xv[1]-xv[0]))-yp1);
|
||
}
|
||
for (uint64_t i = 1; i < n-1; ++i){ // This is the decomposition loop of the tridiagonal algorithm
|
||
sig = (xv[i]-xv[i-1])/(xv[i+1]-xv[i-1]);
|
||
p = sig*y2[i-1]+T{2};
|
||
y2[i] = (sig - T{1})/p; // y2 and u are used for temporary storage of the decomposed factors.
|
||
u[i]=((yv[i+1]-yv[i])/(xv[i+1]-xv[i])) - ((yv[i]-yv[i-1])/(xv[i]-xv[i-1]));
|
||
u[i]=((T{6}*u[i]/(xv[i+1]-xv[i-1])) - sig*u[i-1])/p;
|
||
}
|
||
if (ypn > static_cast<T>(0.99e99)){ // The upper boundary condition is set either to be “natural”
|
||
qn = un = T{0};
|
||
}else{ // or else to have a specified first derivative.
|
||
qn = T{0.5};
|
||
un = (T{3}/(xv[n-1]-xv[n-2]))*(ypn-((yv[n-1]-yv[n-2])/(xv[n-1]-xv[n-2])));
|
||
}
|
||
y2[n-1] = (un-(qn*u[n-2]))/((qn*y2[n-2])+T{1});
|
||
for (int64_t k = n-2; k >= 0; --k){
|
||
y2[k] = y2[k] * y2[k+1]+u[k];
|
||
}
|
||
}
|
||
|
||
|
||
T rawinterp(int64_t jl, T x) override{
|
||
|
||
int64_t klo=jl, khi=jl+1;
|
||
T y, h, b, a;
|
||
|
||
h = xx[khi] - xx[klo];
|
||
if (h == T{0}){ // The xa’s must be distinct.
|
||
throw std::invalid_argument("interp_cubic_spline: Bad input to routine splint");
|
||
}
|
||
|
||
a = (xx[khi] - x)/h; // Cubic spline polynomial is now evaluated.
|
||
b = (x - xx[klo])/h;
|
||
y = a*yy[klo] + b*yy[khi] + ( ((a*a*a) - a)*y2[klo] + ((b*b*b) - b)*y2[khi] ) * (h*h) / T{6};
|
||
|
||
|
||
return y;
|
||
}
|
||
|
||
};
|
||
|
||
} // namespace numerics
|
||
|