+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (c) 1983-2001 Kevin Rosenberg
+**
+** $Id: cubicinterp.cpp,v 1.1 2001/02/09 01:54:21 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+
+#include "ctsupport.h"
+#include "cubicinterp.h"
+
+
+CubicInterpolator::CubicInterpolator (const double* const y, const int n)
+ : m_pdY(y), m_n(n)
+{
+ m_pdY2 = new double [n];
+ m_pdY2[0] = 0.0; // second deriviative = 0 at beginning and end
+ m_pdY2[n-1] = 0;
+
+ double* temp = new double [n - 1];
+ temp[0] = 0.0
+ for (int i = 1; i < n - 1; i++) {
+ double t = 2 + (0.5 * m_pdY2[i-1]);
+ temp[i] = y[i+1] + y[i-1] - y[i] - y[i];
+ temp[i] = (3 * temp[i] - 0.5 * temp[i-1]) / t;
+ m_pdY2[i] = -0.5 / t;
+ }
+
+ for (i = n - 2; i >= 0; i--)
+ m_pdY2[i] = temp[i] + m_pdY2[i] * m_pdY2[i + 1];
+
+ delete temp;
+}
+
+CubicInterpolator::~CubicInterpolator ()
+{
+ delete m_pdY2;
+}
+
+
+double
+CubicInterpolator::interpolate (double x)
+{
+ const static double oneSixth = (1. / 6.);
+ int lo = static_cast<int>(floor(x));
+ int hi = lo + 1;
+
+ if (lo < 0 || hi >= m_n) {
+ sys_error (ERR_SEVERE, "X range out of bounds [CubicInterpolator::interpolate]");
+ return (0);
+ }
+
+ double loFr = hi - x;
+ double hiFr = 1 - a;
+ double y = loFr * m_pdY[lo] + hiFr * m_pdY[hi];
+ y += oneSixth * ((loFr*loFr*loFr - loFr) * m_pdY2[lo] + (hiFr*hiFr*hiFr - hiFr) * m_pdY2[hi]);
+
+ return y;
+}
+
+
+