/* interpolation/interp2d.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 <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>

/**
 * Triggers a GSL error if the argument is not equal to GSL_SUCCESS.
 * If the argument is GSL_SUCCESS, this does nothing.
 */
#define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s),  GSL_NAN); }

#define IDX2D(i, j, w) ((j) * ((w)->xsize) + (i))

gsl_interp2d *
gsl_interp2d_alloc(const gsl_interp2d_type * T, const size_t xsize,
                   const size_t ysize)
{
  gsl_interp2d * interp;

  if (xsize < T->min_size || ysize < T->min_size)
    {
      GSL_ERROR_NULL ("insufficient number of points for interpolation type",
                      GSL_EINVAL);
    }

  interp = (gsl_interp2d *) calloc(1, sizeof(gsl_interp2d));
  if (interp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for gsl_interp2d struct",
                      GSL_ENOMEM);
    }

  interp->type = T;
  interp->xsize = xsize;
  interp->ysize = ysize;

  if (interp->type->alloc == NULL)
    {
      interp->state = NULL;
      return interp;
    }

  interp->state = interp->type->alloc(xsize, ysize);
  if (interp->state == NULL)
    {
      free(interp);
      GSL_ERROR_NULL ("failed to allocate space for gsl_interp2d state",
                      GSL_ENOMEM);
    }

  return interp;
} /* gsl_interp2d_alloc() */

void
gsl_interp2d_free (gsl_interp2d * interp)
{
  RETURN_IF_NULL(interp);

  if (interp->type->free)
    interp->type->free(interp->state);

  free(interp);
} /* gsl_interp2d_free() */

int
gsl_interp2d_init (gsl_interp2d * interp, const double xarr[], const double yarr[],
                   const double zarr[], const size_t xsize, const size_t ysize)
{
  size_t i;

  if (xsize != interp->xsize || ysize != interp->ysize)
    {
      GSL_ERROR("data must match size of interpolation object", GSL_EINVAL);
    }

  for (i = 1; i < xsize; i++)
    {
      if (xarr[i-1] >= xarr[i])
        {
          GSL_ERROR("x values must be strictly increasing", GSL_EINVAL);
        }
    }

  for (i = 1; i < ysize; i++)
    {
      if (yarr[i-1] >= yarr[i])
        {
          GSL_ERROR("y values must be strictly increasing", GSL_EINVAL);
        }
    }

  interp->xmin = xarr[0];
  interp->xmax = xarr[xsize - 1];
  interp->ymin = yarr[0];
  interp->ymax = yarr[ysize - 1];

  {
    int status = interp->type->init(interp->state, xarr, yarr, zarr,
                                    xsize, ysize);
    return status;
  }
} /* gsl_interp2d_init() */

/*
 * A wrapper function that checks boundary conditions, calls an evaluator
 * which implements the actual calculation of the function value or 
 * derivative etc., and checks the return status.
 */
static int
interp2d_eval(int (*evaluator)(const void *, const double xa[], const double ya[],
                               const double za[], size_t xsize, size_t ysize,
                               double x, double y, gsl_interp_accel *,
                               gsl_interp_accel *, double * z),
              const gsl_interp2d * interp, const double xarr[],
              const double yarr[], const double zarr[],
              const double x, const double y,
              gsl_interp_accel * xa, gsl_interp_accel * ya,
              double * result)
{
  if (x < interp->xmin || x > interp->xmax)
    {
      GSL_ERROR ("interpolation x value out of range", GSL_EDOM);
    }
  else if (y < interp->ymin || y > interp->ymax)
    {
      GSL_ERROR ("interpolation y value out of range", GSL_EDOM);
    }

  return evaluator(interp->state, xarr, yarr, zarr,
                   interp->xsize, interp->ysize,
                   x, y, xa, ya, result);
}

/*
 * Another wrapper function that serves as a drop-in replacement for
 * interp2d_eval but does not check the bounds. This can be used
 * for extrapolation.
 */
static int
interp2d_eval_extrap(int (*evaluator)(const void *, const double xa[],
                                      const double ya[], const double za[],
                                      size_t xsize, size_t ysize,
                                      double x, double y,
                                      gsl_interp_accel *,
                                      gsl_interp_accel *, double * z),
                     const gsl_interp2d * interp, const double xarr[],
                     const double yarr[], const double zarr[],
                     const double x, const double y,
                     gsl_interp_accel * xa, gsl_interp_accel * ya,
                     double * result)
{
  return evaluator(interp->state, xarr, yarr, zarr,
                   interp->xsize, interp->ysize, x, y, xa, ya, result);
}

double
gsl_interp2d_eval (const gsl_interp2d * interp, const double xarr[],
                   const double yarr[], const double zarr[],
                   const double x, const double y,
                   gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
} /* gsl_interp2d_eval() */

