/* interpolation/bilinear.c

 * 

 * Copyright 2012 David Zaslavsky

 * 

 * 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; either version 3 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.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

#include <config.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_interp.h>

#include <gsl/gsl_interp2d.h>

#define IDX2D(i, j, xsize, ysize) ((j) * (xsize) + (i))

static int

bilinear_init(void * state, const double xa[], const double ya[],

              const double za[], size_t xsize, size_t ysize)

{

  return GSL_SUCCESS;

}

static int

bilinear_eval(const void * state, const double xarr[], const double yarr[],

              const double zarr[], size_t xsize, size_t ysize,

              double x, double y, gsl_interp_accel * xa,

              gsl_interp_accel * ya, double * z)

{

  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;

  double dx, dy;

  double t, u;

  size_t xi, yi;

  if (xa != NULL)

    xi = gsl_interp_accel_find(xa, xarr, xsize, x);

  else

    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)

    yi = gsl_interp_accel_find(ya, yarr, ysize, y);

  else

    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];

  xmax = xarr[xi + 1];

  ymin = yarr[yi];

  ymax = yarr[yi + 1];

  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];

  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];

  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];

  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];

  dx = xmax - xmin;

  dy = ymax - ymin;

  t = (x - xmin)/dx;

  u = (y - ymin)/dy;

  *z = (1.-t)*(1.-u)*zminmin + t*(1.-u)*zmaxmin + (1.-t)*u*zminmax + t*u*zmaxmax;

  return GSL_SUCCESS;

}

static int

bilinear_deriv_x(const void * state, const double xarr[],

                 const double yarr[], const double zarr[],

                 size_t xsize, size_t ysize, double x, double y,

                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)

{

  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;

  double dx, dy;

  double dt, u;

  size_t xi, yi;

  if (xa != NULL)

    xi = gsl_interp_accel_find(xa, xarr, xsize, x);

  else

    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)

    yi = gsl_interp_accel_find(ya, yarr, ysize, y);

  else

    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];

  xmax = xarr[xi + 1];

  ymin = yarr[yi];

  ymax = yarr[yi + 1];

  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];

  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];

  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];

  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];

  dx = xmax - xmin;

  dy = ymax - ymin;

  dt = 1./dx; /* partial t / partial x */

  u = (y - ymin)/dy;

  *z_p = dt*(-(1.-u)*zminmin + (1.-u)*zmaxmin - u*zminmax + u*zmaxmax);

  return GSL_SUCCESS;

}

static int

bilinear_deriv_y(const void * state, const double xarr[],

                 const double yarr[], const double zarr[],

                 size_t xsize, size_t ysize, double x, double y,

                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)

{

  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;

  double dx, dy;

  double t, du;

  size_t xi, yi;

  if (xa != NULL)

    xi = gsl_interp_accel_find(xa, xarr, xsize, x);

  else

    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)

    yi = gsl_interp_accel_find(ya, yarr, ysize, y);

  else

    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];

  xmax = xarr[xi + 1];

  ymin = yarr[yi];

  ymax = yarr[yi + 1];

  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];

  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];

  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];

  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];

  dx = xmax - xmin;

  dy = ymax - ymin;

  t = (x - xmin)/dx;

  du = 1./dy; /* partial u / partial y */

  *z_p = du*(-(1.-t)*zminmin - t*zmaxmin + (1.-t)*zminmax + t*zmaxmax);

  return GSL_SUCCESS;

}

static int

bilinear_deriv2(const void * state, const double xarr[],

                const double yarr[], const double zarr[],

                size_t xsize, size_t ysize, double x, double y,

                gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)

{

  *z_pp = 0.0;

  return GSL_SUCCESS;

}

static int

bilinear_derivxy(const void * state, const double xarr[],

                 const double yarr[], const double zarr[],

                 size_t xsize, size_t ysize, double x, double y,

                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)

{

  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;

  double dx, dy;

  double dt, du;

  size_t xi, yi;

  if (xa != NULL)

    xi = gsl_interp_accel_find(xa, xarr, xsize, x);

  else

    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)

    yi = gsl_interp_accel_find(ya, yarr, ysize, y);

  else

    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];

  xmax = xarr[xi + 1];

  ymin = yarr[yi];

  ymax = yarr[yi + 1];

  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];

  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];

  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];

  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];

  dx = xmax - xmin;

  dy = ymax - ymin;

  dt = 1./dx; /* partial t / partial x */

  du = 1./dy; /* partial u / partial y */

  *z_pp = dt*du*(zminmin-zmaxmin-zminmax+zmaxmax);

  return GSL_SUCCESS;

}

static const gsl_interp2d_type bilinear_type = {

  "bilinear",

  2,

  NULL,

  &bilinear_init,

  &bilinear_eval,

  &bilinear_deriv_x,

  &bilinear_deriv_y,

  &bilinear_deriv2,

  &bilinear_derivxy,

  &bilinear_deriv2,

  NULL

};

const gsl_interp2d_type * gsl_interp2d_bilinear = &bilinear_type;

#undef IDX2D