### 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.

*   Declare variables
INTEGER ID                      ! Specifies use value not differential
INTEGER IFAIL                   ! Was the surface successfully created?
INTEGER MXY                     ! Size of the grid
INTEGER ORD                     ! Order of polynomial used

DOUBLE PRECISION DVALUE         ! Interpolated value returned
DOUBLE PRECISION FV1(8,8)       ! Grid Z values
DOUBLE PRECISION X1(8), Y1(8)   ! X,Y grid locations
DOUBLE PRECISION XD, YD         ! X,Y coord for interpolation

DOUBLE PRECISION BCOEF(8,8)     ! Array used by PDA_DB2INK
DOUBLE PRECISION TX(11),TY(11)  ! Array used by PDA_DB2INK
DOUBLE PRECISION WORK(168)      ! Array used by PDA_DB2INK

*   Subroutine initial error values.
IFAIL=0
STATUS=0

*   Order of polynomial used in this example.
ORD=3

*   Use value rather than differential.
ID=0

*   Size of data grid. 8x8 in this instance.
MXY=8

*   Build the surface fit using the grid contents.
CALL PDA_DB2INK(X1,MXY,Y1,MXY,FV1,MXY,
:                ORD,ORD,TX,TY,BCOEF,
:                WORK,IFAIL,STATUS)

*   If IFAIL=1 then okay to interpolate values.
IF (IFAIL.EQ.1) THEN

*       Setup evaluation routine.
IFAIL=0
CALL PDA_DB2VAL(XD,YD,ID,ID,TX,TY,
:                    MXY,MXY,ORD,ORD,BCOEF,WORK,
:                    DVALUE,IFAIL,STATUS)

END IF

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.

*  Local Variables:
INTEGER  ISTAT                   ! Status
INTEGER  MD                      ! Mode
INTEGER  NDP                     ! Grid points
INTEGER  NCP                     ! Not used
INTEGER  NOP                     ! Size of returned array
INTEGER  IWK(2500)               ! Workspace
REAL     WK(640)                 ! Workspace
REAL     XD(80),YD(80),ZD(80)    ! Surface data
REAL     XI(1),YI(1)             ! Extrap points
REAL     ZI(1,1)                 ! Results

*   Set the error flag default values.
ISTAT=0
STATUS=0

*   Set mode and the grid positions (for PDA_IDBVIP).
MD= 1
NCP=2
NOP=1

*   Call interpolation subroutine.
CALL PDA_IDBVIP(MD,NCP,NDP,XD,YD,ZD,NOP,XI,YI,
:                ZI,IWK,WK,ISTAT,STATUS)

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:

INTEGER  ISTAT                   ! Status
INTEGER  MD                      ! Mode
INTEGER  NDP                     ! Grid points
INTEGER  NCP                     ! Not used
INTEGER  NXI,NYI                 ! Output gird size
INTEGER  IWK(2500)               ! Workspace
REAL     XD(80),YD(80),ZD(80)    ! Surface data
REAL     XI(1),YI(1)             ! Extrap points
REAL     ZI(1,1)                 ! Results
REAL     WK(640)                 ! Workspace

*   Set the error flag default values.
ISTAT=0
STATUS=0

*   Set up the grid positions.
NXI=1
NYI=1

*   Set mode and the grid positions (for PDA_IDBVIP).
MD= 1
NCP=2

*         Call the surface fitting subroutine.
CALL PDA_IDSFFT(MD,NCP,NDP,XD,YD,ZD,
:                      NXI,NYI,XI,YI,ZI,IWK,
:                      WK,ISTAT,STATUS)

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.