r521: no message
[ctsim.git] / libctsupport / cubicinterp.cpp
diff --git a/libctsupport/cubicinterp.cpp b/libctsupport/cubicinterp.cpp
deleted file mode 100644 (file)
index 8d0c62c..0000000
+++ /dev/null
@@ -1,81 +0,0 @@
-/*****************************************************************************
-**  This is part of the CTSim program
-**  Copyright (c) 1983-2001 Kevin Rosenberg
-**
-**  $Id: cubicinterp.cpp,v 1.4 2001/02/09 21:27:51 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)
-{
-  // Precalculate 2nd derivative of y and put in m_pdY2
-  // Calculated by solving set of simultaneous cubic spline equations
-  // Only n-2 cubic spline equations, but able to make two more
-  // equations by setting second derivative to 0 at ends 
-
-  m_pdY2 = new double [n];
-  m_pdY2[0] = 0;   // second deriviative = 0 at beginning and end
-  m_pdY2[n-1] = 0;
-
-  double* temp = new double [n - 1];
-  temp[0] = 0;
-  int i;
-  for (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 - loFr;
-  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;
-}
-
-
-