@cindex interpolation
@cindex spline

This chapter describes functions for performing interpolation.  The
library provides a variety of interpolation methods, including Cubic,
Akima, and Steffen splines.  The interpolation types are interchangeable,
allowing different methods to be used without recompiling.
Interpolations can be defined for both normal and periodic boundary
conditions.  Additional functions are available for computing
derivatives and integrals of interpolating functions. Routines
are provided for interpolating both one and two dimensional datasets.

These interpolation methods produce curves that pass through each
datapoint.  To interpolate noisy data with a smoothing curve see
@ref{Basis Splines}.

The functions described in this section are declared in the header files
@file{gsl_interp.h} and @file{gsl_spline.h}.

@menu
* 1D Introduction to Interpolation::  
* 1D Interpolation Functions::     
* 1D Interpolation Types::         
* 1D Index Look-up and Acceleration::  
* 1D Evaluation of Interpolating Functions::  
* 1D Higher-level Interface::      
* 1D Interpolation Example programs::  
* 1D Interpolation References and Further Reading::  
* 2D Introduction to Interpolation::  
* 2D Interpolation Functions::     
* 2D Interpolation Grids::     
* 2D Interpolation Types::         
* 2D Evaluation of Interpolating Functions::  
* 2D Higher-level Interface::      
* 2D Interpolation Example programs::  
@end menu

@node 1D Introduction to Interpolation
@section Introduction to 1D Interpolation

Given a set of data points @math{(x_1, y_1) \dots (x_n, y_n)} the
routines described in this section compute a continuous interpolating
function @math{y(x)} such that @math{y(x_i) = y_i}.  The interpolation
is piecewise smooth, and its behavior at the end-points is determined by
the type of interpolation used.

@node 1D Interpolation Functions
@section 1D Interpolation Functions

The interpolation function for a given dataset is stored in a
@code{gsl_interp} object.  These are created by the following functions.

@deftypefun {gsl_interp *} gsl_interp_alloc (const gsl_interp_type * @var{T}, size_t @var{size})
@tindex gsl_interp
This function returns a pointer to a newly allocated interpolation
object of type @var{T} for @var{size} data-points.
@end deftypefun

@deftypefun int gsl_interp_init (gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], size_t @var{size})
This function initializes the interpolation object @var{interp} for the
data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size
@var{size}.  The interpolation object (@code{gsl_interp}) does not save
the data arrays @var{xa} and @var{ya} and only stores the static state
computed from the data.  The @var{xa} data array is always assumed to be
strictly ordered, with increasing @math{x} values; 
the behavior for other arrangements is not defined.
@end deftypefun

@deftypefun void gsl_interp_free (gsl_interp * @var{interp})
This function frees the interpolation object @var{interp}.
@end deftypefun

@node 1D Interpolation Types
@section 1D Interpolation Types
@tindex gsl_interp_type
The interpolation library provides the following interpolation types:

@deffn {Interpolation Type} gsl_interp_linear
@cindex linear interpolation
Linear interpolation.  This interpolation method does not require any
additional memory.
@end deffn

@deffn {Interpolation Type} gsl_interp_polynomial
@cindex polynomial interpolation
Polynomial interpolation.  This method should only be used for
interpolating small numbers of points because polynomial interpolation
introduces large oscillations, even for well-behaved datasets.  The
number of terms in the interpolating polynomial is equal to the number
of points.
@end deffn

@deffn {Interpolation Type} gsl_interp_cspline
@cindex cubic splines
Cubic spline with natural boundary conditions.  The resulting curve is
piecewise cubic on each interval, with matching first and second
derivatives at the supplied data-points.  The second derivative is
chosen to be zero at the first point and last point.
@end deffn

@deffn {Interpolation Type} gsl_interp_cspline_periodic
Cubic spline with periodic boundary conditions.  The resulting curve
is piecewise cubic on each interval, with matching first and second
derivatives at the supplied data-points.  The derivatives at the first
and last points are also matched.  Note that the last point in the
data must have the same y-value as the first point, otherwise the
resulting periodic interpolation will have a discontinuity at the
boundary.

