spline3.cs

三样条插值函数算法

        while( l!=r-1 )
        {
            m = (l+r)/2;
            if( c[m]>=x )
            {
                r = m;
            }
            else
            {
                l = m;
            }
        }
        
        //
        // Interpolation
        //
        x = x-c[l];
        m = 3+n+4*(l-3);
        result = c[m]+x*(c[m+1]+x*(c[m+2]+x*c[m+3]));
        return result;
    }


    /*************************************************************************
    This subroutine differentiates the spline.

    Input parameters:
        C   -   coefficients table. Built by BuildLinearSpline,
                BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
        X   -   point

    Result:
        S   -   S(x)
        DS  -   S'(x)
        D2S -   S''(x)

      -- ALGLIB PROJECT --
         Copyright 24.06.2007 by Bochkanov Sergey
    *************************************************************************/
    public static void splinedifferentiation(ref double[] c,
        double x,
        ref double s,
        ref double ds,
        ref double d2s)
    {
        int n = 0;
        int l = 0;
        int r = 0;
        int m = 0;

        System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
        n = (int)Math.Round(c[2]);
        
        //
        // Binary search
        //
        l = 3;
        r = 3+n-2+1;
        while( l!=r-1 )
        {
            m = (l+r)/2;
            if( c[m]>=x )
            {
                r = m;
            }
            else
            {
                l = m;
            }
        }
        
        //
        // Differentiation
        //
        x = x-c[l];
        m = 3+n+4*(l-3);
        s = c[m]+x*(c[m+1]+x*(c[m+2]+x*c[m+3]));
        ds = c[m+1]+2*x*c[m+2]+3*AP.Math.Sqr(x)*c[m+3];
        d2s = 2*c[m+2]+6*x*c[m+3];
    }


    /*************************************************************************
    This subroutine makes the copy of the spline.

    Input parameters:
        C   -   coefficients table. Built by BuildLinearSpline,
                BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.

    Result:
        CC  -   spline copy

      -- ALGLIB PROJECT --
         Copyright 29.06.2007 by Bochkanov Sergey
    *************************************************************************/
    public static void splinecopy(ref double[] c,
        ref double[] cc)
    {
        int s = 0;
        int i_ = 0;

        s = (int)Math.Round(c[0]);
        cc = new double[s-1+1];
        for(i_=0; i_<=s-1;i_++)
        {
            cc[i_] = c[i_];
        }
    }


    /*************************************************************************
    This subroutine unpacks the spline into the coefficients table.

    Input parameters:
        C   -   coefficients table. Built by BuildLinearSpline,
                BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
        X   -   point

    Result:
        Tbl -   coefficients table, unpacked format, array[0..N-2, 0..5].
                For I = 0...N-2:
                    Tbl[I,0] = X[i]
                    Tbl[I,1] = X[i+1]
                    Tbl[I,2] = C0
                    Tbl[I,3] = C1
                    Tbl[I,4] = C2
                    Tbl[I,5] = C3
                On [x[i], x[i+1]] spline is equals to:
                    S(x) = C0 + C1*t + C2*t^2 + C3*t^3
                    t = x-x[i]

      -- ALGLIB PROJECT --
         Copyright 29.06.2007 by Bochkanov Sergey
    *************************************************************************/
    public static void splineunpack(ref double[] c,
        ref int n,
        ref double[,] tbl)
    {
        int i = 0;

        System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineUnpack: incorrect C!");
        n = (int)Math.Round(c[2]);
        tbl = new double[n-2+1, 5+1];
        
        //
        // Fill
        //
        for(i=0; i<=n-2; i++)
        {
            tbl[i,0] = c[3+i];
            tbl[i,1] = c[3+i+1];
            tbl[i,2] = c[3+n+4*i];
            tbl[i,3] = c[3+n+4*i+1];
            tbl[i,4] = c[3+n+4*i+2];
            tbl[i,5] = c[3+n+4*i+3];
        }
    }


    /*************************************************************************
    This subroutine performs linear transformation of the spline argument.

    Input parameters:
        C   -   coefficients table. Built by BuildLinearSpline,
                BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
        A, B-   transformation coefficients: x = A*t + B
    Result:
        C   -   transformed spline

      -- ALGLIB PROJECT --
         Copyright 30.06.2007 by Bochkanov Sergey
    *************************************************************************/
    public static void splinelintransx(ref double[] c,
        double a,
        double b)
    {
        int i = 0;
        int n = 0;
        double v = 0;
        double dv = 0;
        double d2v = 0;
        double[] x = new double[0];
        double[] y = new double[0];
        double[] d = new double[0];

        System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
        n = (int)Math.Round(c[2]);
        
        //
        // Special case: A=0
        //
        if( a==0 )
        {
            v = splineinterpolation(ref c, b);
            for(i=0; i<=n-2; i++)
            {
                c[3+n+4*i] = v;
                c[3+n+4*i+1] = 0;
                c[3+n+4*i+2] = 0;
                c[3+n+4*i+3] = 0;
            }
            return;
        }
        
        //
        // General case: A<>0.
        // Unpack, X, Y, dY/dX.
        // Scale and pack again.
        //
        x = new double[n-1+1];
        y = new double[n-1+1];
        d = new double[n-1+1];
        for(i=0; i<=n-1; i++)
        {
            x[i] = c[3+i];
            splinedifferentiation(ref c, x[i], ref v, ref dv, ref d2v);
            x[i] = (x[i]-b)/a;
            y[i] = v;
            d[i] = a*dv;
        }
        buildhermitespline(x, y, d, n, ref c);
    }


    /*************************************************************************
    This subroutine performs linear transformation of the spline.

    Input parameters:
        C   -   coefficients table. Built by BuildLinearSpline,
                BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
        A, B-   transformation coefficients: S2(x) = A*S(x) + B
    Result:
        C   -   transformed spline

      -- ALGLIB PROJECT --
         Copyright 30.06.2007 by Bochkanov Sergey
    *************************************************************************/
    public static void splinelintransy(ref double[] c,
        double a,
        double b)
    {
        int i = 0;
        int n = 0;
        double v = 0;
        double dv = 0;
        double d2v = 0;
        double[] x = new double[0];
        double[] y = new double[0];
        double[] d = new double[0];

        System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
        n = (int)Math.Round(c[2]);
        
        //
        // Special case: A=0
        //
        for(i=0; i<=n-2; i++)
        {
            c[3+n+4*i] = a*c[3+n+4*i]+b;
            c[3+n+4*i+1] = a*c[3+n+4*i+1];
            c[3+n+4*i+2] = a*c[3+n+4*i+2];
            c[3+n+4*i+3] = a*c[3+n+4*i+3];
        }
    }


    /*************************************************************************
    This subroutine integrates the spline.

    Input parameters:
        C   -   coefficients table. Built by BuildLinearSpline,
                BuildHermiteSpline, BuildCubicSpline, BuildAkimaSpline.
        X   -   right bound of the integration interval [a, x]
    Result:
        integral(S(t)dt,a,x)

      -- ALGLIB PROJECT --
         Copyright 23.06.2007 by Bochkanov Sergey
    *************************************************************************/
    public static double splineintegration(ref double[] c,
        double x)
    {
        double result = 0;
        int n = 0;
        int i = 0;
        int l = 0;
        int r = 0;
        int m = 0;
        double w = 0;

        System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineIntegration: incorrect C!");
        n = (int)Math.Round(c[2]);
        
        //
        // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
        //
        l = 3;
        r = 3+n-2+1;
        while( l!=r-1 )
        {
            m = (l+r)/2;
            if( c[m]>=x )
            {
                r = m;
            }
            else
            {
                l = m;
            }
        }
        
        //
        // Integration
        //
        result = 0;
        for(i=3; i<=l-1; i++)
        {
            w = c[i+1]-c[i];
            m = 3+n+4*(i-3);
            result = result+c[m]*w;
            result = result+c[m+1]*AP.Math.Sqr(w)/2;
            result = result+c[m+2]*AP.Math.Sqr(w)*w/3;
            result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
        }
        w = x-c[l];
        m = 3+n+4*(l-3);
        result = result+c[m]*w;
        result = result+c[m+1]*AP.Math.Sqr(w)/2;
        result = result+c[m+2]*AP.Math.Sqr(w)*w/3;
        result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
        return result;
    }


    /*************************************************************************
    Obsolete subroutine, left for backward compatibility.
    *************************************************************************/
    public static void spline3buildtable(int n,
        int diffn,
        double[] x,
        double[] y,
        double boundl,
        double boundr,
        ref double[,] ctbl)
    {
        bool c = new bool();
        int e = 0;
        int g = 0;
        double tmp = 0;
        int nxm1 = 0;
        int i = 0;
        int j = 0;
        double dx = 0;
        double dxj = 0;
        double dyj = 0;
        double dxjp1 = 0;
        double dyjp1 = 0;
        double dxp = 0;
        double dyp = 0;
        double yppa = 0;
        double yppb = 0;
        double pj = 0;
        double b1 = 0;
        double b2 = 0;
        double b3 = 0;
        double b4 = 0;

        x = (double[])x.Clone();
        y = (double[])y.Clone();

        n = n-1;
        g = (n+1)/2;
        do
        {
            i = g;
            do
            {
                j = i-g;
                c = true;
                do
                {
                    if( x[j]<=x[j+g] )
                    {
                        c = false;
                    }
                    else
                    {
                        tmp = x[j];
                        x[j] = x[j+g];
                        x[j+g] = tmp;
                        tmp = y[j];
                        y[j] = y[j+g];
                        y[j+g] = tmp;
                    }
                    j = j-1;
                }
                while( j>=0 & c );
                i = i+1;
            }
            while( i<=n );
            g = g/2;
        }
        while( g>0 );
        ctbl = new double[4+1, n+1];
        n = n+1;
        if( diffn==1 )
        {
            b1 = 1;
            b2 = 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-boundl);
            b3 = 1;
            b4 = 6/(x[n-1]-x[n-2])*(boundr-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
        }
        else
        {
            b1 = 0;
            b2 = 2*boundl;
            b3 = 0;
            b4 = 2*boundr;
        }
        nxm1 = n-1;
        if( n>=2 )
        {
            if( n>2 )
            {
                dxj = x[1]-x[0];
                dyj = y[1]-y[0];
                j = 2;
                while( j<=nxm1 )