double
gsl_interp2d_eval_extrap (const gsl_interp2d * interp,
                          const double xarr[],
                          const double yarr[],
                          const double zarr[],
                          const double x,
                          const double y,
                          gsl_interp_accel * xa,
                          gsl_interp_accel * ya)
{
  double z;
  int status =
    interp2d_eval_extrap(interp->type->eval, interp,
                         xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_e (const gsl_interp2d * interp, const double xarr[],
                     const double yarr[], const double zarr[],
                     const double x, const double y,
                     gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
} /* gsl_interp2d_eval_e() */

int
gsl_interp2d_eval_e_extrap (const gsl_interp2d * interp,
                            const double xarr[], const double yarr[],
                            const double zarr[], const double x,
                            const double y, gsl_interp_accel * xa,
                            gsl_interp_accel * ya, double * z)
{
  return interp2d_eval_extrap(interp->type->eval, interp,
                              xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_x (const gsl_interp2d * interp, const double xarr[],
                           const double yarr[], const double zarr[],
                           const double x, const double y,
                           gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_x_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_x_e (const gsl_interp2d * interp, const double xarr[],
                             const double yarr[], const double zarr[],
                             const double x, const double y,
                             gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_x, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_y (const gsl_interp2d * interp, const double xarr[],
                           const double yarr[], const double zarr[],
                           const double x, const double y,
                           gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_y_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_y_e (const gsl_interp2d * interp, const double xarr[],
                             const double yarr[], const double zarr[],
                             const double x, const double y,
                             gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_y, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_xx (const gsl_interp2d * interp, const double xarr[],
                            const double yarr[], const double zarr[],
                            const double x, const double y,
                            gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_xx_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_xx_e (const gsl_interp2d * interp, const double xarr[],
                              const double yarr[], const double zarr[],
                              const double x, const double y,
                              gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_xx, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_yy (const gsl_interp2d * interp, const double xarr[],
                            const double yarr[], const double zarr[],
                            const double x, const double y,
                            gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_yy_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_yy_e (const gsl_interp2d * interp, const double xarr[],
                              const double yarr[], const double zarr[],
                              const double x, const double y,
                              gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_yy, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

double
gsl_interp2d_eval_deriv_xy (const gsl_interp2d * interp, const double xarr[],
                            const double yarr[], const double zarr[],
                            const double x, const double y,
                            gsl_interp_accel * xa, gsl_interp_accel * ya)
{
  double z;
  int status = gsl_interp2d_eval_deriv_xy_e(interp, xarr, yarr, zarr, x, y, xa, ya, &z);
  DISCARD_STATUS(status)
  return z;
}

int
gsl_interp2d_eval_deriv_xy_e (const gsl_interp2d * interp, const double xarr[],
                              const double yarr[], const double zarr[],
                              const double x, const double y,
                              gsl_interp_accel * xa, gsl_interp_accel * ya, double * z)
{
  return interp2d_eval(interp->type->eval_deriv_xy, interp,
                       xarr, yarr, zarr, x, y, xa, ya, z);
}

size_t
gsl_interp2d_type_min_size(const gsl_interp2d_type * T)
{
  return T->min_size;
}

size_t
gsl_interp2d_min_size(const gsl_interp2d * interp)
{
  return interp->type->min_size;
}

const char *
gsl_interp2d_name(const gsl_interp2d * interp)
{
  return interp->type->name;
}

size_t
gsl_interp2d_idx(const gsl_interp2d * interp,
                 const size_t i, const size_t j)
{
  if (i >= interp->xsize)
    {
      GSL_ERROR_VAL ("x index out of range", GSL_ERANGE, 0);
    }
  else if (j >= interp->ysize)
    {
      GSL_ERROR_VAL ("y index out of range", GSL_ERANGE, 0);
    }
  else
    {
      return IDX2D(i, j, interp);
    }
} /* gsl_interp2d_idx() */

int
gsl_interp2d_set(const gsl_interp2d * interp, double zarr[],
                 const size_t i, const size_t j, const double z)
{
  if (i >= interp->xsize)
    {
      GSL_ERROR ("x index out of range", GSL_ERANGE);
    }
  else if (j >= interp->ysize)
    {
      GSL_ERROR ("y index out of range", GSL_ERANGE);
    }
  else
    {
      zarr[IDX2D(i, j, interp)] = z;
      return GSL_SUCCESS;
    }
} /* gsl_interp2d_set() */

double
gsl_interp2d_get(const gsl_interp2d * interp, const double zarr[],
                 const size_t i, const size_t j)
{
  if (i >= interp->xsize)
    {
      GSL_ERROR_VAL ("x index out of range", GSL_ERANGE, 0);
    }
  else if (j >= interp->ysize)
    {
      GSL_ERROR_VAL ("y index out of range", GSL_ERANGE, 0);
    }
  else
    {
      return zarr[IDX2D(i, j, interp)];
    }
} /* gsl_interp2d_get() */

#undef IDX2D