@end deffn

@deffn {Interpolation Type} gsl_interp_akima
@cindex Akima splines
Non-rounded Akima spline with natural boundary conditions.  This method
uses the non-rounded corner algorithm of Wodicka.
@end deffn

@deffn {Interpolation Type} gsl_interp_akima_periodic
Non-rounded Akima spline with periodic boundary conditions.  This method
uses the non-rounded corner algorithm of Wodicka.
@end deffn

@deffn {Interpolation Type} gsl_interp_steffen
Steffen's method guarantees the monotonicity of the interpolating function
between the given data points. Therefore, minima and maxima can only occur
exactly at the data points, and there can never be spurious oscillations
between data points. The interpolated function is piecewise cubic
in each interval. The resulting curve and its first derivative
are guaranteed to be continuous, but the second derivative may be
discontinuous.
@end deffn

@noindent
The following related functions are available:

@deftypefun {const char *} gsl_interp_name (const gsl_interp * @var{interp})
This function returns the name of the interpolation type used by @var{interp}.
For example,

@example
printf ("interp uses '%s' interpolation.\n", 
        gsl_interp_name (interp));
@end example

@noindent
would print something like,

@example
interp uses 'cspline' interpolation.
@end example
@end deftypefun

@deftypefun {unsigned int} gsl_interp_min_size (const gsl_interp * @var{interp})
@deftypefunx {unsigned int} gsl_interp_type_min_size (const gsl_interp_type * @var{T})
These functions return the minimum number of points required by the
interpolation object @var{interp} or interpolation type @var{T}.  For
example, Akima spline interpolation requires a minimum of 5 points.
@end deftypefun



@node 1D Index Look-up and Acceleration
@section 1D Index Look-up and Acceleration

The state of searches can be stored in a @code{gsl_interp_accel} object,
which is a kind of iterator for interpolation lookups.  It caches the
previous value of an index lookup.  When the subsequent interpolation
point falls in the same interval its index value can be returned
immediately.

@deftypefun size_t gsl_interp_bsearch (const double @var{x_array}[], double @var{x}, size_t @var{index_lo}, size_t @var{index_hi})
This function returns the index @math{i} of the array @var{x_array} such
that @code{x_array[i] <= x < x_array[i+1]}.  The index is searched for
in the range [@var{index_lo},@var{index_hi}].  @inlinefn{}
@end deftypefun

@deftypefun {gsl_interp_accel *} gsl_interp_accel_alloc (void)
@tindex gsl_interp_accel
This function returns a pointer to an accelerator object, which is a
kind of iterator for interpolation lookups.  It tracks the state of
lookups, thus allowing for application of various acceleration
strategies.
@end deftypefun

@deftypefun size_t gsl_interp_accel_find (gsl_interp_accel * @var{a}, const double @var{x_array}[], size_t @var{size}, double @var{x})
This function performs a lookup action on the data array @var{x_array}
of size @var{size}, using the given accelerator @var{a}.  This is how
lookups are performed during evaluation of an interpolation.  The
function returns an index @math{i} such that @code{x_array[i] <= x <
x_array[i+1]}.  @inlinefn{}
@end deftypefun

@deftypefun int gsl_interp_accel_reset (gsl_interp_accel * @var{acc});
This function reinitializes the accelerator object @var{acc}.  It
should be used when the cached information is no longer
applicable---for example, when switching to a new dataset.
@end deftypefun

@deftypefun void gsl_interp_accel_free (gsl_interp_accel* @var{acc})
This function frees the accelerator object @var{acc}.
@end deftypefun

@node 1D Evaluation of Interpolating Functions
@section 1D Evaluation of Interpolating Functions

@deftypefun  double gsl_interp_eval (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_interp_eval_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}, double * @var{y})
These functions return the interpolated value of @var{y} for a given
point @var{x}, using the interpolation object @var{interp}, data
arrays @var{xa} and @var{ya} and the accelerator @var{acc}.  When
@var{x} is outside the range of @var{xa}, the error code
@code{GSL_EDOM} is returned with a value of @code{GSL_NAN} for
@var{y}.
@end deftypefun

