package algorithmrepository;

import jafama.FastMath;
import algorithmrepository.exceptions.NotImplementedException;


/**
 * A class for doing 2D cubic interpolation
 * on uniform grids. The cubic interpolator give continous
 * function values and first derivatives across grid points.
 *
 * Both function values, first derivatives and second order mixed derivatives
 * can be interpolated.
 *
 * The interpolator can be initialised either with function values and
 * derivatives at each point, or just function values, in which case the
 * neccessary derivatives at the grid points are calculated from finite
 * differences on the grid.
 *
 * The static methods of the class can be used directly, but the
 * standard way of using it is through:
 *
 * 	CubicInterpolation2D ip = new CubicInterpolation2D(double[] x, double[] y, double[][] f);
 *  ip.eval(3,4);
 *  ip.eval(new double[] {1,2,3 }, new double[] {2,4,8 });
 *  etc.
 *
 *  // Derivatives (dfdx,dfdy,d2fdxdy) can be interpolated through:
 *  double[] der = nwe double[3];
 *  ip.eval(3,4,der);
 *  System.out.println("dfdx(3,4)="+der[0]);
 *  System.out.println("dfdy(3,4)="+der[1]);
 *  System.out.println("d2fdxdy(3,4)="+der[2]);
 *
 *  See http://www.algorithmrepository.org
 *
 */
public class CubicInterpolation2D implements Interpolation2D {
    private double[] xp, yp;
    private double[][] fp, dfdx, dfdy, d2fdxdy;
    private int numx, numy;
    private double minx, maxx;
    private double miny, maxy;
    private double dx, dy;
    private double[][][][] coeffs;
    private double[] f_cell = new double[4], dfdx_cell=new double[4], dfdy_cell=new double[4], d2fdxdy_cell=new double[4];
    private double extrapValue = Double.NaN;
  
    
    /**
     * Initialises bicubic interpolation for a uniform grid.
     *
     * Neccessary derivatives are calculated using finite differences
     * from the given function values on the grid.
     * 
     * Sets up to return NaN outside grid
     *
     * @param xp The x-coordinates of the grid points
     * @param yp The y-coordinates of the grid points.
     * @param fp The function values at the grid points
     * fp[i][j] is the function value at the ith x-coordinate and jth y-coordinate.
     */
    public CubicInterpolation2D(double[] xp, double[] yp, double[][] fp) {
        this(xp,yp,fp,null,null,null,Double.NaN);
    }

    public CubicInterpolation2D(double[] xp, double[] yp, double[][] fp, double[][] dfdx, double[][] dfdy, double[][] d2fdxdy) {
        throw new NotImplementedException(); // <oliford 16/5/11> What was supposed to happen here? Does anyone know?
    }
    
    /**
     * Initialises bicubic interpolation for a uniform grid when
     * derivatives are of the funcion are known.
     *
     * @param xp The x-coordinates of the grid points
     * @param yp The y-coordinates of the grid points.
     * @param fp The function values at the grid points
     * fp[i][j] is the function value at the ith x-coordinate and jth y-coordinate.
     * @param dfdx dfdx at the grid points.
     * @param dfdy  dfdy at the grid points.
     * @param d2fdxdy  d2fdxdy at the grid points.
     * @param extrapValue Value to return if outside grid (extrapolation)
     */
    public CubicInterpolation2D(double[] xp, double[] yp, double[][] fp, double[][] dfdx, double[][] dfdy, double[][] d2fdxdy, double extrapValue) {
        this.xp = xp;
        this.yp = yp;
        this.fp = fp;
        this.dfdx = dfdx;
        this.dfdy = dfdy;
        this.d2fdxdy = d2fdxdy;
        this.numx = xp.length;
        this.numy = yp.length;
        this.minx = xp[0];
        this.maxx = xp[numx-1];
        this.dx = xp[1]-xp[0];
        this.dy = yp[1]-yp[0];
        this.miny = yp[0];
        this.maxy = yp[numy-1];
        this.coeffs = new double[numx-1][numy-1][][];
        this.extrapValue = extrapValue;
        
        if (fp.length != xp.length || fp[0].length != yp.length) {
        	throw new RuntimeException("Matrix of function values to CubicInterpolation wrong dimension, has ["+fp.length+"]["+fp[0].length+"], but should be ["+xp.length+"]["+yp.length+"]");
        }
    }


