import java.util.*; import java.awt.*; // a cubic spline to be drawn through a vector of VPoint's // // Author: Alexander Bogomolny, CTK Software, Inc. // URL: http://www.cut-the-knot.com // Date: November 30, 2000 // Copyright: A. Bogomolny // Permission to use and modify the file is therefore granted // as long as this comment remains unchanged. Do this at your // own risk. // class CubicSpline { int N; Vector P; // N+1 points to interpolate through (screen coordinates) public void SetMoveables(Vector mp) { N = mp.size()-1; P = mp; } public Vector GetMoveables() { return P; } double[] M; // (N-1) second derivatives at "inner" points double[] In; // N integrals over N intervals private int I = 0; // points to the current interval where value is computed boolean Shade = false; public void SetShade(boolean set) { Shade = set; } public CubicSpline(Vector mp) { try { N = mp.size()-1; P = mp; } catch (Exception e) { N = 0; P = new Vector(); } } // solves for the vector of // second derivatives // public void Resolve() { if (N < 2) return; double[] h = new double[N]; for (int i = 0; i < N; i++) h[i] = ((VPoint)P.elementAt(i+1)).ScaledX() - ((VPoint)P.elementAt(i)).ScaledX(); double[] d = new double[N-1]; for (int i = 0; i < N-1; i++) d[i] = (h[i] + h[i+1])/3; double[] m = new double[N-2]; double[] p = new double[N-2]; for (int i = 0; i < N-2; i++) m[i] = p[i] = h[i+1]/6; Tridiagonal tr = new Tridiagonal(N-1, m, d, p); double[] dy = new double[N]; for (int i = 0; i < N; i++) dy[i] = ((VPoint)P.elementAt(i+1)).ScaledY() - ((VPoint)P.elementAt(i)).ScaledY(); double[] r = new double[N-1]; for (int i = 0; i < N-1; i++) r[i] = dy[i+1]/h[i+1] - dy[i]/h[i]; M = tr.Solve(r); } public double Value(double x) { I = FindInterval(x); return VPoint.ReverseY(Compute(I, VPoint.ScaledX(x))); } private int FindInterval(double x) { int Lo = 0; int Hi = N; do { int i = (Lo + Hi)/2; double X = ((VPoint)P.elementAt(i)).x; if (X > x) Hi = i; else Lo = i; } while (Hi - Lo > 1); return Lo; } // x comes in screen coordinates // private double Compute(int index, double x) { double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX(); double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h; double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h; double yA = ((VPoint)P.elementAt(index)).ScaledY(); double yB = ((VPoint)P.elementAt(index+1)).ScaledY(); double mA = index == 0 ? 0 : M[index-1]; double mB = index > N-2 ? 0 : M[index]; return A*yA + B*yB + ((A*A*A - A)*mA + (B*B*B - B)*mB)*h*h/6; } public double Derivative(double x) { I = FindInterval(x); return VPoint.ReverseY(ComputeDerivative(I, VPoint.ScaledX(x))); } private double ComputeDerivative(int index, double x) { double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX(); double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h; double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h; double yA = ((VPoint)P.elementAt(index)).ScaledY(); double yB = ((VPoint)P.elementAt(index+1)).ScaledY(); double mA = index == 0 ? 0 : M[index-1]; double mB = index > N-2 ? 0 : M[index]; return (yB - yA)/h - (mB - mA)*h/6 + (mB*B*B - mA*A*A)*h/2; } public void Draw(Graphics g, int xFrom, int xTo, int dx) { int x = xFrom; Point p = new Point(x, (int)Math.round(Value(x))); int y = (int)Math.round(VPoint.ReverseY(0)); if (Shade) g.drawLine(x, y, x, p.y); while (x < xTo) { x += dx; Point q = new Point(x, (int)Math.round(Value(x))); g.drawLine(p.x, p.y, q.x, q.y); if (Shade) g.drawLine(x, y, x, q.y); p = q; } } public void DrawDerivative(Graphics g, int xFrom, int xTo, int dx) { int x = xFrom; Point p = new Point(x, (int)Math.round(Derivative(x))); while (x < xTo) { x += dx; Point q = new Point(x, (int)Math.round(Derivative(x))); g.drawLine(p.x, p.y, q.x, q.y); p = q; } } // saves time when computing integrals // for a variable x // public void FillIntegralArray() { double[] h = new double[N]; for (int i = 0; i < N; i++) h[i] = ((VPoint)P.elementAt(i+1)).ScaledX() - ((VPoint)P.elementAt(i)).ScaledX(); In = new double[N]; for (int i = 0; i < N; i++) { double mA = i == 0 ? 0 : M[i-1]; double mB = i > N-2 ? 0 : M[i]; double yA = ((VPoint)P.elementAt(i)).ScaledY(); double yB = ((VPoint)P.elementAt(i+1)).ScaledY(); double G = h[i]*(yB+yA)/2 - h[i]*h[i]*h[i]*(mB+mA)/24; In[i] = i > 0 ? In[i-1] + G : G; } } // definite integral from the first // data point through x // public double Integral(double x) { I = FindInterval(x); return VPoint.ReverseY(ComputeIntegral(I, VPoint.ScaledX(x))); } private double ComputeIntegral(int index, double x) { if (index >= N) return In[N-1]; double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX(); double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h; double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h; double A2 = A*A; double B2 = B*B; double mA = index == 0 ? 0 : M[index-1]; double mB = index > N-2 ? 0 : M[index]; double yA = ((VPoint)P.elementAt(index)).ScaledY(); double yB = ((VPoint)P.elementAt(index+1)).ScaledY(); double f = (index > 0 ? In[index-1] : 0) + (mB*B2*B2 - mA*A2*A2)*h*h*h/24 + (yB*B2 - yA*A2)*h/2 - (mB*B2 - mA*A2)*h*h*h/12 + h*yA/2 - h*h*h*mA/24; return f; } public void DrawIntegral(Graphics g, int xFrom, int xTo, int dx) { int x = xFrom; Point p = new Point(x, (int)Math.round(Integral(x))); while (x < xTo) { x += dx; Point q = new Point(x, (int)Math.round(Integral(x))); g.drawLine(p.x, p.y, q.x, q.y); p = q; } } // Checks whether a mouse press has occured // close to one of the interpolation points. // If not creates such a point and adds it // to the array // public Point CloseEnough(int x, int y) { int i = FindInterval(x); double Y = VPoint.ReverseY(Compute(i, VPoint.ScaledX(x))); if (Math.abs(Y - y) < 5) { VPoint p = new VPoint(x, y); p.SetRadius(2); P.insertElementAt(p, i+1); N++; return p; } return null; } public double SecondDerivatve(double x) { I = FindInterval(x); return VPoint.ReverseY(ComputeSecondDerivative(I, VPoint.ScaledX(x))); } private double ComputeSecondDerivative(int index, double x) { double h = ((VPoint)P.elementAt(index+1)).ScaledX() - ((VPoint)P.elementAt(index)).ScaledX(); double A = (((VPoint)P.elementAt(index+1)).ScaledX() - x)/h; double B = (x - ((VPoint)P.elementAt(index)).ScaledX())/h; double mA = index == 0 ? 0 : M[index-1]; double mB = index > N-2 ? 0 : M[index]; return (A*mA + B*mB)*h; } }