Changeset 2849


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/CubicSplineInterpolator1D.tcc

    r2756 r2849  
    9898
    9999  U *u = new U[ny2_-1];
     100  unsigned int *idx = new unsigned int[this->n_];
    100101
    101102  // Natural cubic spline.
     
    104105  u[0] = 0.0;
    105106
     107  if (this->x_[0] < this->x_[this->n_-1]) {
     108    // ascending
     109    for (unsigned int i = 0; i < this->n_; ++i)
     110      idx[i] = i;
     111  }
     112  else {
     113    // descending
     114    for (unsigned int i = 0; i < this->n_; ++i)
     115      idx[i] = this->n_ - 1 - i;
     116  }
     117
     118
    106119  // Solve tridiagonal system.
    107120  // Here, tridiagonal matrix is decomposed to upper triangular
     
    109122  // right-hand side vector. The diagonal elements are normalized
    110123  // to 1.
    111   T a1 = this->x_[1] - this->x_[0];
     124  T a1 = this->x_[idx[1]] - this->x_[idx[0]];
    112125  T a2, bi;
    113126  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
     127    a2 = this->x_[idx[i+1]] - this->x_[idx[i]];
     128    bi = 1.0 / (this->x_[idx[i+1]] - this->x_[idx[i-1]]);
     129    y2_[i] = 3.0 * bi * ((this->y_[idx[i+1]] - this->y_[idx[i]]) / a2
     130                         - (this->y_[idx[i]] - this->y_[idx[i-1]]) / a1
    118131                         - y2_[i-1] * 0.5 * a1);
    119132    a1 = 1.0 / (1.0 - u[i-1] * 0.5 * a1 * bi);
     
    127140  for (int k = ny2_ - 2; k >= 0; k--)
    128141    y2_[k] -= u[k] * y2_[k+1];
    129 
     142 
     143  delete[] idx;
    130144  delete[] u;
    131145}
     
    134148U CubicSplineInterpolator1D<T, U>::dospline(T x, unsigned int i)
    135149{
    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;
     150  unsigned int index_lower;
     151  unsigned int index_higher;
     152  unsigned int index_lower_correct;
     153  unsigned int index_higher_correct;
     154  if (this->x_[0] < this->x_[this->n_-1]) {
     155    index_lower = i - 1;
     156    index_higher = i;
     157    index_lower_correct = index_lower;
     158    index_higher_correct = index_higher;
     159  }
     160  else {
     161    index_lower = i;
     162    index_higher = i - 1;
     163    index_lower_correct = this->n_ - 1 - index_lower;
     164    index_higher_correct = this->n_ - 1 - index_higher;
     165  }
     166  T dx = this->x_[index_higher] - this->x_[index_lower];
     167  T a = (this->x_[index_higher] - x) / dx;
     168  T b = (x - this->x_[index_lower]) / dx;
     169  U y = a * this->y_[index_lower] + b * this->y_[index_higher] +
     170    ((a * a * a - a) * y2_[index_lower_correct] + (b * b * b - b) * y2_[index_higher_correct]) * (dx * dx) / 6.0;
    142171  return y;
    143172}
Note: See TracChangeset for help on using the changeset viewer.