cubic_spline.cxx

Go to the documentation of this file.

#include <cxxtls/cubic_spline.h>
#include <vector>

namespace cxxtls
{


// this routine was originally fortran, I wrote it in the following
// goofy way to keep the algorithm the same.  My arrays don't use
// location 0 and the indexing thus runs from 1 to n like fortran.

// This implementation of cubic spline came from the book
//
//  Numerical Computing:  an introduction
//  by Shampine and Allen, 1973.

// yes this file needs to be rewrite not to have to copy these arrays...

void 
cubic_spline::
compute_coefficients()
{
  using std::vector;

  dirty_ = false;

  size_t n = rep_.size();

  // copy the map to arrays, doh!

  vector<double> xn(n+1), fn(n+1), rho(n+1), tau(n+1), s(n+1); 

    {
      const_iterator first = rep_.begin(),
                                      last  = rep_.end();

      vector<double>::iterator xout = xn.begin() + 1;
      vector<double>::iterator yout = fn.begin() + 1;

      while(first != last)
      {
                 *xout++ = first->first;
                 *yout++ = first->second.y_;

                 ++first;
      }

    }

  // the two vectors, xn and fn are have undefined values in location 0

  size_t nm1 = n-1;
  
  rho[2] = 0;
  tau[2] = 0;

  int i;

  for(i=2; i <= static_cast<int>(nm1); ++i)
  {
    size_t iim1 = i-1;
    size_t ii   = i;
    size_t iip1 = i+1;

    double him1 = xn[ii] - xn[iim1];
    double hi   = xn[iip1]-xn[ii];
    double temp = (him1/hi)*(rho[i]+2)+2;
                                   
    rho[i+1]=-1.0/temp;

    double d = 6.0 *((fn[iip1]-fn[ii])/hi-(fn[ii]-fn[iim1])/him1)/hi;

    tau[i+1] = (d-him1*tau[i]/hi)/temp;

  }

  s[1] = 0;
  s[n] = 0;

  int nm2 = n-2;

  for(i=1; i <= nm2; ++i)
  {
    size_t ib = n-i;
    s[ib] = rho[ib+1] * s[ib+1] +tau[ib+1];
  }

  {
    vector<double>::iterator first = s.begin() + 1,
                                                      last  = s.end();

    rep_t::iterator output = rep_.begin();

    while(first != last)
    {
      output->second.s_ = *first++;
      ++output;
    }

  }


}

double 
cubic_spline::
evaluate(double x) const
{
  using std::vector;

  if(dirty_)
  {
    cubic_spline* me = const_cast<cubic_spline*>(this);
    me->compute_coefficients();
  }

  size_t n = rep_.size();
  
  if(n == 0)
    return 0;

  if(rep_.size() == 1)
  {
    return rep_.begin()->second.y_;
  }

  if( x < rep_.begin()->first )
  {
    const_iterator lb = rep_.begin();
    const_iterator ub = lb; ++ub;

    double m = (ub->second.y_ - lb->second.y_) / 
                        (ub->first     - lb->first );

    double deltax = lb->first  - x;
    double deltay = m * deltax;

    return lb->second.y_ - deltay;
  }

  rep_t::const_iterator end = rep_.end(); --end;

  if( x >= end->first)
  {
    const_iterator ub = rep_.end(); --ub;
    const_iterator lb = ub; --lb;

    double m = (ub->second.y_ - lb->second.y_) / 
                        (ub->first     - lb->first );

    double deltax = x - ub->first ;
    double deltay = m * deltax;

    return ub->second.y_ + deltay;
  }

  const_iterator ilp1 = ceil(x);
  const_iterator il   = floor(x);

  double a = ilp1->first - x;
  double b = x - il->first;
  double h = ilp1->first - il->first;

  return  a * il->second.s_ * (a * a / h - h)/6.0 +
                   b * ilp1->second.s_ * (b * b / h - h)/6.0 +
                   ( a * il->second.y_ + b * ilp1->second.y_)/h;

}

cubic_spline::const_iterator
cubic_spline::
floor(double x) const
  //

{
  const_iterator ub    = rep_.upper_bound(x);
  const_iterator begin = rep_.begin();

  while(ub != begin)
  {
    --ub;                    // at most one iteration is realistic

    if(ub->first <= x)
      return ub;

  }

  return rep_.end();

}

} // namespace cxxtls