    /**
     * Single bicubic interpolation. The bicubic interpolation
     * uses for each point in a 2x2 grid: the function value,
     * the partial derivative of the function in the x-direction,
     * the partial derivative in the y-direction, the mixed second order
     * derivative d2f/dxdy.
     *
     * @param x1 left x-coordinate
     * @param x2 right x-coordinate
     * @param y1 lower y-coordinate
     * @param y2 upper y-coordinate
     * @param fp The four function values corresponding to the given
     * grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param dfdx The four first derivatives in the x-direction corresponding
     * to the given grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param dfdy The four first derivatives in the y-direction corresponding
     * to the given grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param d2fdxdy The four mixed derivatives corresponding
     * to the given grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param x The x-coordinate of the point that should be interpolated.
     * @param y The y-coordinate of the point that should be interpolated.
     * @param der An array of derivatives of length 3 that has to be allocated before the
     * call if used:
     * der[0] = dfdx
     * der[1] = dfdy
     * der[2] = d2fdxdy
     *
     * @return The interpolated value of the function at the given point.
     */
    public final static double eval(double x1, double x2, double y1, double y2, double[] fp, double[] dfdx, double[] dfdy, double[] d2fdxdy, double x, double y, double[] der) {
        double dx = (x2-x1);
        double dy = (y2-y1);
        double u = (x-x1)/dx;
        double v = (y-y1)/dy;

        double[][] c = cubicInterpolationCoeffs2D(dx, dy, fp, dfdx, dfdy, d2fdxdy);
        double ret = interp(u,v,dx,dy,c,der);

        return ret;
    }


    /**
     *
     * Does a bicubic interpolation on uniform grid given in the constructor.
     *
     * Extrapolation: Doesn't do extrapolation, returns 'extrapValue' if outside grid.
     *
     * @param x The x-coordinates where the function should be interpolated
     * @param y The y-coordinates where the function should be interpolated
     *
     * @return The interpolated values of the function.
     */
    public final double[] eval(double[] x, double[] y) {
        return eval(x,y,null);
    }


    /**
     *
     * Does a bicubic interpolation on uniform grid given in the constructor.
     *
     * Extrapolation: Doesn't do extrapolation, returns 'extrapValue' if outside grid.
     *
     * @param x The x-coordinate where the function should be interpolated
     * @param y The y-coordinate where the function should be interpolated
     *
     * @return The interpolated value of the function.
     */
    public final double eval(double x, double y) {
        return eval(x,y,null);
    }
    
    /**
     * Evaluates the function on a 2D grid with given knots in the x and y-dimension.
     * 
     * @param xKnots The knots along the x dimension.
     * @param yKnots The knots along the y dimension.
     * @return ret[i][j] is the function value for know xi, yj.
     */
    public double[][] evalGrid(double[] xKnots, double[] yKnots) {
        double[][] ret = new double[xKnots.length][yKnots.length];
        
        for(int i=0;i<ret.length;++i) {
            for(int j=0;j<ret[0].length;++j) {
                ret[i][j] = eval(xKnots[i], yKnots[j]);
            }
        }
        
        return ret;
    }
    

