libctsupport/cubicinterp.cpp

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (c) 1983-2001 Kevin Rosenberg
   4 **
   5 **  $Id: cubicinterp.cpp,v 1.4 2001/02/09 21:27:51 kevin Exp $
   6 **
   7 **  This program is free software; you can redistribute it and/or modify
   8 **  it under the terms of the GNU General Public License (version 2) as
   9 **  published by the Free Software Foundation.
  10 **
  11 **  This program is distributed in the hope that it will be useful,
  12 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 **  GNU General Public License for more details.
  15 **
  16 **  You should have received a copy of the GNU General Public License
  17 **  along with this program; if not, write to the Free Software
  18 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19 ******************************************************************************/
  20
  21
  22 #include "ctsupport.h"
  23 #include "cubicinterp.h"
  24
  25
  26 CubicInterpolator::CubicInterpolator (const double* const y, const int n)
  27   : m_pdY(y), m_n(n)
  28 {
  29   // Precalculate 2nd derivative of y and put in m_pdY2
  30   // Calculated by solving set of simultaneous cubic spline equations
  31   // Only n-2 cubic spline equations, but able to make two more
  32   // equations by setting second derivative to 0 at ends
  33
  34   m_pdY2 = new double [n];
  35   m_pdY2[0] = 0;   // second deriviative = 0 at beginning and end
  36   m_pdY2[n-1] = 0;
  37
  38   double* temp = new double [n - 1];
  39   temp[0] = 0;
  40   int i;
  41   for (i = 1; i < n - 1; i++) {
  42     double t = 2 + (0.5 * m_pdY2[i-1]);
  43     temp[i] = y[i+1] + y[i-1] - y[i] - y[i];
  44     temp[i] = (3 * temp[i] - 0.5 * temp[i-1]) / t;
  45     m_pdY2[i] = -0.5 / t;
  46   }
  47
  48   for (i = n - 2; i >= 0; i--)
  49     m_pdY2[i] = temp[i] + m_pdY2[i] * m_pdY2[i + 1];
  50
  51   delete temp;
  52 }
  53
  54 CubicInterpolator::~CubicInterpolator ()
  55 {
  56   delete m_pdY2;
  57 }
  58
  59
  60 double
  61 CubicInterpolator::interpolate (double x)
  62 {
  63   const static double oneSixth = (1. / 6.);
  64   int lo = static_cast<int>(floor(x));
  65   int hi = lo + 1;
  66
  67   if (lo < 0 || hi >= m_n) {
  68     sys_error (ERR_SEVERE, "X range out of bounds [CubicInterpolator::interpolate]");
  69     return (0);
  70   }
  71
  72   double loFr = hi - x;
  73   double hiFr = 1 - loFr;
  74   double y = loFr * m_pdY[lo] + hiFr * m_pdY[hi];
  75   y += oneSixth * ((loFr*loFr*loFr - loFr) * m_pdY2[lo] + (hiFr*hiFr*hiFr - hiFr) * m_pdY2[hi]);
  76
  77   return y;
  78 }
  79
  80
  81