mirror of
https://github.com/gwm17/catima.git
synced 2024-11-23 02:38:51 -05:00
94 lines
2.9 KiB
C++
94 lines
2.9 KiB
C++
/*
|
|
*
|
|
* This is modification of Tino Kluge tk spline
|
|
* calculation is optimized for tridiagonal matrices
|
|
*
|
|
* Copyright(C) 2017
|
|
* This program is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU Affero General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Affero General Public License for more details.
|
|
|
|
* You should have received a copy of the GNU Affero General Public License
|
|
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#ifndef TK_SPLINE_H
|
|
#define TK_SPLINE_H
|
|
|
|
#include <cstdio>
|
|
#include <cassert>
|
|
#include <vector>
|
|
#include <algorithm>
|
|
|
|
namespace catima
|
|
{
|
|
|
|
// band matrix solver
|
|
class band_matrix
|
|
{
|
|
private:
|
|
std::vector<double> a;
|
|
std::vector<double> d;
|
|
std::vector<double> c;
|
|
std::vector<double> save;
|
|
public:
|
|
band_matrix() {}; // constructor
|
|
band_matrix(int dim); // constructor
|
|
~band_matrix() {}; // destructor
|
|
void resize(int dim); // init with dim,n_u,n_l
|
|
int dim() const; // matrix dimension
|
|
|
|
// access operator
|
|
double & operator () (int i, int j); // write
|
|
double operator () (int i, int j) const; // read
|
|
// we can store an additional diogonal (in m_lower)
|
|
std::vector<double> trig_solve(const std::vector<double>& b) const;
|
|
};
|
|
|
|
|
|
// spline interpolation
|
|
class spline
|
|
{
|
|
public:
|
|
enum class bd_type {
|
|
first_deriv = 1,
|
|
second_deriv = 2
|
|
};
|
|
|
|
private:
|
|
const double *m_x, *m_y; // x,y coordinates of points
|
|
size_t n=0;
|
|
// interpolation parameters
|
|
// f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
|
|
std::vector<double> m_a,m_b,m_c; // spline coefficients
|
|
double m_b0, m_c0; // for left extrapol
|
|
bd_type m_left = bd_type::second_deriv;
|
|
bd_type m_right = bd_type::second_deriv;
|
|
double m_left_value = 0.0;
|
|
double m_right_value = 0.0;
|
|
bool m_force_linear_extrapolation = false;
|
|
|
|
public:
|
|
// set default boundary condition to be zero curvature at both ends
|
|
spline(){}
|
|
|
|
// optional, but if called it has to come be before set_points()
|
|
void set_boundary(bd_type left, double left_value,
|
|
bd_type right, double right_value,
|
|
bool force_linear_extrapolation=false);
|
|
void set_points(const double *x,
|
|
const double *y,
|
|
const size_t num);
|
|
double operator() (double x) const;
|
|
double deriv(int order, double x) const;
|
|
};
|
|
|
|
} // namespace tk
|
|
|
|
#endif /* TK_SPLINE_H */
|