source: trunk/src/PolynomialInterpolator1D.tcc@ 2736

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

New Development: No

JIRA Issue: Yes CAS-4770

Ready for Test: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Redefined Interpolator1D and derived classes as template class.


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