    /**
     *
     * Does a bicubic interpolation on uniform grid given in the constructor.
     *
     * Extrapolation: Doesn't do extrapolation, returns 'extrapValue' if outside grid.
     *
     * @param x The x-coordinate where the function should be interpolated
     * @param y The y-coordinate where the function should be interpolated
     * @param der On return der contains dfdx, dfdy, d2fdxdy for the interpolation point.
     * If this functionality is sought, der has to be allocated by the caller to der[3].
     *
     * @return The interpolated value of the function.
     */
    public final double eval(double x, double y, double[] der) {
        double[][] c;
        
        //int ix = (int) Math.floor((numx-1)*(x-minx)/(maxx-minx));
        //int iy = (int) Math.floor((numy-1)*(y-miny)/(maxy-miny));

        
        //maybe it was faster to do this after the ix= and iy=, in
        // order to do the tests in integer arithmetic
        // but then we have to use Math.floor() rather than just (int)
        if(x < minx || y < miny || x > maxx || y > maxy)
            return extrapValue;
        
        int ix = (int)((numx-1)*(x-minx)/(maxx-minx));
        int iy = (int)((numy-1)*(y-miny)/(maxy-miny));

        //we should be able to hit the end exactly, but no futher
        if (x == maxx) ix -= 1;
        if (y == maxy) iy -= 1;
        
        if (coeffs[ix][iy] == null) { // Recalculate parameters for cell
            if (dfdx == null) {
                gridCellDerivatives(xp, yp, fp, ix, iy, f_cell, dfdx_cell, dfdy_cell, d2fdxdy_cell);
            }
            else {
                f_cell[0] = fp[ix][iy];
                f_cell[1] = fp[ix+1][iy];
                f_cell[2] = fp[ix+1][iy+1];
                f_cell[3] = fp[ix][iy+1];

                dfdx_cell[0] = dfdx[ix][iy];
                dfdx_cell[1] = dfdx[ix+1][iy];
                dfdx_cell[2] = dfdx[ix+1][iy+1];
                dfdx_cell[3] = dfdx[ix][iy+1];

                dfdy_cell[0] = dfdy[ix][iy];
                dfdy_cell[1] = dfdy[ix+1][iy];
                dfdy_cell[2] = dfdy[ix+1][iy+1];
                dfdy_cell[3] = dfdy[ix][iy+1];

                d2fdxdy_cell[0] = d2fdxdy[ix][iy];
                d2fdxdy_cell[1] = d2fdxdy[ix+1][iy];
                d2fdxdy_cell[2] = d2fdxdy[ix+1][iy+1];
                d2fdxdy_cell[3] = d2fdxdy[ix][iy+1];
            }

            c = cubicInterpolationCoeffs2D(dx,dy,f_cell,dfdx_cell,dfdy_cell,d2fdxdy_cell);
            coeffs[ix][iy] = c;
        }
        else {
            c = coeffs[ix][iy];
        }

        double u = (x-xp[ix])/dx;
        double v = (y-yp[iy])/dy;
        double ret = interp(u,v,dx,dy,c,der);

        return ret;
    }




    /**
     *
     * Does a bicubic interpolation on uniform grid given in the constructor.
     *
     * Extrapolation: Doesn't do extrapolation, returns 'extrapValue' if outside grid.
     *
     * @param x The x-coordinates where the function should be interpolated
     * @param y The y-coordinates where the function should be interpolated
     * @param der On return der contains dfdx, dfdy and optionally d2fdxdy for each interpolation point.
     * If this functionality is sought, der has to be allocated by the caller to der[numPoints][2 or 3],
     * where numPoints is the length of the x-vector above.
     *
     * @return The interpolated values of the function.
     */
    public final double[] eval(double[] x, double[] y, double[][] der) {
        double[] ret = new double[x.length];

        if (der == null) {
            for(int i=0;i<x.length;++i) {
                ret[i] = eval(x[i],y[i],null);
            }
        }
        else {
            for(int i=0;i<x.length;++i) {
                ret[i] = eval(x[i],y[i],der[i]);
            }
        }

        return ret;
    }


