1 /*****************************************************************************
2 ** This is part of the CTSim program
3 ** Copyright (c) 1983-2001 Kevin Rosenberg
5 ** $Id: interpolator.cpp,v 1.1 2001/02/11 21:57:08 kevin Exp $
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.
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.
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 ******************************************************************************/
22 #include "ctsupport.h"
23 #include "interpolator.h"
26 CubicPolyInterpolator::CubicPolyInterpolator (const double* const y, const int n)
30 sys_error (ERR_SEVERE, "Too few points (%d) in CubicPolyInterpolator", m_n);
33 CubicPolyInterpolator::~CubicPolyInterpolator ()
39 CubicPolyInterpolator::interpolate (double x)
41 int lo = static_cast<int>(floor(x)) - 1;
45 sys_error (ERR_WARNING, "x=%f, out of range [CubicPolyInterpolator]", x);
47 } else if (lo == -1) // linear interpolate at between x = 0 & 1
48 return m_pdY[0] + x * (m_pdY[1] - m_pdY[0]);
51 sys_error (ERR_WARNING, "x=%f, out of range [CubicPolyInterpolator]", x);
53 } else if (hi == m_n) {// linear interpolate between x = (n-2) and (n-1)
54 double frac = x - (lo + 1);
55 return m_pdY[m_n - 2] + frac * (m_pdY[m_n - 1] - m_pdY[m_n - 2]);
58 // Lagrange formula for N=4 (cubic)
61 double xd_1 = x - (lo + 1);
62 double xd_2 = x - (lo + 2);
63 double xd_3 = x - (lo + 3);
65 static double oneSixth = (1. / 6.);
67 double y = xd_1 * xd_2 * xd_3 * -oneSixth * m_pdY[lo];
68 y += xd_0 * xd_2 * xd_3 * 0.5 * m_pdY[lo+1];
69 y += xd_0 * xd_1 * xd_3 * -0.5 * m_pdY[lo+2];
70 y += xd_0 * xd_1 * xd_2 * oneSixth * m_pdY[lo+3];
77 CubicSplineInterpolator::CubicSplineInterpolator (const double* const y, const int n)
80 // Precalculate 2nd derivative of y and put in m_pdY2
81 // Calculated by solving set of simultaneous CubicSpline spline equations
82 // Only n-2 CubicSpline spline equations, but able to make two more
83 // equations by setting second derivative to 0 at ends
85 m_pdY2 = new double [n];
86 m_pdY2[0] = 0; // second deriviative = 0 at beginning and end
89 double* temp = new double [n - 1];
92 for (i = 1; i < n - 1; i++) {
93 double t = 2 + (0.5 * m_pdY2[i-1]);
94 temp[i] = y[i+1] + y[i-1] - y[i] - y[i];
95 temp[i] = (3 * temp[i] - 0.5 * temp[i-1]) / t;
99 for (i = n - 2; i >= 0; i--)
100 m_pdY2[i] = temp[i] + m_pdY2[i] * m_pdY2[i + 1];
105 CubicSplineInterpolator::~CubicSplineInterpolator ()
112 CubicSplineInterpolator::interpolate (double x)
114 const static double oneSixth = (1. / 6.);
115 int lo = static_cast<int>(floor(x));
118 if (lo < 0 || hi >= m_n) {
119 sys_error (ERR_SEVERE, "X range out of bounds [CubicSplineInterpolator::interpolate]");
123 double loFr = hi - x;
124 double hiFr = 1 - loFr;
125 double y = loFr * m_pdY[lo] + hiFr * m_pdY[hi];
126 y += oneSixth * ((loFr*loFr*loFr - loFr) * m_pdY2[lo] + (hiFr*hiFr*hiFr - hiFr) * m_pdY2[hi]);