Changeset 2849

Timestamp:

09/05/13 16:25:53 (11 years ago)

Author:

Takeshi Nakazato

Message:

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: No

Module(s): sd

Description: Describe your changes here...

Fixed a bug in spline interpolation that causes slight error for
LSB spectral window.

File:

• trunk/src/CubicSplineInterpolator1D.tcc (modified) (5 diffs)

trunk/src/CubicSplineInterpolator1D.tcc

@@ -98 +98 @@
 
   U *u = new U[ny2_-1];
+  unsigned int *idx = new unsigned int[this->n_];
 
   // Natural cubic spline.
@@ -104 +105 @@
   u[0] = 0.0;
 
+  if (this->x_[0] < this->x_[this->n_-1]) {
+    // ascending
+    for (unsigned int i = 0; i < this->n_; ++i)
+      idx[i] = i;
+  }
+  else {
+    // descending
+    for (unsigned int i = 0; i < this->n_; ++i)
+      idx[i] = this->n_ - 1 - i;
+  }
+
+
   // Solve tridiagonal system.
   // Here, tridiagonal matrix is decomposed to upper triangular 
@@ -109 +122 @@
   // right-hand side vector. The diagonal elements are normalized 
   // to 1.
-  T a1 = this->x_[1] - this->x_[0];
+  T a1 = this->x_[idx[1]] - this->x_[idx[0]];
   T a2, bi;
   for (unsigned int i = 1; i < ny2_ - 1; i++) {
-    a2 = this->x_[i+1] - this->x_[i];
-    bi = 1.0 / (this->x_[i+1] - this->x_[i-1]);
-    y2_[i] = 3.0 * bi * ((this->y_[i+1] - this->y_[i]) / a2 
-                         - (this->y_[i] - this->y_[i-1]) / a1 
+    a2 = this->x_[idx[i+1]] - this->x_[idx[i]];
+    bi = 1.0 / (this->x_[idx[i+1]] - this->x_[idx[i-1]]);
+    y2_[i] = 3.0 * bi * ((this->y_[idx[i+1]] - this->y_[idx[i]]) / a2 
+                         - (this->y_[idx[i]] - this->y_[idx[i-1]]) / a1 
                          - y2_[i-1] * 0.5 * a1);
     a1 = 1.0 / (1.0 - u[i-1] * 0.5 * a1 * bi);
@@ -127 +140 @@
   for (int k = ny2_ - 2; k >= 0; k--)
     y2_[k] -= u[k] * y2_[k+1];
-
+  
+  delete[] idx;
   delete[] u;
 }
@@ -134 +148 @@
 U CubicSplineInterpolator1D<T, U>::dospline(T x, unsigned int i)
 {
-  unsigned int j = i - 1;
-  T h = this->x_[i] - this->x_[j];
-  T a = (this->x_[i] - x) / h;
-  T b = (x - this->x_[j]) / h;
-  U y = a * this->y_[j] + b * this->y_[i] + 
-    ((a * a * a - a) * y2_[j] + (b * b * b - b) * y2_[i]) * (h * h) / 6.0;
+  unsigned int index_lower;
+  unsigned int index_higher;
+  unsigned int index_lower_correct;
+  unsigned int index_higher_correct;
+  if (this->x_[0] < this->x_[this->n_-1]) {
+    index_lower = i - 1;
+    index_higher = i;
+    index_lower_correct = index_lower;
+    index_higher_correct = index_higher;
+  }
+  else {
+    index_lower = i;
+    index_higher = i - 1;
+    index_lower_correct = this->n_ - 1 - index_lower;
+    index_higher_correct = this->n_ - 1 - index_higher;
+  }
+  T dx = this->x_[index_higher] - this->x_[index_lower];
+  T a = (this->x_[index_higher] - x) / dx;
+  T b = (x - this->x_[index_lower]) / dx;
+  U y = a * this->y_[index_lower] + b * this->y_[index_higher] + 
+    ((a * a * a - a) * y2_[index_lower_correct] + (b * b * b - b) * y2_[index_higher_correct]) * (dx * dx) / 6.0;
   return y;
 }