@deftypefun double gsl_interp_eval_deriv (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_interp_eval_deriv_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d})
These functions return the derivative @var{d} of an interpolated
function for a given point @var{x}, using the interpolation object
@var{interp}, data arrays @var{xa} and @var{ya} and the accelerator
@var{acc}. 
@end deftypefun

@deftypefun double gsl_interp_eval_deriv2 (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_interp_eval_deriv2_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d2})
These functions return the second derivative @var{d2} of an interpolated
function for a given point @var{x}, using the interpolation object
@var{interp}, data arrays @var{xa} and @var{ya} and the accelerator
@var{acc}. 
@end deftypefun

@deftypefun double gsl_interp_eval_integ (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{a}, double @var{b}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_interp_eval_integ_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{a}, double @var{b}, gsl_interp_accel * @var{acc}, double * @var{result})
These functions return the numerical integral @var{result} of an
interpolated function over the range [@var{a}, @var{b}], using the
interpolation object @var{interp}, data arrays @var{xa} and @var{ya} and
the accelerator @var{acc}.
@end deftypefun

@node 1D Higher-level Interface
@section 1D Higher-level Interface

The functions described in the previous sections required the user to
supply pointers to the @math{x} and @math{y} arrays on each call.  The
following functions are equivalent to the corresponding
@code{gsl_interp} functions but maintain a copy of this data in the
@code{gsl_spline} object.  This removes the need to pass both @var{xa}
and @var{ya} as arguments on each evaluation. These functions are
defined in the header file @file{gsl_spline.h}.

@deftypefun {gsl_spline *} gsl_spline_alloc (const gsl_interp_type * @var{T}, size_t @var{size})
@tindex gsl_spline
@end deftypefun

@deftypefun int gsl_spline_init (gsl_spline * @var{spline}, const double @var{xa}[], const double @var{ya}[], size_t @var{size})
@end deftypefun

@deftypefun void gsl_spline_free (gsl_spline * @var{spline})
@end deftypefun

@deftypefun {const char *} gsl_spline_name (const gsl_spline * @var{spline})
@end deftypefun

@deftypefun {unsigned int} gsl_spline_min_size (const gsl_spline * @var{spline})
@end deftypefun

@deftypefun double gsl_spline_eval (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_spline_eval_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{y})
@end deftypefun

@deftypefun double gsl_spline_eval_deriv (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_spline_eval_deriv_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d})
@end deftypefun

@deftypefun double gsl_spline_eval_deriv2 (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_spline_eval_deriv2_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d2})
@end deftypefun

@deftypefun double gsl_spline_eval_integ (const gsl_spline * @var{spline}, double @var{a}, double @var{b}, gsl_interp_accel * @var{acc})
@deftypefunx int gsl_spline_eval_integ_e (const gsl_spline * @var{spline}, double @var{a}, double @var{b}, gsl_interp_accel * @var{acc}, double * @var{result})
@end deftypefun

@node 1D Interpolation Example programs
@section Examples of 1D Interpolation

The following program demonstrates the use of the interpolation and
spline functions.  It computes a cubic spline interpolation of the
10-point dataset @math{(x_i, y_i)} where @math{x_i = i + \sin(i)/2} and
@math{y_i = i + \cos(i^2)} for @math{i = 0 \dots 9}.

@example
@verbatiminclude examples/interp.c
@end example

@noindent
The output is designed to be used with the @sc{gnu} plotutils
@code{graph} program,

@example
$ ./a.out > interp.dat
$ graph -T ps < interp.dat > interp.ps
@end example

@iftex
@sp 1
@center @image{interp2,3.4in}
@end iftex

@noindent
The result shows a smooth interpolation of the original points.  The
interpolation method can be changed simply by varying the first argument of
@code{gsl_spline_alloc}.

The next program demonstrates a periodic cubic spline with 4 data
points.  Note that the first and last points must be supplied with 
the same y-value for a periodic spline.

@example
@verbatiminclude examples/interpp.c
@end example

@noindent

