source: trunk/src/PolynomialInterpolator1D.tcc@ 2859

Last change on this file since 2859 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: 2.7 KB
Line 
1//
2// C++ Implementation: PolynomialInterpolator1D
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#include <math.h>
14#include <iostream>
15using namespace std;
16
17#include <casa/Exceptions/Error.h>
18#include <casa/Utilities/Assert.h>
19
20#include "PolynomialInterpolator1D.h"
21
22namespace asap {
23
24template <class T, class U>
25PolynomialInterpolator1D<T, U>::PolynomialInterpolator1D()
26 : Interpolator1D<T, U>()
27{}
28
29template <class T, class U>
30PolynomialInterpolator1D<T, U>::~PolynomialInterpolator1D()
31{}
32
33template <class T, class U>
34U PolynomialInterpolator1D<T, U>::interpolate(T x)
35{
36 //casa::AlwaysAssert((this->isready()),(casa::AipsError));
37 casa::assert_<casa::AipsError>(this->isready(), "object is not ready to process.");
38 if (this->n_ == 1)
39 return this->y_[0];
40
41 unsigned int i = this->locator_->locate(x);
42
43 // do not perform extrapolation
44 if (i == 0) {
45 return this->y_[i];
46 }
47 else if (i == this->n_) {
48 return this->y_[i-1];
49 }
50
51 // polynomial interpolation
52 U y;
53 if (this->order_ >= this->n_ - 1) {
54 // use full region
55 y = dopoly(x, 0, this->n_);
56 }
57 else {
58 // use sub-region
59 int j = i - 1 - this->order_ / 2;
60 unsigned int m = this->n_ - 1 - this->order_;
61 unsigned int k = (unsigned int)((j > 0) ? j : 0);
62 k = ((k > m) ? m : k);
63 y = dopoly(x, k, this->order_ + 1);
64 }
65
66 return y;
67}
68
69template <class T, class U>
70U PolynomialInterpolator1D<T, U>::dopoly(T x, unsigned int left,
71 unsigned int n)
72{
73 T *xa = &this->x_[left];
74 U *ya = &this->y_[left];
75
76 // storage for C and D in Neville's algorithm
77 U *c = new U[n];
78 U *d = new U[n];
79 for (unsigned int i = 0; i < n; i++) {
80 c[i] = ya[i];
81 d[i] = ya[i];
82 }
83
84 // Neville's algorithm
85 U y = c[0];
86 for (unsigned int m = 1; m < n; m++) {
87 // Evaluate Cm1, Cm2, Cm3, ... Cm[n-m] and Dm1, Dm2, Dm3, ... Dm[n-m].
88 // Those are stored to c[0], c[1], ..., c[n-m-1] and d[0], d[1], ...,
89 // d[n-m-1].
90 for (unsigned int i = 0; i < n - m; i++) {
91 U cd = c[i+1] - d[i];
92 T dx = xa[i] - xa[i+m];
93 try {
94 cd /= (U)dx;
95 }
96 catch (...) {
97 delete[] c;
98 delete[] d;
99 throw casa::AipsError("x_ has duplicate elements");
100 }
101 c[i] = (xa[i] - x) * cd;
102 d[i] = (xa[i+m] - x) * cd;
103 }
104
105 // In each step, c[0] holds Cm1 which is a correction between
106 // P12...m and P12...[m+1]. Thus, the following repeated update
107 // corresponds to the route P1 -> P12 -> P123 ->...-> P123...n.
108 y += c[0];
109 }
110
111 delete[] c;
112 delete[] d;
113
114 return y;
115}
116}
Note: See TracBrowser for help on using the repository browser.