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;
    }
}