The output can be plotted with @sc{gnu} @code{graph}.

@example
$ ./a.out > interp.dat
$ graph -T ps < interp.dat > interp.ps
@end example

@iftex
@sp 1
@center @image{interpp2,3.4in}
@end iftex

@noindent
The result shows a periodic interpolation of the original points. The
slope of the fitted curve is the same at the beginning and end of the
data, and the second derivative is also.

The next program illustrates the difference between the cubic spline,
Akima, and Steffen interpolation types on a difficult dataset.

@example
@verbatiminclude examples/interp_compare.c
@end example

@iftex
@sp 1
@center @image{interp_compare,6in}
@end iftex

@noindent
The cubic method exhibits a local maxima between the 6th and 7th data points
and continues oscillating for the rest of the data. Akima also shows a
local maxima but recovers and follows the data well after the 7th grid point.
Steffen preserves monotonicity in all intervals and does not exhibit oscillations,
at the expense of having a discontinuous second derivative.

@node 1D Interpolation References and Further Reading
@section References and Further Reading

Descriptions of the interpolation algorithms and further references can
be found in the following publications:

@itemize @w{}
@item C.W. Ueberhuber,
@cite{Numerical Computation (Volume 1), Chapter 9 ``Interpolation''},
Springer (1997), ISBN 3-540-62058-3.

@item D.M. Young, R.T. Gregory,
@cite{A Survey of Numerical Mathematics (Volume 1), Chapter 6.8},
Dover (1988), ISBN 0-486-65691-8.

@item M. Steffen,
@cite{A simple method for monotonic interpolation in one dimension},
Astron. Astrophys. 239, 443-450, 1990.
@end itemize

@node 2D Introduction to Interpolation
@section Introduction to 2D Interpolation

Given a set of @math{x} coordinates @math{x_1,...,x_m} and a set of
@math{y} coordinates @math{y_1,...,y_n}, each in increasing order,
plus a set of function values @math{z_{ij}}
for each grid point @math{(x_i,y_j)}, the routines described in this
section compute a continuous interpolation function @math{z(x,y)} such
that @math{z(x_i,y_j) = z_{ij}}.

@node 2D Interpolation Functions
@section 2D Interpolation Functions

The interpolation function for a given dataset is stored in a
@code{gsl_interp2d} object. These are created by the following functions.

@deftypefun {gsl_interp2d *} gsl_interp2d_alloc (const gsl_interp2d_type * @var{T}, const size_t @var{xsize}, const size_t @var{ysize})
@tindex gsl_interp
This function returns a pointer to a newly allocated interpolation
object of type @var{T} for @var{xsize} grid points in the @math{x}
direction and @var{ysize} grid points in the @math{y} direction.
@end deftypefun

@deftypefun int gsl_interp2d_init (gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const size_t @var{xsize}, const size_t @var{ysize})
This function initializes the interpolation object @var{interp} for the
data (@var{xa},@var{ya},@var{za}) where @var{xa} and @var{ya} are arrays of
the @math{x} and @math{y} grid points of size @var{xsize} and @var{ysize}
respectively, and @var{za} is an array of function values of size
@var{xsize}*@var{ysize}.  The interpolation object (@code{gsl_interp2d}) does
not save the data arrays @var{xa}, @var{ya}, and @var{za} and only stores the
static state computed from the data. The @var{xa} and @var{ya} data arrays
are always assumed to be strictly ordered, with increasing @math{x,y} values; 
the behavior for other arrangements is not defined.
@end deftypefun

@deftypefun void gsl_interp2d_free (gsl_interp2d * @var{interp})
This function frees the interpolation object @var{interp}.
@end deftypefun

@node 2D Interpolation Grids
@section 2D Interpolation Grids

The 2D interpolation routines access the function values @math{z_{ij}}
with the following ordering:
@tex
\beforedisplay
$$
z_{ij} = za[j*xsize + i]
$$
\afterdisplay
@end tex
@ifinfo

@example
z_ij = za[j*xsize + i]
@end example

