Changeset 2730 for trunk/src/PolynomialInterpolator1D.cpp

Timestamp:

01/16/13 16:00:28 (11 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: Yes CAS-4770

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): Module Names change impacts.

Description: Describe your changes here...

Rewrite implementations for locator and interpolator.
Documentation (doxygen format) is added to header files.

File:

• trunk/src/PolynomialInterpolator1D.cpp (modified) (4 diffs)

trunk/src/PolynomialInterpolator1D.cpp

@@ -12 +12 @@
 #include <assert.h>
 #include <math.h>
+#include <iostream>
+using namespace std;
+
+#include <casa/Exceptions/Error.h>
 
 #include "PolynomialInterpolator1D.h"
@@ -43 +47 @@
   float y;
   if (order_ >= n_ - 1) {
-    y = polint(x, i, 0, n_);
+    // use full region
+    y = dopoly(x, 0, n_);
   }
   else {
+    // use partial region
     int j = i - 1 - order_ / 2;
     unsigned int m = n_ - 1 - order_;
     unsigned int k = (unsigned int)((j > 0) ? j : 0);
     k = ((k > m) ? m : k);
-    y = polint(x, i, k, order_ + 1);
+    y = dopoly(x, k, order_ + 1);
   }
 
@@ -56 +62 @@
 }
 
-float PolynomialInterpolator1D::polint(double x, unsigned int loc, 
-                                       unsigned int left, unsigned int n)
+float PolynomialInterpolator1D::dopoly(double x, unsigned int left,
+                                       unsigned int n)
 {
-  int ns = loc - left;
-  if (fabs(x - x_[loc]) >= fabs(x - x_[loc-1])) {
-    ns--;
-  }
-
   double *xa = &x_[left];
   float *ya = &y_[left];
 
-  float y = ya[ns];
-  
-  // c stores C11, C21, C31, ..., C[n-1]1
-  // d stores D11, D21, D31, ..., D[n-1]1
+  // storage for C and D in Neville's algorithm
   float *c = new float[n];
   float *d = new float[n];
@@ -78 +76 @@
   }
 
+  // Neville's algorithm
+  float y = c[0];
   for (unsigned int m = 1; m < n; m++) {
+    // Evaluate Cm1, Cm2, Cm3, ... Cm[n-m] and Dm1, Dm2, Dm3, ... Dm[n-m].
+    // Those are stored to c[0], c[1], ..., c[n-m-1] and d[0], d[1], ..., 
+    // d[n-m-1].
     for (unsigned int i = 0; i < n - m; i++) {
-      double ho = xa[i] - x;
-      double hp = xa[i+m] - x;
-      float w = c[i+1] - d[i];
-      double denom = ho - hp;
-      if (denom == 0.0) {
+      float cd = c[i+1] - d[i];
+      double dx = xa[i] - xa[i+m];
+      try {
+        cd /= dx;
+      }
+      catch (...) {
         delete[] c;
         delete[] d;
-        assert(denom != 0.0);
+        throw casa::AipsError("x_ has duplicate elements");
       }
-      denom = w / denom;
-      
-      d[i] = hp * denom;
-      c[i] = ho * denom;
+      c[i] = (xa[i] - x) * cd;
+      d[i] = (xa[i+m] - x) * cd;
     }
 
-    float dy;
-    if (2 * ns < (int)(n - m)) {
-      dy = c[ns];
-    }
-    else {
-      dy = d[ns-1];
-      ns--;
-    }
-    y += dy;
+    // In each step, c[0] holds Cm1 which is a correction between 
+    // P12...m and P12...[m+1]. Thus, the following repeated update 
+    // corresponds to the route P1 -> P12 -> P123 ->...-> P123...n.
+    y += c[0];
   }
-  
+
   delete[] c;
   delete[] d;