    /**
     * Interpolates a single value withing a cell given precalculated
     * coefficients for the polynomial (from (@link cubicInterpolationCoeffs2D()}.
     *
     * @param u Normalised [0,1] x-coordinate within the cell
     * @param v Normalised [0,1] y-coordinate within the cell
     * @param dx Size of cell in x-direction
     * @param dy Size of cell in y-direction
     * @param c 4x4 coefficient matrix (@see cubicInterpolationCoeffs2D()}.
     * @param der If allocated (double[3]) returns the derivatives:
     * der[0] = dfdx
     * der[1] = dfdy
     * if(der.length > 2) der[2] = d2fdxdy
     * @return The interpolated value.
     */
    protected final static double interp(double u, double v, double dx, double dy, double[][] c, double[] der) {
        double val = c[0][0]             + c[0][1]  *v           + c[0][2]  *v*v         + c[0][3]  *v*v*v
                    + c[1][0] *u        + c[1][1] *u *v         + c[1][2] *u *v*v       + c[1][3] *u *v*v*v
                    + c[2][0] *u*u      + c[2][1] *u*u *v       + c[2][2] *u*u *v*v     + c[2][3] *u*u *v*v*v
                    + c[3][0] *u*u*u    + c[3][1] *u*u*u *v     + c[3][2] *u*u*u *v*v   + c[3][3] *u*u*u *v*v*v;
        
        if (der != null) {
            der[0] =  (1*c[1][0]   + 1*c[1][1]  *v + 1*c[1][2]  *v*v + 1*c[1][3]  *v*v*v
                        + 2*c[2][0] *u  + 2*c[2][1] *u *v + 2*c[2][2] *u *v*v + 2*c[2][3] *u *v*v*v
                        + 3*c[3][0] *u*u  + 3*c[3][1] *u*u *v + 3*c[3][2] *u*u *v*v + 3*c[3][3] *u*u *v*v*v) / dx;

        
            
            der[1] =  (1*c[0][1]   + 2*c[0][2]  *v + 3*c[0][3]  *v*v
                        + 1*c[1][1] *u  + 2*c[1][2] *u *v + 3*c[1][3] *u *v*v
                        + 1*c[2][1] *u*u  + 2*c[2][2] *u*u *v + 3*c[2][3] *u*u *v*v
                        + 1*c[3][1] *u*u*u  + 2*c[3][2] *u*u*u *v + 3*c[3][3] *u*u*u *v*v) / dy;
            
            if(der.length > 2){
                der[2] = (1*1*c[1][1]   + 1*2*c[1][2]  *v + 1*3*c[1][3]  *v*v
                            + 2*1*c[2][1] *u  + 2*2*c[2][2] *u *v + 2*3*c[2][3] *u *v*v
                            + 3*1*c[3][1] *u*u  + 3*2*c[3][2] *u*u *v + 3*3*c[3][3] *u*u *v*v) / (dx*dy);
            }
                           
        }

        return val;
    }



    /**
     * Precalculates coefficients neccessary for cubic interpolation.
     * This is a convenience method to get the inverted
     * coefficient matrix for the linear equation system in the 16 free
     * parameters of a 3rd degree polynomial in 2 dimensions, from
     * values of the function and its derivatives:
     *
     *    a*c=g
     *
     *    where g = (f00,f01,f11,f01,fu00,fu01,fu11,fu01,fv00,fv01,fv11,fv01,fuv00,fuv01,fuv11,fuv01)'
     *    fij is function value at the i,j corner point of the cell
     *    fuij is dfdu at the i,j corner point of the cell (u is normalised x within the cell)
     *    fvij is dfdv at the i,j corner point of the cell (v is normalised y within the cell)
     *    fuvij is d2fdudv at the i,j corner point of the cell
     *
     *    The coefficients c are then
     *
     *    c = ainv*g
     *
     *    where ainv is the inverse of a.
     *
     * In the return value from this method c is arranged as a 4x4 matrix so that
     * the function values and its derivatives can be calculated as:
     *
     * (u,v are normalised ([0,1]) x,y-coordinates within the cell)
     *
     * f(u,v) = sum_i( sum_j( c_ij*u^(i)*v^(j) ) )
     * dfdu(u,v) = sum_i( sum_j( i*c_ij*u^(i-1)*v^(j) ) )
     * dfdv(u,v) = sum_i( sum_j( j*c_ij*u^(i)*v^(j-1) ) )
     * d2fdudv(u,v) = sum_i( sum_j( i*j*c_ij*u^(i-1)*v^(j-1) ) )
     *
     * Remember: to go from derivatives in (u,v) to derivatives in (x,y) do the
     * following transformation:
     *
     * dfdx = dx*dfdu
     * dfdy = dy*dfdv
     * d2fdxdy = dx*dy*d2fdudv
     *
     * @param dx Extent of grid cell in x-direction
     * @param dy Extent of grid cell in y-direction
     * @param fp The four function values corresponding to the given
     * grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param dfdx The four first derivatives in the x-direction corresponding
     * to the given grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param dfdy The four first derivatives in the y-direction corresponding
     * to the given grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param d2fdxdy The four mixed derivatives corresponding
     * to the given grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @return A 4x4 matrix of values from which interpolated
     * values and interpolated derivatives can be calculated within the given cell by
     * the formulas above.
     *
     */
    public final static double[][] cubicInterpolationCoeffs2D(double dx, double dy, double[] fp, double[] dfdx, double[] dfdy, double[] d2fdxdy) {
        double[][] c = new double[4][4];
        int p,l;
        double dxdy = dx*dy;

        for(int i=0;i<4;++i) {
            for(int j=0;j<4;++j) {
                c[i][j] = 0;
                p = i*4+j;
                l = 0;
                for(int k=0;k<4;++k,++l) {
                    c[i][j] += ainv[p][l]*fp[k];
                }
                for(int k=0;k<4;++k,++l) {
                    c[i][j] += ainv[p][l]*dfdx[k]*dx;
                }
                for(int k=0;k<4;++k,++l) {
                    c[i][j] += ainv[p][l]*dfdy[k]*dy;
                }
                for(int k=0;k<4;++k,++l) {
                    c[i][j] += ainv[p][l]*d2fdxdy[k]*dxdy;
                }
            }
        }
        return c;
    }

