10 Two-dimensional Interpolation and Fitting

The routines for this sort of application are:

• PDA_CHE2D PDA_CHE2D evaluates a 2-dimensional Chebyshev polynomial. A single precision version (PDA_CHE2R) is available.
• PDA_BISPEV PDA_BISPEV evaluates the bivariate spline approximation found by PDA_SURFIT.
• PDA_DB2INK PDA_DB2INK determines a piecewise polynomial function that interpolates the two-dimensional gridded data. Users specify the polynomial order (degree+1) of the interpolant and (optionally) the knot sequence.

  The interpolating function is a piecewise polynomial represented as a tensor product of one-dimensional B-splines.

• PDA_DB2VAL PDA_DB2VAL evaluates the tensor product piecewise polynomial interpolant constructed by the routine PDA_DB2INK, or, alternatively evaluates one of its derivatives, at a given point. Function values returned are double precision.
• PDA_IDBVIP PDA_IDBVIP performs bivariate interpolation when the the data points are irregularly distributed in the x-y plane. Function values returned are single precision.
• PDA_IDSFFT PDA_IDSFFT performs smooth surface fitting when the data points are irregularly distributed in the x-y plane. Function values returned are single precision.
• PDA_SURFIT PDA_SURFIT determines a smooth bivariate spline approximation for irregularly distributed data.

E02DAF calculates a 2-D bi-cubic spline interpolating surface for points from a regular grid. The routines E02DEF/E02DFF may then be used to compute values of the spline at the required location. This approach has has been replaced by using the routine pair PDA_DB2INK and PDA_DB2VAL. No significant change in program performance was found.

E02SAF provides a 2-D surface fit for data on an irregular spaced grid. The method employed is that of Renka and Cline where the grid is used to construct a set of suitably weighed equiangular triangles that (with appropriate weighting) describe the surface. The interpolated value of the gridded data at any point within the grid is extracted by E02SBF. This pair of routines has been replaced using PDA_IDBVIP with no significant change in program performance. During incomplete trials the routine PDA_IDSFFT also appeared to give sensible results.

10.1 Replacing calls to E02DAF

E02DAF can be replaced with the following code.

In the example shown the value of the surface fitted is calculated at a point defined by XD and YD. The surface constructed is a 3rd order, polynomial and the value returned by PDA_DB2VAL is the surface value not its differential. The example assumes that the array FV1() already contains values for the data at each of the grid points defined in arrays X1 and Y1. The size of the work arrays is determined by the size of the grid required. The variable DVALUE contains the interpolated value.

10.2 Replacing calls to E02SAF

It has been found in tests that the following code examples adequately replace calls to E02SAF.

In the example shown, the value of the interpolated surface at XI(1),YI(1) is returned in the variable ZI(1,1). The surface is constructed using co-ordinate information from the arrays of XD() and YD(), and surface values contained in the array ZD(). The number of data points available to define the surface is NDP.

If several values are required from locations within the irregular grid, the mode variable MD should for the first point be set to 1 but may subsequently be 2. This significantly increases the execution speed. The size of the work arrays is determined by the size of the grid required.

Alternatively, the routine PDA_IDSFFT may be used thus:

In the example shown the surface is constructed from the irregular grid information contained in arrays XD()/YD() (location) and ZD() (surface value). Interpolated values calculated from the fitted surface are returned in the ZI() array at the co-ordinates specified by the data in the XI() and YI() arrays.

If several values are required from locations within the irregular grid, the mode variable MD should for the first point be set to 1 but may subsequently be 2. This significantly increases the execution speed. The size of the work arrays is determined by the size of the grid required.

10.3 Replacing calls to E02CBF

The routine PDA_CHE2D (or its single precision equivalent PDA_CHE2R) can be used to replace E02CBF. The replacement is straightforward with most of the arguments being the same, albeit in a slightly different order.