@end ifinfo
with @math{i = 0,...,xsize-1} and @math{j = 0,...,ysize-1}. However,
for ease of use, the following functions are provided to add and retrieve
elements from the function grid without requiring knowledge of the
internal ordering.

@deftypefun int gsl_interp2d_set (const gsl_interp2d * @var{interp}, double @var{za}[], const size_t @var{i}, const size_t @var{j}, const double @var{z})
This function sets the value @math{z_{ij}} for grid point
(@var{i},@var{j}) of the array @var{za} to @var{z}.
@end deftypefun

@deftypefun double gsl_interp2d_get (const gsl_interp2d * @var{interp}, const double @var{za}[], const size_t @var{i}, const size_t @var{j})
This function returns the value @math{z_{ij}} for grid point
(@var{i},@var{j}) stored in the array @var{za}.
@end deftypefun

@deftypefun size_t gsl_interp2d_idx (const gsl_interp2d * @var{interp}, const size_t @var{i}, const size_t @var{j})
This function returns the index corresponding to the grid point
(@var{i},@var{j}). The index is given by @math{j*xsize + i}.
@end deftypefun

@node 2D Interpolation Types
@section 2D Interpolation Types
@tindex gsl_interp2d_type
The interpolation library provides the following 2D interpolation types:

@deffn {Interpolation Type} gsl_interp2d_bilinear
@cindex bilinear interpolation
Bilinear interpolation.  This interpolation method does not require any
additional memory.
@end deffn

@deffn {Interpolation Type} gsl_interp2d_bicubic
@cindex bicubic interpolation
Bicubic interpolation.
@end deffn

@deftypefun {const char *} gsl_interp2d_name (const gsl_interp2d * @var{interp})
This function returns the name of the interpolation type used by @var{interp}.
For example,

@example
printf ("interp uses '%s' interpolation.\n", 
        gsl_interp2d_name (interp));
@end example

@noindent
would print something like,

@example
interp uses 'bilinear' interpolation.
@end example
@end deftypefun

@deftypefun {unsigned int} gsl_interp2d_min_size (const gsl_interp2d * @var{interp})
@deftypefunx {unsigned int} gsl_interp2d_type_min_size (const gsl_interp2d_type * @var{T})
These functions return the minimum number of points required by the
interpolation object @var{interp} or interpolation type @var{T}.  For
example, bicubic interpolation requires a minimum of 4 points.
@end deftypefun

@node 2D Evaluation of Interpolating Functions
@section 2D Evaluation of Interpolating Functions

@deftypefun double gsl_interp2d_eval (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{z})
These functions return the interpolated value of @var{z} for a given
point (@var{x},@var{y}), using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}.  When @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, the error code
@code{GSL_EDOM} is returned.
@end deftypefun

@deftypefun double gsl_interp2d_eval_extrap (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_extrap_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{z})
These functions return the interpolated value of @var{z} for a given
point (@var{x},@var{y}), using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}. The functions perform no bounds checking, so
when @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, extrapolation is performed.
@end deftypefun

@deftypefun double gsl_interp2d_eval_deriv_x (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_deriv_x_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
These functions return the interpolated value @var{d}
@math{= \partial z / \partial x} for a given point (@var{x},@var{y}),
using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}.  When @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, the error code
@code{GSL_EDOM} is returned.
@end deftypefun

@deftypefun double gsl_interp2d_eval_deriv_y (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_deriv_y_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
These functions return the interpolated value @var{d}
@math{= \partial z / \partial y} for a given point (@var{x},@var{y}),
using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}.  When @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, the error code
@code{GSL_EDOM} is returned.
@end deftypefun

@deftypefun double gsl_interp2d_eval_deriv_xx (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_deriv_xx_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
These functions return the interpolated value @var{d}
@math{= \partial^2 z / \partial x^2} for a given point (@var{x},@var{y}),
using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}.  When @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, the error code
@code{GSL_EDOM} is returned.
@end deftypefun

@deftypefun double gsl_interp2d_eval_deriv_yy (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_deriv_yy_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
These functions return the interpolated value @var{d}
@math{= \partial^2 z / \partial y^2} for a given point (@var{x},@var{y}),
using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}.  When @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, the error code
@code{GSL_EDOM} is returned.
@end deftypefun