    private final static int[][] ainv = {
        { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
        { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
        { -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0 },
        { 2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0 },
        { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
        { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
        { 0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1 },
        { 0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1 },
        { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
        { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
        { 9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2 },
        { -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2 },
        { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
        { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
        { -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1 },
        { 4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1 } };



    /**
     * Calculates dfdx, dfdy, d2fdxdy on each corner point of a grid cell in
     * a uniform grid.
     *
     * Derivatives for internal points use central differences.
     * Derivatives for boundary points use forward and backward differences.
     *
     * @param xp The x-coordinates of the grid points
     * @param yp The y-coordinates of the grid points.
     * @param fp The function values at the grid points
     * fp[i][j] is the function value at the ith x-coordinate and jth y-coordinate.
     * @param ix Lower x-index of grid cell
     * @param iy Lower y-index of grid cell
     * @param fv On return contains the function values for the corner points with (0,0) at (ix,iy)
     * @param dfdx On return contains dfdx for the corner points with (0,0) at (ix,iy)
     * @param dfdy On return contains dfdy for the corner points with (0,0) at (ix,iy)
     * @param d2fdxdy On return contains d2fdxdy for the corner points with (0,0) at (ix,iy)
     */
    public final static void gridCellDerivatives(double[] xp, double[] yp, double[][] fp, int ix, int iy, double[] fv, double[] dfdx, double[] dfdy, double[] d2fdxdy) {
        double[] der = new double[3];

        gridPointDerivatives(xp,yp,fp,ix,iy,der);
        fv[0] = fp[ix][iy];
        dfdx[0] = der[0];
        dfdy[0] = der[1];
        d2fdxdy[0] = der[2];

        gridPointDerivatives(xp,yp,fp,ix+1,iy,der);
        fv[1] = fp[ix+1][iy];
        dfdx[1] = der[0];
        dfdy[1] = der[1];
        d2fdxdy[1] = der[2];

        gridPointDerivatives(xp,yp,fp,ix+1,iy+1,der);
        fv[2] = fp[ix+1][iy+1];
        dfdx[2] = der[0];
        dfdy[2] = der[1];
        d2fdxdy[2] = der[2];

        gridPointDerivatives(xp,yp,fp,ix,iy+1,der);
        fv[3] = fp[ix][iy+1];
        dfdx[3] = der[0];
        dfdy[3] = der[1];
        d2fdxdy[3] = der[2];
    }


    /**
     * Calculates dfdx, dfdy, d2fdxdy for the specified grid point on
     * a uniform grid.
     *
     * Derivatives for internal points use central differences.
     * Derivatives for boundary points use forward and backward differences.
     *
     * @param xp The x-coordinates of the grid points
     * @param yp The y-coordinates of the grid points.
     * @param fp The function values at the grid points
     * fp[i][j] is the function value at the ith x-coordinate and jth y-coordinate.
     * @param ix x-index of grid point
     * @param iy y-index of grid point
     * @param der On return (must be allocated by caller) contains:
     * der[0] = dfdx
     * der[1] = dfdy
     * der[2] = d2fdxdy
     */
    public final static void gridPointDerivatives(double[] xp, double[] yp, double[][] fp, int ix, int iy, double[] der) {
        int numx = xp.length;
        int numy = yp.length;
        double dx = xp[1]-xp[0];
        double dy = yp[1]-yp[0];
        double dfdx,dfdy,d2fdxdy;

        // Do edge derivatives using forward difference
        // and inner derivatives using central difference

        // dfdx
        if (ix == 0) {
            dfdx = (fp[ix+1][iy] - fp[ix][iy])/dx;
        }
        else if (ix == numx-1) {
            dfdx = (fp[ix][iy] - fp[ix-1][iy])/dx;

        } else {
            dfdx = (fp[ix+1][iy] - fp[ix-1][iy])/(2*dx);
        }

        // dfdy
        if (iy == 0) {
            dfdy = (fp[ix][iy+1] - fp[ix][iy])/dy;
        }
        else if (iy == numy-1) {
            dfdy = (fp[ix][iy] - fp[ix][iy-1])/dy;

        } else {
            dfdy = (fp[ix][iy+1] - fp[ix][iy-1])/(2*dy);
        }

        // d2fdxdy
        if (ix == 0 && iy != numy-1) {
            d2fdxdy = (fp[ix+1][iy+1] - fp[ix][iy+1] - fp[ix+1][iy] + fp[ix][iy])/(dx*dy);
        }
        else if (ix == 0 && iy == numy -1) {
            d2fdxdy = (fp[ix+1][iy] - fp[ix][iy] - fp[ix+1][iy-1] + fp[ix][iy-1])/(dx*dy);
        }
        else if (ix == numx-1 && iy != numy-1) {
            d2fdxdy = (fp[ix][iy+1] - fp[ix-1][iy+1] - fp[ix][iy] + fp[ix-1][iy])/(dx*dy);
        }
        else if (ix == numx-1 && iy == numy-1) {
            d2fdxdy = (fp[ix][iy] - fp[ix-1][iy] - fp[ix][iy-1] + fp[ix-1][iy-1])/(dx*dy);
        }
        else if (iy == 0) {
            d2fdxdy = (fp[ix+1][iy+1] - fp[ix][iy+1] - fp[ix+1][iy] + fp[ix][iy])/(dx*dy);
        }
        else if (iy == numy-1) {
            d2fdxdy = (fp[ix+1][iy] - fp[ix][iy] - fp[ix+1][iy-1] + fp[ix][iy-1])/(dx*dy);
        }
        else {
            d2fdxdy = (fp[ix+1][iy+1] - fp[ix-1][iy+1] - fp[ix+1][iy-1] + fp[ix-1][iy-1])/(4*dx*dy);
        }

        der[0] = dfdx;
        der[1] = dfdy;
        der[2] = d2fdxdy;
    }
    

    /**
     * Returns the vector of x-coordinates from which the grid is 
     * defined.
     * 
     * @return The vector of x-coordinates defining the grid division
     * in the x-direction.
     */
    public final double[] getxp() {
        return xp;
    }

    
    /**
     * Returns the vector of y-coordinates from which the grid is 
     * defined.
     * 
     * @return The vector of y-coordinates defining the grid division
     * in the y-direction.
     */    
    public final double[] getyp() {
        return yp;
    }
    
    
    /**
     * Returns the grid values. 
     * 
     * @return The matrix of grid values given in the constructor.
     * ret[i][j] is the function value corresponding to the grid cell with
     * lower left coordinates x[i],y[j].
     */    
    public final double[][] getfp() {
        return fp;
    }
    
    public void setfp(double fp[][]) {
        this.fp = fp;
        coeffs = new double[numx-1][numy-1][][];
    }
}