source: trunk/src/CubicSplineInterpolator1D.tcc@ 2804

Last change on this file since 2804 was 2756, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-4770

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdcal2

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Replace assert with casa::assert_ to avoid aborting casapy session.


File size: 3.2 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
101 // Natural cubic spline.
102 y2_[0] = 0.0;
103 y2_[ny2_-1] = 0.0;
104 u[0] = 0.0;
105
106 // Solve tridiagonal system.
107 // Here, tridiagonal matrix is decomposed to upper triangular
108 // matrix. u stores upper triangular components while y2_ stores
109 // right-hand side vector. The diagonal elements are normalized
110 // to 1.
111 T a1 = this->x_[1] - this->x_[0];
112 T a2, bi;
113 for (unsigned int i = 1; i < ny2_ - 1; i++) {
114 a2 = this->x_[i+1] - this->x_[i];
115 bi = 1.0 / (this->x_[i+1] - this->x_[i-1]);
116 y2_[i] = 3.0 * bi * ((this->y_[i+1] - this->y_[i]) / a2
117 - (this->y_[i] - this->y_[i-1]) / a1
118 - y2_[i-1] * 0.5 * a1);
119 a1 = 1.0 / (1.0 - u[i-1] * 0.5 * a1 * bi);
120 y2_[i] *= a1;
121 u[i] = 0.5 * a2 * bi * a1;
122 a1 = a2;
123 }
124
125 // Then, solve the system by backsubstitution and store solution
126 // vector to y2_.
127 for (int k = ny2_ - 2; k >= 0; k--)
128 y2_[k] -= u[k] * y2_[k+1];
129
130 delete[] u;
131}
132
133template <class T, class U>
134U CubicSplineInterpolator1D<T, U>::dospline(T x, unsigned int i)
135{
136 unsigned int j = i - 1;
137 T h = this->x_[i] - this->x_[j];
138 T a = (this->x_[i] - x) / h;
139 T b = (x - this->x_[j]) / h;
140 U y = a * this->y_[j] + b * this->y_[i] +
141 ((a * a * a - a) * y2_[j] + (b * b * b - b) * y2_[i]) * (h * h) / 6.0;
142 return y;
143}
144
145}
Note: See TracBrowser for help on using the repository browser.