source: trunk/src/CubicSplineInterpolator1D.tcc@ 2891

Last change on this file since 2891 was 2850, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): sd

Description: Describe your changes here...

Quite minor (and non-effective) bug fix.


File size: 4.1 KB
Line 
1//
2// C++ Implementation: CubicSplineInterpolator1D
3//
4// Description:
5//
6//
7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2012
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <assert.h>
13
14#include <casa/Exceptions/Error.h>
15#include <casa/Utilities/Assert.h>
16
17#include <iostream>
18using namespace std;
19
20#include "CubicSplineInterpolator1D.h"
21
22namespace asap {
23
24template <class T, class U>
25CubicSplineInterpolator1D<T, U>::CubicSplineInterpolator1D()
26 : Interpolator1D<T, U>(),
27 y2_(0),
28 ny2_(0),
29 reusable_(false)
30{}
31
32template <class T, class U>
33CubicSplineInterpolator1D<T, U>::~CubicSplineInterpolator1D()
34{
35 if (y2_)
36 delete[] y2_;
37}
38
39template <class T, class U>
40void CubicSplineInterpolator1D<T, U>::setData(T *x, U *y, unsigned int n)
41{
42 Interpolator1D<T, U>::setData(x, y, n);
43 reusable_ = false;
44}
45
46template <class T, class U>
47void CubicSplineInterpolator1D<T, U>::setX(T *x, unsigned int n)
48{
49 Interpolator1D<T, U>::setX(x, n);
50 reusable_ = false;
51}
52
53template <class T, class U>
54void CubicSplineInterpolator1D<T, U>::setY(U *y, unsigned int n)
55{
56 Interpolator1D<T, U>::setY(y, n);
57 reusable_ = false;
58}
59
60template <class T, class U>
61U CubicSplineInterpolator1D<T, U>::interpolate(T x)
62{
63 //assert(this->isready());
64 assert_<casa::AipsError>(this->isready(), "object is not ready to process.");
65 if (this->n_ == 1)
66 return this->y_[0];
67
68 unsigned int i = this->locator_->locate(x);
69
70 // do not perform extrapolation
71 if (i == 0) {
72 return this->y_[i];
73 }
74 else if (i == this->n_) {
75 return this->y_[i-1];
76 }
77
78 // determine second derivative of each point
79 if (!reusable_) {
80 evaly2();
81 reusable_ = true;
82 }
83
84 // cubic spline interpolation
85 float y = dospline(x, i);
86 return y;
87}
88
89template <class T, class U>
90void CubicSplineInterpolator1D<T, U>::evaly2()
91{
92 if (this->n_ > ny2_) {
93 if (y2_)
94 delete[] y2_;
95 y2_ = new U[this->n_];
96 ny2_ = this->n_;
97 }
98
99 U *u = new U[ny2_-1];
100 unsigned int *idx = new unsigned int[this->n_];
101
102 // Natural cubic spline.
103 y2_[0] = 0.0;
104 y2_[ny2_-1] = 0.0;
105 u[0] = 0.0;
106
107 if (this->x_[0] < this->x_[this->n_-1]) {
108 // ascending
109 for (unsigned int i = 0; i < this->n_; ++i)
110 idx[i] = i;
111 }
112 else {
113 // descending
114 for (unsigned int i = 0; i < this->n_; ++i)
115 idx[i] = this->n_ - 1 - i;
116 }
117
118
119 // Solve tridiagonal system.
120 // Here, tridiagonal matrix is decomposed to upper triangular
121 // matrix. u stores upper triangular components while y2_ stores
122 // right-hand side vector. The diagonal elements are normalized
123 // to 1.
124 T a1 = this->x_[idx[1]] - this->x_[idx[0]];
125 T a2, bi;
126 for (unsigned int i = 1; i < ny2_ - 1; i++) {
127 a2 = this->x_[idx[i+1]] - this->x_[idx[i]];
128 bi = 1.0 / (this->x_[idx[i+1]] - this->x_[idx[i-1]]);
129 y2_[i] = 3.0 * bi * ((this->y_[idx[i+1]] - this->y_[idx[i]]) / a2
130 - (this->y_[idx[i]] - this->y_[idx[i-1]]) / a1
131 - y2_[i-1] * 0.5 * a1);
132 a1 = 1.0 / (1.0 - u[i-1] * 0.5 * a1 * bi);
133 y2_[i] *= a1;
134 u[i] = 0.5 * a2 * bi * a1;
135 a1 = a2;
136 }
137
138 // Then, solve the system by backsubstitution and store solution
139 // vector to y2_.
140 for (int k = ny2_ - 2; k >= 1; k--)
141 y2_[k] -= u[k] * y2_[k+1];
142
143 delete[] idx;
144 delete[] u;
145}
146
147template <class T, class U>
148U CubicSplineInterpolator1D<T, U>::dospline(T x, unsigned int i)
149{
150 unsigned int index_lower;
151 unsigned int index_higher;
152 unsigned int index_lower_correct;
153 unsigned int index_higher_correct;
154 if (this->x_[0] < this->x_[this->n_-1]) {
155 index_lower = i - 1;
156 index_higher = i;
157 index_lower_correct = index_lower;
158 index_higher_correct = index_higher;
159 }
160 else {
161 index_lower = i;
162 index_higher = i - 1;
163 index_lower_correct = this->n_ - 1 - index_lower;
164 index_higher_correct = this->n_ - 1 - index_higher;
165 }
166 T dx = this->x_[index_higher] - this->x_[index_lower];
167 T a = (this->x_[index_higher] - x) / dx;
168 T b = (x - this->x_[index_lower]) / dx;
169 U y = a * this->y_[index_lower] + b * this->y_[index_higher] +
170 ((a * a * a - a) * y2_[index_lower_correct] + (b * b * b - b) * y2_[index_higher_correct]) * (dx * dx) / 6.0;
171 return y;
172}
173
174}
Note: See TracBrowser for help on using the repository browser.