1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119 | /*
* spline.h
*
* simple cubic spline interpolation library without external
* dependencies
*
* ---------------------------------------------------------------------
* Copyright (C) 2011, 2014 Tino Kluge (ttk448 at gmail.com)
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* 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 General Public License for more details.
*
* You should have received a copy of the GNU 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>
// unnamed namespace only because the implementation is in this
// header file and we don't want to export symbols to the obj files
//namespace
//{
namespace tk
{
// band matrix solver
class band_matrix
{
private:
std::vector< std::vector<double> > m_upper; // upper band
std::vector< std::vector<double> > m_lower; // lower band
public:
band_matrix() {}; // constructor
band_matrix(int dim, int n_u, int n_l); // constructor
~band_matrix() {}; // destructor
void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l
int dim() const; // matrix dimension
int num_upper() const
{
return static_cast<int>(m_upper.size())-1;
}
int num_lower() const
{
return static_cast<int>(m_lower.size())-1;
}
// 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)
double& saved_diag(int i);
double saved_diag(int i) const;
void lu_decompose();
std::vector<double> r_solve(const std::vector<double>& b) const;
std::vector<double> l_solve(const std::vector<double>& b) const;
std::vector<double> lu_solve(const std::vector<double>& b,
bool is_lu_decomposed=false);
};
// spline interpolation
class spline
{
public:
enum bd_type {
first_deriv = 1,
second_deriv = 2
};
private:
std::vector<double> m_x,m_y; // x,y coordinates of points
// 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, m_right;
double m_left_value, m_right_value;
bool m_force_linear_extrapolation;
public:
// set default boundary condition to be zero curvature at both ends
spline(): m_left(second_deriv), m_right(second_deriv),<--- Member variable 'spline::m_b0' is not initialized in the constructor.<--- Member variable 'spline::m_c0' is not initialized in the constructor.
m_left_value(0.0), m_right_value(0.0),
m_force_linear_extrapolation(false)
{
;
}
// 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 std::vector<double>& x,
const std::vector<double>& y, bool cubic_spline=true);
double operator() (double x) const;
double deriv(int order, double x) const;
};
} // namespace tk
//}
#endif /* TK_SPLINE_H */
|