/*****************************************************************************
** This is part of the CTSim program
-** Copyright (c) 1983-2001 Kevin Rosenberg
-**
-** $Id: interpolator.cpp,v 1.1 2001/02/11 21:57:08 kevin Exp $
+** Copyright (c) 1983-2009 Kevin Rosenberg
**
** 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
int hi = lo + 3;
if (lo < -1) {
+#ifdef DEBUG
sys_error (ERR_WARNING, "x=%f, out of range [CubicPolyInterpolator]", x);
+#endif
return (0);
} else if (lo == -1) // linear interpolate at between x = 0 & 1
return m_pdY[0] + x * (m_pdY[1] - m_pdY[0]);
if (hi > m_n) {
+#ifdef DEBUG
sys_error (ERR_WARNING, "x=%f, out of range [CubicPolyInterpolator]", x);
+#endif
return (0);
} else if (hi == m_n) {// linear interpolate between x = (n-2) and (n-1)
double frac = x - (lo + 1);
}
// Lagrange formula for N=4 (cubic)
-
+
double xd_0 = x - lo;
double xd_1 = x - (lo + 1);
double xd_2 = x - (lo + 2);
// Precalculate 2nd derivative of y and put in m_pdY2
// Calculated by solving set of simultaneous CubicSpline spline equations
// Only n-2 CubicSpline spline equations, but able to make two more
- // equations by setting second derivative to 0 at ends
+ // 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
double* temp = new double [n - 1];
temp[0] = 0;
int i;
- for (i = 1; i < n - 1; 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--)
+ for (i = n - 2; i >= 0; i--)
m_pdY2[i] = temp[i] + m_pdY2[i] * m_pdY2[i + 1];
delete temp;
int hi = lo + 1;
if (lo < 0 || hi >= m_n) {
- sys_error (ERR_SEVERE, "X range out of bounds [CubicSplineInterpolator::interpolate]");
+#ifdef DEBUG
+ sys_error (ERR_SEVERE, "x out of bounds [CubicSplineInterpolator::interpolate]");
+#endif
return (0);
}
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;
}
+