projects
/
ctsim.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
r514: no message
[ctsim.git]
/
libctsupport
/
cubicinterp.cpp
diff --git
a/libctsupport/cubicinterp.cpp
b/libctsupport/cubicinterp.cpp
index 8392dec898d93412c4bcc20758ae4f37fb1f9846..8d0c62cbc9dd28bbac59a5d43db85ba9c8d94d9b 100644
(file)
--- a/
libctsupport/cubicinterp.cpp
+++ b/
libctsupport/cubicinterp.cpp
@@
-2,7
+2,7
@@
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: cubicinterp.cpp,v 1.
1 2001/02/09 01:54:2
1 kevin Exp $
+** $Id: cubicinterp.cpp,v 1.
4 2001/02/09 21:27:5
1 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
**
** 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
@@
-26,13
+26,19
@@
CubicInterpolator::CubicInterpolator (const double* const y, const int n)
: m_pdY(y), m_n(n)
{
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 = new double [n];
- m_pdY2[0] = 0
.0
; // second deriviative = 0 at beginning and end
+ m_pdY2[0] = 0; // second deriviative = 0 at beginning and end
m_pdY2[n-1] = 0;
double* temp = new double [n - 1];
m_pdY2[n-1] = 0;
double* temp = new double [n - 1];
- temp[0] = 0.0
- for (int i = 1; i < n - 1; i++) {
+ 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;
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;
@@
-64,7
+70,7
@@
CubicInterpolator::interpolate (double x)
}
double loFr = hi - x;
}
double loFr = hi - x;
- double hiFr = 1 -
a
;
+ 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]);
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]);