@deftypefun double gsl_interp2d_eval_deriv_xy (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_interp2d_eval_deriv_xy_e (const gsl_interp2d * @var{interp}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
These functions return the interpolated value @var{d}
@math{= \partial^2 z / \partial x \partial y} for a given point (@var{x},@var{y}),
using the interpolation object @var{interp}, data
arrays @var{xa}, @var{ya}, and @var{za} and the accelerators @var{xacc}
and @var{yacc}.  When @var{x} is outside the range of @var{xa} or @var{y}
is outside the range of @var{ya}, the error code
@code{GSL_EDOM} is returned.
@end deftypefun

@node 2D Higher-level Interface
@section 2D Higher-level Interface

The functions described in the previous sections required the user to
supply pointers to the @math{x}, @math{y}, and @math{z} arrays on each call.
The following functions are equivalent to the corresponding
@code{gsl_interp2d} functions but maintain a copy of this data in the
@code{gsl_spline2d} object.  This removes the need to pass @var{xa},
@var{ya}, and @var{za} as arguments on each evaluation. These functions are
defined in the header file @file{gsl_spline2d.h}.

@deftypefun {gsl_spline2d *} gsl_spline2d_alloc (const gsl_interp2d_type * @var{T}, size_t @var{xsize}, size_t @var{ysize})
@tindex gsl_spline2d
@end deftypefun

@deftypefun int gsl_spline2d_init (gsl_spline2d * @var{spline}, const double @var{xa}[], const double @var{ya}[], const double @var{za}[], size_t @var{xsize}, size_t @var{ysize})
@end deftypefun

@deftypefun void gsl_spline2d_free (gsl_spline2d * @var{spline})
@end deftypefun

@deftypefun {const char *} gsl_spline2d_name (const gsl_spline2d * @var{spline})
@end deftypefun

@deftypefun {unsigned int} gsl_spline2d_min_size (const gsl_spline2d * @var{spline})
@end deftypefun

@deftypefun double gsl_spline2d_eval (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_spline2d_eval_e (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{z})
@end deftypefun

@deftypefun double gsl_spline2d_eval_deriv_x (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_spline2d_eval_deriv_x_e (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
@end deftypefun

@deftypefun double gsl_spline2d_eval_deriv_y (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_spline2d_eval_deriv_y_e (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
@end deftypefun

@deftypefun double gsl_spline2d_eval_deriv_xx (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_spline2d_eval_deriv_xx_e (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
@end deftypefun

@deftypefun double gsl_spline2d_eval_deriv_yy (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_spline2d_eval_deriv_yy_e (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
@end deftypefun

@deftypefun double gsl_spline2d_eval_deriv_xy (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc})
@deftypefunx int gsl_spline2d_eval_deriv_xy_e (const gsl_spline2d * @var{spline}, const double @var{x}, const double @var{y}, gsl_interp_accel * @var{xacc}, gsl_interp_accel * @var{yacc}, double * @var{d})
@end deftypefun

@deftypefun int gsl_spline2d_set (const gsl_spline2d * @var{spline}, double @var{za}[], const size_t @var{i}, const size_t @var{j}, const double @var{z})
@end deftypefun

@deftypefun double gsl_spline2d_get (const gsl_spline2d * @var{spline}, const double @var{za}[], const size_t @var{i}, const size_t @var{j})
This function returns the value @math{z_{ij}} for grid point
(@var{i},@var{j}) stored in the array @var{za}.
@end deftypefun

@node 2D Interpolation Example programs
@section 2D Interpolation Example programs

The following example performs bilinear interpolation on the unit
square, using @math{z} values of @math{(0,1,0.5,1)} going clockwise
around the square.

@example
@verbatiminclude examples/interp2d.c
@end example

@noindent
The results of the interpolation are shown in the following plot,
where the corners are labeled with their fixed @math{z} values.

@iftex
@sp 1
@center @image{interp2d,6in}
@end iftex