source: trunk/sources/ALGLIB/spline2d.cs @ 2643

/*************************************************************************
Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

using System;

namespace alglib
{
    public class spline2d
    {
        /*************************************************************************
        2-dimensional spline inteprolant
        *************************************************************************/
        public struct spline2dinterpolant
        {
            public int k;
            public double[] c;
        };




        public const int spline2dvnum = 12;


        /*************************************************************************
        This subroutine builds bilinear spline coefficients table.

        Input parameters:
            X   -   spline abscissas, array[0..N-1]
            Y   -   spline ordinates, array[0..M-1]
            F   -   function values, array[0..M-1,0..N-1]
            M,N -   grid size, M>=2, N>=2

        Output parameters:
            C   -   spline interpolant

          -- ALGLIB PROJECT --
             Copyright 05.07.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dbuildbilinear(double[] x,
            double[] y,
            double[,] f,
            int m,
            int n,
            ref spline2dinterpolant c)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int tblsize = 0;
            int shift = 0;
            double t = 0;
            double[,] dx = new double[0,0];
            double[,] dy = new double[0,0];
            double[,] dxy = new double[0,0];

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

            System.Diagnostics.Debug.Assert(n>=2 & m>=2, "Spline2DBuildBilinear: N<2 or M<2!");
            
            //
            // Sort points
            //
            for(j=0; j<=n-1; j++)
            {
                k = j;
                for(i=j+1; i<=n-1; i++)
                {
                    if( (double)(x[i])<(double)(x[k]) )
                    {
                        k = i;
                    }
                }
                if( k!=j )
                {
                    for(i=0; i<=m-1; i++)
                    {
                        t = f[i,j];
                        f[i,j] = f[i,k];
                        f[i,k] = t;
                    }
                    t = x[j];
                    x[j] = x[k];
                    x[k] = t;
                }
            }
            for(i=0; i<=m-1; i++)
            {
                k = i;
                for(j=i+1; j<=m-1; j++)
                {
                    if( (double)(y[j])<(double)(y[k]) )
                    {
                        k = j;
                    }
                }
                if( k!=i )
                {
                    for(j=0; j<=n-1; j++)
                    {
                        t = f[i,j];
                        f[i,j] = f[k,j];
                        f[k,j] = t;
                    }
                    t = y[i];
                    y[i] = y[k];
                    y[k] = t;
                }
            }
            
            //
            // Fill C:
            //  C[0]            -   length(C)
            //  C[1]            -   type(C):
            //                      -1 = bilinear interpolant
            //                      -3 = general cubic spline
            //                           (see BuildBicubicSpline)
            //  C[2]:
            //      N (x count)
            //  C[3]:
            //      M (y count)
            //  C[4]...C[4+N-1]:
            //      x[i], i = 0...N-1
            //  C[4+N]...C[4+N+M-1]:
            //      y[i], i = 0...M-1
            //  C[4+N+M]...C[4+N+M+(N*M-1)]:
            //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
            //
            c.k = 1;
            tblsize = 4+n+m+n*m;
            c.c = new double[tblsize-1+1];
            c.c[0] = tblsize;
            c.c[1] = -1;
            c.c[2] = n;
            c.c[3] = m;
            for(i=0; i<=n-1; i++)
            {
                c.c[4+i] = x[i];
            }
            for(i=0; i<=m-1; i++)
            {
                c.c[4+n+i] = y[i];
            }
            for(i=0; i<=m-1; i++)
            {
                for(j=0; j<=n-1; j++)
                {
                    shift = i*n+j;
                    c.c[4+n+m+shift] = f[i,j];
                }
            }
        }


        /*************************************************************************
        This subroutine builds bicubic spline coefficients table.

        Input parameters:
            X   -   spline abscissas, array[0..N-1]
            Y   -   spline ordinates, array[0..M-1]
            F   -   function values, array[0..M-1,0..N-1]
            M,N -   grid size, M>=2, N>=2

        Output parameters:
            C   -   spline interpolant

          -- ALGLIB PROJECT --
             Copyright 05.07.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dbuildbicubic(double[] x,
            double[] y,
            double[,] f,
            int m,
            int n,
            ref spline2dinterpolant c)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int tblsize = 0;
            int shift = 0;
            double t = 0;
            double[,] dx = new double[0,0];
            double[,] dy = new double[0,0];
            double[,] dxy = new double[0,0];

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

            System.Diagnostics.Debug.Assert(n>=2 & m>=2, "BuildBicubicSpline: N<2 or M<2!");
            
            //
            // Sort points
            //
            for(j=0; j<=n-1; j++)
            {
                k = j;
                for(i=j+1; i<=n-1; i++)
                {
                    if( (double)(x[i])<(double)(x[k]) )
                    {
                        k = i;
                    }
                }
                if( k!=j )
                {
                    for(i=0; i<=m-1; i++)
                    {
                        t = f[i,j];
                        f[i,j] = f[i,k];
                        f[i,k] = t;
                    }
                    t = x[j];
                    x[j] = x[k];
                    x[k] = t;
                }
            }
            for(i=0; i<=m-1; i++)
            {
                k = i;
                for(j=i+1; j<=m-1; j++)
                {
                    if( (double)(y[j])<(double)(y[k]) )
                    {
                        k = j;
                    }
                }
                if( k!=i )
                {
                    for(j=0; j<=n-1; j++)
                    {
                        t = f[i,j];
                        f[i,j] = f[k,j];
                        f[k,j] = t;
                    }
                    t = y[i];
                    y[i] = y[k];
                    y[k] = t;
                }
            }
            
            //
            // Fill C:
            //  C[0]            -   length(C)
            //  C[1]            -   type(C):
            //                      -1 = bilinear interpolant
            //                           (see BuildBilinearInterpolant)
            //                      -3 = general cubic spline
            //  C[2]:
            //      N (x count)
            //  C[3]:
            //      M (y count)
            //  C[4]...C[4+N-1]:
            //      x[i], i = 0...N-1
            //  C[4+N]...C[4+N+M-1]:
            //      y[i], i = 0...M-1
            //  C[4+N+M]...C[4+N+M+(N*M-1)]:
            //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
            //  C[4+N+M+N*M]...C[4+N+M+(2*N*M-1)]:
            //      df(i,j)/dx table.
            //  C[4+N+M+2*N*M]...C[4+N+M+(3*N*M-1)]:
            //      df(i,j)/dy table.
            //  C[4+N+M+3*N*M]...C[4+N+M+(4*N*M-1)]:
            //      d2f(i,j)/dxdy table.
            //
            c.k = 3;
            tblsize = 4+n+m+4*n*m;
            c.c = new double[tblsize-1+1];
            c.c[0] = tblsize;
            c.c[1] = -3;
            c.c[2] = n;
            c.c[3] = m;
            for(i=0; i<=n-1; i++)
            {
                c.c[4+i] = x[i];
            }
            for(i=0; i<=m-1; i++)
            {
                c.c[4+n+i] = y[i];
            }
            bicubiccalcderivatives(ref f, ref x, ref y, m, n, ref dx, ref dy, ref dxy);
            for(i=0; i<=m-1; i++)
            {
                for(j=0; j<=n-1; j++)
                {
                    shift = i*n+j;
                    c.c[4+n+m+shift] = f[i,j];
                    c.c[4+n+m+n*m+shift] = dx[i,j];
                    c.c[4+n+m+2*n*m+shift] = dy[i,j];
                    c.c[4+n+m+3*n*m+shift] = dxy[i,j];
                }
            }
        }


        /*************************************************************************
        This subroutine calculates the value of the bilinear or bicubic spline  at
        the given point X.

        Input parameters:
            C   -   coefficients table.
                    Built by BuildBilinearSpline or BuildBicubicSpline.
            X, Y-   point

        Result:
            S(x,y)

          -- ALGLIB PROJECT --
             Copyright 05.07.2007 by Bochkanov Sergey
        *************************************************************************/
        public static double spline2dcalc(ref spline2dinterpolant c,
            double x,
            double y)
        {
            double result = 0;
            double v = 0;
            double vx = 0;
            double vy = 0;
            double vxy = 0;

            spline2ddiff(ref c, x, y, ref v, ref vx, ref vy, ref vxy);
            result = v;
            return result;
        }


        /*************************************************************************
        This subroutine calculates the value of the bilinear or bicubic spline  at
        the given point X and its derivatives.

        Input parameters:
            C   -   spline interpolant.
            X, Y-   point

        Output parameters:
            F   -   S(x,y)
            FX  -   dS(x,y)/dX
            FY  -   dS(x,y)/dY
            FXY -   d2S(x,y)/dXdY

          -- ALGLIB PROJECT --
             Copyright 05.07.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2ddiff(ref spline2dinterpolant c,
            double x,
            double y,
            ref double f,
            ref double fx,
            ref double fy,
            ref double fxy)
        {
            int n = 0;
            int m = 0;
            double t = 0;
            double dt = 0;
            double u = 0;
            double du = 0;
            int i = 0;
            int j = 0;
            int ix = 0;
            int iy = 0;
            int l = 0;
            int r = 0;
            int h = 0;
            int shift1 = 0;
            int s1 = 0;
            int s2 = 0;
            int s3 = 0;
            int s4 = 0;
            int sf = 0;
            int sfx = 0;
            int sfy = 0;
            int sfxy = 0;
            double y1 = 0;
            double y2 = 0;
            double y3 = 0;
            double y4 = 0;
            double v = 0;
            double t0 = 0;
            double t1 = 0;
            double t2 = 0;
            double t3 = 0;
            double u0 = 0;
            double u1 = 0;
            double u2 = 0;
            double u3 = 0;

            System.Diagnostics.Debug.Assert((int)Math.Round(c.c[1])==-1 | (int)Math.Round(c.c[1])==-3, "Spline2DDiff: incorrect C!");
            n = (int)Math.Round(c.c[2]);
            m = (int)Math.Round(c.c[3]);
            
            //
            // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
            //
            l = 4;
            r = 4+n-2+1;
            while( l!=r-1 )
            {
                h = (l+r)/2;
                if( (double)(c.c[h])>=(double)(x) )
                {
                    r = h;
                }
                else
                {
                    l = h;
                }
            }
            t = (x-c.c[l])/(c.c[l+1]-c.c[l]);
            dt = 1.0/(c.c[l+1]-c.c[l]);
            ix = l-4;
            
            //
            // Binary search in the [ y[0], ..., y[m-2] ] (y[m-1] is not included)
            //
            l = 4+n;
            r = 4+n+(m-2)+1;
            while( l!=r-1 )
            {
                h = (l+r)/2;
                if( (double)(c.c[h])>=(double)(y) )
                {
                    r = h;
                }
                else
                {
                    l = h;
                }
            }
            u = (y-c.c[l])/(c.c[l+1]-c.c[l]);
            du = 1.0/(c.c[l+1]-c.c[l]);
            iy = l-(4+n);
            
            //
            // Prepare F, dF/dX, dF/dY, d2F/dXdY
            //
            f = 0;
            fx = 0;
            fy = 0;
            fxy = 0;
            
            //
            // Bilinear interpolation
            //
            if( (int)Math.Round(c.c[1])==-1 )
            {
                shift1 = 4+n+m;
                y1 = c.c[shift1+n*iy+ix];
                y2 = c.c[shift1+n*iy+(ix+1)];
                y3 = c.c[shift1+n*(iy+1)+(ix+1)];
                y4 = c.c[shift1+n*(iy+1)+ix];
                f = (1-t)*(1-u)*y1+t*(1-u)*y2+t*u*y3+(1-t)*u*y4;
                fx = (-((1-u)*y1)+(1-u)*y2+u*y3-u*y4)*dt;
                fy = (-((1-t)*y1)-t*y2+t*y3+(1-t)*y4)*du;
                fxy = (y1-y2+y3-y4)*du*dt;
                return;
            }
            
            //
            // Bicubic interpolation
            //
            if( (int)Math.Round(c.c[1])==-3 )
            {
                
                //
                // Prepare info
                //
                t0 = 1;
                t1 = t;
                t2 = AP.Math.Sqr(t);
                t3 = t*t2;
                u0 = 1;
                u1 = u;
                u2 = AP.Math.Sqr(u);
                u3 = u*u2;
                sf = 4+n+m;
                sfx = 4+n+m+n*m;
                sfy = 4+n+m+2*n*m;
                sfxy = 4+n+m+3*n*m;
                s1 = n*iy+ix;
                s2 = n*iy+(ix+1);
                s3 = n*(iy+1)+(ix+1);
                s4 = n*(iy+1)+ix;
                
                //
                // Calculate
                //
                v = +(1*c.c[sf+s1]);
                f = f+v*t0*u0;
                v = +(1*c.c[sfy+s1]/du);
                f = f+v*t0*u1;
                fy = fy+1*v*t0*u0*du;
                v = -(3*c.c[sf+s1])+3*c.c[sf+s4]-2*c.c[sfy+s1]/du-1*c.c[sfy+s4]/du;
                f = f+v*t0*u2;
                fy = fy+2*v*t0*u1*du;
                v = +(2*c.c[sf+s1])-2*c.c[sf+s4]+1*c.c[sfy+s1]/du+1*c.c[sfy+s4]/du;
                f = f+v*t0*u3;
                fy = fy+3*v*t0*u2*du;
                v = +(1*c.c[sfx+s1]/dt);
                f = f+v*t1*u0;
                fx = fx+1*v*t0*u0*dt;
                v = +(1*c.c[sfxy+s1]/(dt*du));
                f = f+v*t1*u1;
                fx = fx+1*v*t0*u1*dt;
                fy = fy+1*v*t1*u0*du;
                fxy = fxy+1*v*t0*u0*dt*du;
                v = -(3*c.c[sfx+s1]/dt)+3*c.c[sfx+s4]/dt-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
                f = f+v*t1*u2;
                fx = fx+1*v*t0*u2*dt;
                fy = fy+2*v*t1*u1*du;
                fxy = fxy+2*v*t0*u1*dt*du;
                v = +(2*c.c[sfx+s1]/dt)-2*c.c[sfx+s4]/dt+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
                f = f+v*t1*u3;
                fx = fx+1*v*t0*u3*dt;
                fy = fy+3*v*t1*u2*du;
                fxy = fxy+3*v*t0*u2*dt*du;
                v = -(3*c.c[sf+s1])+3*c.c[sf+s2]-2*c.c[sfx+s1]/dt-1*c.c[sfx+s2]/dt;
                f = f+v*t2*u0;
                fx = fx+2*v*t1*u0*dt;
                v = -(3*c.c[sfy+s1]/du)+3*c.c[sfy+s2]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du);
                f = f+v*t2*u1;
                fx = fx+2*v*t1*u1*dt;
                fy = fy+1*v*t2*u0*du;
                fxy = fxy+2*v*t1*u0*dt*du;
                v = +(9*c.c[sf+s1])-9*c.c[sf+s2]+9*c.c[sf+s3]-9*c.c[sf+s4]+6*c.c[sfx+s1]/dt+3*c.c[sfx+s2]/dt-3*c.c[sfx+s3]/dt-6*c.c[sfx+s4]/dt+6*c.c[sfy+s1]/du-6*c.c[sfy+s2]/du-3*c.c[sfy+s3]/du+3*c.c[sfy+s4]/du+4*c.c[sfxy+s1]/(dt*du)+2*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+2*c.c[sfxy+s4]/(dt*du);
                f = f+v*t2*u2;
                fx = fx+2*v*t1*u2*dt;
                fy = fy+2*v*t2*u1*du;
                fxy = fxy+4*v*t1*u1*dt*du;
                v = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-4*c.c[sfx+s1]/dt-2*c.c[sfx+s2]/dt+2*c.c[sfx+s3]/dt+4*c.c[sfx+s4]/dt-3*c.c[sfy+s1]/du+3*c.c[sfy+s2]/du+3*c.c[sfy+s3]/du-3*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-2*c.c[sfxy+s4]/(dt*du);
                f = f+v*t2*u3;
                fx = fx+2*v*t1*u3*dt;
                fy = fy+3*v*t2*u2*du;
                fxy = fxy+6*v*t1*u2*dt*du;
                v = +(2*c.c[sf+s1])-2*c.c[sf+s2]+1*c.c[sfx+s1]/dt+1*c.c[sfx+s2]/dt;
                f = f+v*t3*u0;
                fx = fx+3*v*t2*u0*dt;
                v = +(2*c.c[sfy+s1]/du)-2*c.c[sfy+s2]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du);
                f = f+v*t3*u1;
                fx = fx+3*v*t2*u1*dt;
                fy = fy+1*v*t3*u0*du;
                fxy = fxy+3*v*t2*u0*dt*du;
                v = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-3*c.c[sfx+s1]/dt-3*c.c[sfx+s2]/dt+3*c.c[sfx+s3]/dt+3*c.c[sfx+s4]/dt-4*c.c[sfy+s1]/du+4*c.c[sfy+s2]/du+2*c.c[sfy+s3]/du-2*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-2*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
                f = f+v*t3*u2;
                fx = fx+3*v*t2*u2*dt;
                fy = fy+2*v*t3*u1*du;
                fxy = fxy+6*v*t2*u1*dt*du;
                v = +(4*c.c[sf+s1])-4*c.c[sf+s2]+4*c.c[sf+s3]-4*c.c[sf+s4]+2*c.c[sfx+s1]/dt+2*c.c[sfx+s2]/dt-2*c.c[sfx+s3]/dt-2*c.c[sfx+s4]/dt+2*c.c[sfy+s1]/du-2*c.c[sfy+s2]/du-2*c.c[sfy+s3]/du+2*c.c[sfy+s4]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
                f = f+v*t3*u3;
                fx = fx+3*v*t2*u3*dt;
                fy = fy+3*v*t3*u2*du;
                fxy = fxy+9*v*t2*u2*dt*du;
                return;
            }
        }


        /*************************************************************************
        This subroutine unpacks two-dimensional spline into the coefficients table

        Input parameters:
            C   -   spline interpolant.

        Result:
            M, N-   grid size (x-axis and y-axis)
            Tbl -   coefficients table, unpacked format,
                    [0..(N-1)*(M-1)-1, 0..19].
                    For I = 0...M-2, J=0..N-2:
                        K =  I*(N-1)+J
                        Tbl[K,0] = X[j]
                        Tbl[K,1] = X[j+1]
                        Tbl[K,2] = Y[i]
                        Tbl[K,3] = Y[i+1]
                        Tbl[K,4] = C00
                        Tbl[K,5] = C01
                        Tbl[K,6] = C02
                        Tbl[K,7] = C03
                        Tbl[K,8] = C10
                        Tbl[K,9] = C11
                        ...
                        Tbl[K,19] = C33
                    On each grid square spline is equals to:
                        S(x) = SUM(c[i,j]*(x^i)*(y^j), i=0..3, j=0..3)
                        t = x-x[j]
                        u = y-y[i]

          -- ALGLIB PROJECT --
             Copyright 29.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dunpack(ref spline2dinterpolant c,
            ref int m,
            ref int n,
            ref double[,] tbl)
        {
            int i = 0;
            int j = 0;
            int ci = 0;
            int cj = 0;
            int k = 0;
            int p = 0;
            int shift = 0;
            int s1 = 0;
            int s2 = 0;
            int s3 = 0;
            int s4 = 0;
            int sf = 0;
            int sfx = 0;
            int sfy = 0;
            int sfxy = 0;
            double y1 = 0;
            double y2 = 0;
            double y3 = 0;
            double y4 = 0;
            double dt = 0;
            double du = 0;

            System.Diagnostics.Debug.Assert((int)Math.Round(c.c[1])==-3 | (int)Math.Round(c.c[1])==-1, "SplineUnpack2D: incorrect C!");
            n = (int)Math.Round(c.c[2]);
            m = (int)Math.Round(c.c[3]);
            tbl = new double[(n-1)*(m-1)-1+1, 19+1];
            
            //
            // Fill
            //
            for(i=0; i<=m-2; i++)
            {
                for(j=0; j<=n-2; j++)
                {
                    p = i*(n-1)+j;
                    tbl[p,0] = c.c[4+j];
                    tbl[p,1] = c.c[4+j+1];
                    tbl[p,2] = c.c[4+n+i];
                    tbl[p,3] = c.c[4+n+i+1];
                    dt = 1/(tbl[p,1]-tbl[p,0]);
                    du = 1/(tbl[p,3]-tbl[p,2]);
                    
                    //
                    // Bilinear interpolation
                    //
                    if( (int)Math.Round(c.c[1])==-1 )
                    {
                        for(k=4; k<=19; k++)
                        {
                            tbl[p,k] = 0;
                        }
                        shift = 4+n+m;
                        y1 = c.c[shift+n*i+j];
                        y2 = c.c[shift+n*i+(j+1)];
                        y3 = c.c[shift+n*(i+1)+(j+1)];
                        y4 = c.c[shift+n*(i+1)+j];
                        tbl[p,4] = y1;
                        tbl[p,4+1*4+0] = y2-y1;
                        tbl[p,4+0*4+1] = y4-y1;
                        tbl[p,4+1*4+1] = y3-y2-y4+y1;
                    }
                    
                    //
                    // Bicubic interpolation
                    //
                    if( (int)Math.Round(c.c[1])==-3 )
                    {
                        sf = 4+n+m;
                        sfx = 4+n+m+n*m;
                        sfy = 4+n+m+2*n*m;
                        sfxy = 4+n+m+3*n*m;
                        s1 = n*i+j;
                        s2 = n*i+(j+1);
                        s3 = n*(i+1)+(j+1);
                        s4 = n*(i+1)+j;
                        tbl[p,4+0*4+0] = +(1*c.c[sf+s1]);
                        tbl[p,4+0*4+1] = +(1*c.c[sfy+s1]/du);
                        tbl[p,4+0*4+2] = -(3*c.c[sf+s1])+3*c.c[sf+s4]-2*c.c[sfy+s1]/du-1*c.c[sfy+s4]/du;
                        tbl[p,4+0*4+3] = +(2*c.c[sf+s1])-2*c.c[sf+s4]+1*c.c[sfy+s1]/du+1*c.c[sfy+s4]/du;
                        tbl[p,4+1*4+0] = +(1*c.c[sfx+s1]/dt);
                        tbl[p,4+1*4+1] = +(1*c.c[sfxy+s1]/(dt*du));
                        tbl[p,4+1*4+2] = -(3*c.c[sfx+s1]/dt)+3*c.c[sfx+s4]/dt-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
                        tbl[p,4+1*4+3] = +(2*c.c[sfx+s1]/dt)-2*c.c[sfx+s4]/dt+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
                        tbl[p,4+2*4+0] = -(3*c.c[sf+s1])+3*c.c[sf+s2]-2*c.c[sfx+s1]/dt-1*c.c[sfx+s2]/dt;
                        tbl[p,4+2*4+1] = -(3*c.c[sfy+s1]/du)+3*c.c[sfy+s2]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du);
                        tbl[p,4+2*4+2] = +(9*c.c[sf+s1])-9*c.c[sf+s2]+9*c.c[sf+s3]-9*c.c[sf+s4]+6*c.c[sfx+s1]/dt+3*c.c[sfx+s2]/dt-3*c.c[sfx+s3]/dt-6*c.c[sfx+s4]/dt+6*c.c[sfy+s1]/du-6*c.c[sfy+s2]/du-3*c.c[sfy+s3]/du+3*c.c[sfy+s4]/du+4*c.c[sfxy+s1]/(dt*du)+2*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+2*c.c[sfxy+s4]/(dt*du);
                        tbl[p,4+2*4+3] = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-4*c.c[sfx+s1]/dt-2*c.c[sfx+s2]/dt+2*c.c[sfx+s3]/dt+4*c.c[sfx+s4]/dt-3*c.c[sfy+s1]/du+3*c.c[sfy+s2]/du+3*c.c[sfy+s3]/du-3*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-2*c.c[sfxy+s4]/(dt*du);
                        tbl[p,4+3*4+0] = +(2*c.c[sf+s1])-2*c.c[sf+s2]+1*c.c[sfx+s1]/dt+1*c.c[sfx+s2]/dt;
                        tbl[p,4+3*4+1] = +(2*c.c[sfy+s1]/du)-2*c.c[sfy+s2]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du);
                        tbl[p,4+3*4+2] = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-3*c.c[sfx+s1]/dt-3*c.c[sfx+s2]/dt+3*c.c[sfx+s3]/dt+3*c.c[sfx+s4]/dt-4*c.c[sfy+s1]/du+4*c.c[sfy+s2]/du+2*c.c[sfy+s3]/du-2*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-2*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
                        tbl[p,4+3*4+3] = +(4*c.c[sf+s1])-4*c.c[sf+s2]+4*c.c[sf+s3]-4*c.c[sf+s4]+2*c.c[sfx+s1]/dt+2*c.c[sfx+s2]/dt-2*c.c[sfx+s3]/dt-2*c.c[sfx+s4]/dt+2*c.c[sfy+s1]/du-2*c.c[sfy+s2]/du-2*c.c[sfy+s3]/du+2*c.c[sfy+s4]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
                    }
                    
                    //
                    // Rescale Cij
                    //
                    for(ci=0; ci<=3; ci++)
                    {
                        for(cj=0; cj<=3; cj++)
                        {
                            tbl[p,4+ci*4+cj] = tbl[p,4+ci*4+cj]*Math.Pow(dt, ci)*Math.Pow(du, cj);
                        }
                    }
                }
            }
        }


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

        Input parameters:
            C       -   spline interpolant
            AX, BX  -   transformation coefficients: x = A*t + B
            AY, BY  -   transformation coefficients: y = A*u + B
        Result:
            C   -   transformed spline

          -- ALGLIB PROJECT --
             Copyright 30.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dlintransxy(ref spline2dinterpolant c,
            double ax,
            double bx,
            double ay,
            double by)
        {
            int i = 0;
            int j = 0;
            int n = 0;
            int m = 0;
            double v = 0;
            double[] x = new double[0];
            double[] y = new double[0];
            double[,] f = new double[0,0];
            int typec = 0;

            typec = (int)Math.Round(c.c[1]);
            System.Diagnostics.Debug.Assert(typec==-3 | typec==-1, "Spline2DLinTransXY: incorrect C!");
            n = (int)Math.Round(c.c[2]);
            m = (int)Math.Round(c.c[3]);
            x = new double[n-1+1];
            y = new double[m-1+1];
            f = new double[m-1+1, n-1+1];
            for(j=0; j<=n-1; j++)
            {
                x[j] = c.c[4+j];
            }
            for(i=0; i<=m-1; i++)
            {
                y[i] = c.c[4+n+i];
            }
            for(i=0; i<=m-1; i++)
            {
                for(j=0; j<=n-1; j++)
                {
                    f[i,j] = c.c[4+n+m+i*n+j];
                }
            }
            
            //
            // Special case: AX=0 or AY=0
            //
            if( (double)(ax)==(double)(0) )
            {
                for(i=0; i<=m-1; i++)
                {
                    v = spline2dcalc(ref c, bx, y[i]);
                    for(j=0; j<=n-1; j++)
                    {
                        f[i,j] = v;
                    }
                }
                if( typec==-3 )
                {
                    spline2dbuildbicubic(x, y, f, m, n, ref c);
                }
                if( typec==-1 )
                {
                    spline2dbuildbilinear(x, y, f, m, n, ref c);
                }
                ax = 1;
                bx = 0;
            }
            if( (double)(ay)==(double)(0) )
            {
                for(j=0; j<=n-1; j++)
                {
                    v = spline2dcalc(ref c, x[j], by);
                    for(i=0; i<=m-1; i++)
                    {
                        f[i,j] = v;
                    }
                }
                if( typec==-3 )
                {
                    spline2dbuildbicubic(x, y, f, m, n, ref c);
                }
                if( typec==-1 )
                {
                    spline2dbuildbilinear(x, y, f, m, n, ref c);
                }
                ay = 1;
                by = 0;
            }
            
            //
            // General case: AX<>0, AY<>0
            // Unpack, scale and pack again.
            //
            for(j=0; j<=n-1; j++)
            {
                x[j] = (x[j]-bx)/ax;
            }
            for(i=0; i<=m-1; i++)
            {
                y[i] = (y[i]-by)/ay;
            }
            if( typec==-3 )
            {
                spline2dbuildbicubic(x, y, f, m, n, ref c);
            }
            if( typec==-1 )
            {
                spline2dbuildbilinear(x, y, f, m, n, ref c);
            }
        }


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

        Input parameters:
            C   -   spline interpolant.
            A, B-   transformation coefficients: S2(x,y) = A*S(x,y) + B
            
        Output parameters:
            C   -   transformed spline

          -- ALGLIB PROJECT --
             Copyright 30.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dlintransf(ref spline2dinterpolant c,
            double a,
            double b)
        {
            int i = 0;
            int j = 0;
            int n = 0;
            int m = 0;
            double v = 0;
            double[] x = new double[0];
            double[] y = new double[0];
            double[,] f = new double[0,0];
            int typec = 0;

            typec = (int)Math.Round(c.c[1]);
            System.Diagnostics.Debug.Assert(typec==-3 | typec==-1, "Spline2DLinTransXY: incorrect C!");
            n = (int)Math.Round(c.c[2]);
            m = (int)Math.Round(c.c[3]);
            x = new double[n-1+1];
            y = new double[m-1+1];
            f = new double[m-1+1, n-1+1];
            for(j=0; j<=n-1; j++)
            {
                x[j] = c.c[4+j];
            }
            for(i=0; i<=m-1; i++)
            {
                y[i] = c.c[4+n+i];
            }
            for(i=0; i<=m-1; i++)
            {
                for(j=0; j<=n-1; j++)
                {
                    f[i,j] = a*c.c[4+n+m+i*n+j]+b;
                }
            }
            if( typec==-3 )
            {
                spline2dbuildbicubic(x, y, f, m, n, ref c);
            }
            if( typec==-1 )
            {
                spline2dbuildbilinear(x, y, f, m, n, ref c);
            }
        }


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

        Input parameters:
            C   -   spline interpolant

        Output parameters:
            CC  -   spline copy

          -- ALGLIB PROJECT --
             Copyright 29.06.2007 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dcopy(ref spline2dinterpolant c,
            ref spline2dinterpolant cc)
        {
            int n = 0;
            int i_ = 0;

            System.Diagnostics.Debug.Assert(c.k==1 | c.k==3, "Spline2DCopy: incorrect C!");
            cc.k = c.k;
            n = (int)Math.Round(c.c[0]);
            cc.c = new double[n];
            for(i_=0; i_<=n-1;i_++)
            {
                cc.c[i_] = c.c[i_];
            }
        }


        /*************************************************************************
        Serialization of the spline interpolant

        INPUT PARAMETERS:
            B   -   spline interpolant

        OUTPUT PARAMETERS:
            RA      -   array of real numbers which contains interpolant,
                        array[0..RLen-1]
            RLen    -   RA lenght

          -- ALGLIB --
             Copyright 17.08.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dserialize(ref spline2dinterpolant c,
            ref double[] ra,
            ref int ralen)
        {
            int clen = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert(c.k==1 | c.k==3, "Spline2DSerialize: incorrect C!");
            clen = (int)Math.Round(c.c[0]);
            ralen = 3+clen;
            ra = new double[ralen];
            ra[0] = ralen;
            ra[1] = spline2dvnum;
            ra[2] = c.k;
            i1_ = (0) - (3);
            for(i_=3; i_<=3+clen-1;i_++)
            {
                ra[i_] = c.c[i_+i1_];
            }
        }


        /*************************************************************************
        Unserialization of the spline interpolant

        INPUT PARAMETERS:
            RA  -   array of real numbers which contains interpolant,

        OUTPUT PARAMETERS:
            B   -   spline interpolant

          -- ALGLIB --
             Copyright 17.08.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dunserialize(ref double[] ra,
            ref spline2dinterpolant c)
        {
            int clen = 0;
            int i_ = 0;
            int i1_ = 0;

            System.Diagnostics.Debug.Assert((int)Math.Round(ra[1])==spline2dvnum, "Spline2DUnserialize: corrupted array!");
            c.k = (int)Math.Round(ra[2]);
            clen = (int)Math.Round(ra[3]);
            c.c = new double[clen];
            i1_ = (3) - (0);
            for(i_=0; i_<=clen-1;i_++)
            {
                c.c[i_] = ra[i_+i1_];
            }
        }


        /*************************************************************************
        Bicubic spline resampling

        Input parameters:
            A           -   function values at the old grid,
                            array[0..OldHeight-1, 0..OldWidth-1]
            OldHeight   -   old grid height, OldHeight>1
            OldWidth    -   old grid width, OldWidth>1
            NewHeight   -   new grid height, NewHeight>1
            NewWidth    -   new grid width, NewWidth>1
            
        Output parameters:
            B           -   function values at the new grid,
                            array[0..NewHeight-1, 0..NewWidth-1]

          -- ALGLIB routine --
             15 May, 2007
             Copyright by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dresamplebicubic(ref double[,] a,
            int oldheight,
            int oldwidth,
            ref double[,] b,
            int newheight,
            int newwidth)
        {
            double[,] buf = new double[0,0];
            double[] x = new double[0];
            double[] y = new double[0];
            spline1d.spline1dinterpolant c = new spline1d.spline1dinterpolant();
            int i = 0;
            int j = 0;
            int mw = 0;
            int mh = 0;

            System.Diagnostics.Debug.Assert(oldwidth>1 & oldheight>1, "Spline2DResampleBicubic: width/height less than 1");
            System.Diagnostics.Debug.Assert(newwidth>1 & newheight>1, "Spline2DResampleBicubic: width/height less than 1");
            
            //
            // Prepare
            //
            mw = Math.Max(oldwidth, newwidth);
            mh = Math.Max(oldheight, newheight);
            b = new double[newheight-1+1, newwidth-1+1];
            buf = new double[oldheight-1+1, newwidth-1+1];
            x = new double[Math.Max(mw, mh)-1+1];
            y = new double[Math.Max(mw, mh)-1+1];
            
            //
            // Horizontal interpolation
            //
            for(i=0; i<=oldheight-1; i++)
            {
                
                //
                // Fill X, Y
                //
                for(j=0; j<=oldwidth-1; j++)
                {
                    x[j] = (double)(j)/((double)(oldwidth-1));
                    y[j] = a[i,j];
                }
                
                //
                // Interpolate and place result into temporary matrix
                //
                spline1d.spline1dbuildcubic(x, y, oldwidth, 0, 0.0, 0, 0.0, ref c);
                for(j=0; j<=newwidth-1; j++)
                {
                    buf[i,j] = spline1d.spline1dcalc(ref c, (double)(j)/((double)(newwidth-1)));
                }
            }
            
            //
            // Vertical interpolation
            //
            for(j=0; j<=newwidth-1; j++)
            {
                
                //
                // Fill X, Y
                //
                for(i=0; i<=oldheight-1; i++)
                {
                    x[i] = (double)(i)/((double)(oldheight-1));
                    y[i] = buf[i,j];
                }
                
                //
                // Interpolate and place result into B
                //
                spline1d.spline1dbuildcubic(x, y, oldheight, 0, 0.0, 0, 0.0, ref c);
                for(i=0; i<=newheight-1; i++)
                {
                    b[i,j] = spline1d.spline1dcalc(ref c, (double)(i)/((double)(newheight-1)));
                }
            }
        }


        /*************************************************************************
        Bilinear spline resampling

        Input parameters:
            A           -   function values at the old grid,
                            array[0..OldHeight-1, 0..OldWidth-1]
            OldHeight   -   old grid height, OldHeight>1
            OldWidth    -   old grid width, OldWidth>1
            NewHeight   -   new grid height, NewHeight>1
            NewWidth    -   new grid width, NewWidth>1

        Output parameters:
            B           -   function values at the new grid,
                            array[0..NewHeight-1, 0..NewWidth-1]

          -- ALGLIB routine --
             09.07.2007
             Copyright by Bochkanov Sergey
        *************************************************************************/
        public static void spline2dresamplebilinear(ref double[,] a,
            int oldheight,
            int oldwidth,
            ref double[,] b,
            int newheight,
            int newwidth)
        {
            int i = 0;
            int j = 0;
            int l = 0;
            int c = 0;
            double t = 0;
            double u = 0;

            b = new double[newheight-1+1, newwidth-1+1];
            for(i=0; i<=newheight-1; i++)
            {
                for(j=0; j<=newwidth-1; j++)
                {
                    l = i*(oldheight-1)/(newheight-1);
                    if( l==oldheight-1 )
                    {
                        l = oldheight-2;
                    }
                    u = (double)(i)/((double)(newheight-1))*(oldheight-1)-l;
                    c = j*(oldwidth-1)/(newwidth-1);
                    if( c==oldwidth-1 )
                    {
                        c = oldwidth-2;
                    }
                    t = (double)(j*(oldwidth-1))/((double)(newwidth-1))-c;
                    b[i,j] = (1-t)*(1-u)*a[l,c]+t*(1-u)*a[l,c+1]+t*u*a[l+1,c+1]+(1-t)*u*a[l+1,c];
                }
            }
        }


        /*************************************************************************
        Internal subroutine.
        Calculation of the first derivatives and the cross-derivative.
        *************************************************************************/
        private static void bicubiccalcderivatives(ref double[,] a,
            ref double[] x,
            ref double[] y,
            int m,
            int n,
            ref double[,] dx,
            ref double[,] dy,
            ref double[,] dxy)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            double[] xt = new double[0];
            double[] ft = new double[0];
            double[] c = new double[0];
            double s = 0;
            double ds = 0;
            double d2s = 0;
            double v = 0;

            dx = new double[m-1+1, n-1+1];
            dy = new double[m-1+1, n-1+1];
            dxy = new double[m-1+1, n-1+1];
            
            //
            // dF/dX
            //
            xt = new double[n-1+1];
            ft = new double[n-1+1];
            for(i=0; i<=m-1; i++)
            {
                for(j=0; j<=n-1; j++)
                {
                    xt[j] = x[j];
                    ft[j] = a[i,j];
                }
                spline3.buildcubicspline(xt, ft, n, 0, 0.0, 0, 0.0, ref c);
                for(j=0; j<=n-1; j++)
                {
                    spline3.splinedifferentiation(ref c, x[j], ref s, ref ds, ref d2s);
                    dx[i,j] = ds;
                }
            }
            
            //
            // dF/dY
            //
            xt = new double[m-1+1];
            ft = new double[m-1+1];
            for(j=0; j<=n-1; j++)
            {
                for(i=0; i<=m-1; i++)
                {
                    xt[i] = y[i];
                    ft[i] = a[i,j];
                }
                spline3.buildcubicspline(xt, ft, m, 0, 0.0, 0, 0.0, ref c);
                for(i=0; i<=m-1; i++)
                {
                    spline3.splinedifferentiation(ref c, y[i], ref s, ref ds, ref d2s);
                    dy[i,j] = ds;
                }
            }
            
            //
            // d2F/dXdY
            //
            xt = new double[n-1+1];
            ft = new double[n-1+1];
            for(i=0; i<=m-1; i++)
            {
                for(j=0; j<=n-1; j++)
                {
                    xt[j] = x[j];
                    ft[j] = dy[i,j];
                }
                spline3.buildcubicspline(xt, ft, n, 0, 0.0, 0, 0.0, ref c);
                for(j=0; j<=n-1; j++)
                {
                    spline3.splinedifferentiation(ref c, x[j], ref s, ref ds, ref d2s);
                    dxy[i,j] = ds;
                }
            }
        }
    }
}