001 package edu.nrao.sss.math;
002
003 import java.awt.geom.Point2D;
004 import java.util.ArrayList;
005 import java.util.List;
006
007 /**
008 * An interpolated value calculator that uses cubic spline interpolation.
009 * <p>
010 * <b>Version Info:</b>
011 * <table style="margin-left:2em">
012 * <tr><td>$Revision: 568 $</td></tr>
013 * <tr><td>$Date: 2007-04-27 16:53:09 -0600 (Fri, 27 Apr 2007) $</td></tr>
014 * <tr><td>$Author: dharland $</td></tr>
015 * </table></p>
016 *
017 * @author David M. Harland
018 * @since 2007-04-24
019 */
020 public class CubicSpline
021 extends AbsInterp
022 {
023 private List<Cubic> splines;
024
025 /** Creates a new instance. */
026 public CubicSpline()
027 {
028 super();
029
030 splines = new ArrayList<Cubic>();
031 }
032
033 /* (non-Javadoc)
034 * @see edu.nrao.sss.math.AbsInterp#getComputedValueFor(double)
035 */
036 @Override
037 public double getComputedValueFor(double x)
038 {
039 return getSplineFor(x).evaluate(x);
040 }
041
042 /** Returns the appropriate spline for x. */
043 private Cubic getSplineFor(double x)
044 {
045 //This code assumes:
046 // 1. x is in domain
047 // 2. splines are sorted
048
049 Cubic spline = null;
050
051 for (int s=splines.size()-1; s >= 0; s--)
052 {
053 Cubic candidate = splines.get(s);
054
055 if (x >= candidate.X)
056 {
057 spline = candidate;
058 break;
059 }
060 }
061
062 return spline;
063 }
064
065 /* (non-Javadoc)
066 * @see edu.nrao.sss.math.AbsInterp#compute()
067 */
068 @Override
069 void runAlgorithm()
070 {
071 //TODO Think about also coding clamped spline. In order to calc a clamped
072 // spline, we need the value of the derivate at x[0] and x[n].
073 // If the points are evenly spaced, there are ways to estimate
074 // these values.
075 computeNaturalSplines();
076 }
077
078 /**
079 * Calculates the natural splines for this interpolator.
080 * "Natural" means the second derivatives at the endpoints are taken to be 0.
081 *
082 * This algorithm was take from section 3.6 of <u>Numerical Analysis, Fourth
083 * Edition</u> by Burden & Faires, p 131.
084 * The variable names in this method are, as much as possible, the same as
085 * those in the book.
086 */
087 private void computeNaturalSplines()
088 {
089 int n = this.size() - 1;
090
091 double[] a = new double[n+1]; //term of order 0, also = f(x[i])
092 double[] b = new double[n]; //term of order 1
093 double[] c = new double[n+1]; //term of order 2 (see step 5 for why n+1)
094 double[] d = new double[n]; //term of order 3
095 double[] x = new double[n+1];
096
097 //INPUT
098 for (int i=0; i < n+1; i++)
099 {
100 Point2D point = getPoint(i);
101
102 x[i] = point.getX();
103 a[i] = point.getY();
104 }
105
106 //STEP 1
107 double[] h = new double[n];
108 for (int i=0; i < n; i++)
109 h[i] = x[i+1] - x[i];
110
111 //STEP 2
112 double[] alpha = new double[n];
113 alpha[0] = 0.0;
114 for (int i=1; i < n; i++)
115 {
116 alpha[i] =
117 3 * (a[i+1]*h[i-1] - a[i]*(x[i+1]-x[i-1]) + a[i-1]*h[i]) / (h[i-1]*h[i]);
118 }
119
120 //STEP 3
121 double[] l = new double[n+1];
122 double[] mu = new double[n+1];
123 double[] z = new double[n+1];
124
125 l[0] = 1.0; mu[0] = 0.0; z[0] = 0.0;
126
127 //STEP 4
128 for (int i=1; i < n; i++)
129 {
130 l[i] = 2*(x[i+1]-x[i-1]) - h[i-1]*mu[i-1];
131 mu[i] = h[i] / l[i];
132 z[i] = (alpha[i] - h[i-1]*z[i-1]) / l[i];
133 }
134
135 //STEP 5
136 l[n] = 1.0; z[n] = 0.0; c[n] = 0;
137
138 //STEP 6
139 for (int j=n-1; j >= 0; j--)
140 {
141 c[j] = z[j] - mu[j]*c[j+1];
142 b[j] = (a[j+1]-a[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3.0;
143 d[j] = (c[j+1]-c[j])/(3.0*h[j]);
144 }
145
146 //STEP 7
147 splines.clear();
148 for (int j=0; j < n; j++)
149 {
150 Cubic spline = new Cubic();
151
152 spline.X = x[j];
153 spline.a = a[j];
154 spline.b = b[j];
155 spline.c = c[j];
156 spline.d = d[j];
157
158 splines.add(spline);
159 }
160 }
161
162 /**
163 * A single cubic equation of form
164 * f(x) = a + b*(x-X) + c*(x-X)^2 + d*(x-X)^3.
165 */
166 private class Cubic
167 {
168 double a, b, c, d, X;
169
170 double evaluate(double x)
171 {
172 double delta = x - X;
173
174 return a + b * delta + c * delta * delta + d * delta * delta * delta;
175 }
176 }
177 }