/* 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