SURF
Computes a smooth bivariate interpolant to scattered data that is locally a quintic polynomial in two variables.
Required Arguments
XYDATA — A 2 by NDATA array containing the coordinates of the interpolation points. (Input)
These points must be distinct. The x-coordinate of the I-th data point is stored in XYDATA(1, I) and the y-coordinate of the I-th data point is stored in XYDATA(2, I).
FDATA — Array of length NDATA containing the interpolation values. (Input) FDATA(I) contains the value at (XYDATA(1, I), XYDATA(2, I)).
XOUT — Array of length NXOUT containing an increasing sequence of points. (Input)
These points are the x-coordinates of a grid on which the interpolated surface is to be evaluated.
YOUT — Array of length NYOUT containing an increasing sequence of points. (Input)
These points are the y-coordinates of a grid on which the interpolated surface is to be evaluated.
SUR — Matrix of size NXOUT by NYOUT. (Output)
This matrix contains the values of the surface on the XOUT by YOUT grid, i.e. SUR(I, J) contains the interpolated value at (XOUT(I), YOUT(J)).
Optional Arguments
NDATA — Number of data points. (Input)
NDATA must be at least four.
Default: NDATA = size (FDATA,1).
NXOUT — The number of elements in XOUT. (Input)
Default: NXOUT = size (XOUT,1).
NYOUT — The number of elements in YOUT. (Input)
Default: NYOUT = size (YOUT,1).
LDSUR — Leading dimension of SUR exactly as specified in the dimension statement of the calling program. (Input)
LDSUR must be at least as large as NXOUT.
Default: LDSUR = size (SUR,1).
FORTRAN 90 Interface
Generic: CALL SURF (XYDATA, FDATA, XOUT, YOUT, SUR [, …])
Specific: The specific interface names are S_SURF and D_SURF.
FORTRAN 77 Interface
Single: CALL SURF (NDATA, XYDATA, FDATA, NXOUT, NYOUT, XOUT, YOUT, SUR, LDSUR)
Double: The double precision name is DSURF.
Description
This routine is designed to compute a C 1 interpolant to scattered data in the plane. Given the data points
SURF returns (in SUR, the user-specified grid) the values of the interpolant s. The computation of s is as follows: First the Delaunay triangulation of the points
is computed. On each triangle T in this triangulation, s has the form
Thus, s is a bivariate quintic polynomial on each triangle of the triangulation. In addition, we have
s(xi, yi) = fi for i = 1, …, N
and s is continuously differentiable across the boundaries of neighboring triangles. These conditions do not exhaust the freedom implied by the above representation. This additional freedom is exploited in an attempt to produce an interpolant that is faithful to the global shape properties implied by the data. For more information on this routine, we refer the reader to the article by Akima (1978). The grid is specified by the two integer variables NXOUT, NYOUT that represent, respectively, the number of grid points in the first (second) variable and by two real vectors that represent, respectively, the first (second) coordinates of the grid.
Comments
1. Workspace may be explicitly provided, if desired, by use of S2RF/DS2RF. The reference is:
CALL S2RF (NDATA, XYDATA, FDATA, NXOUT, NYOUT, XOUT, YOUT, SUR, LDSUR, IWK, WK)
The additional arguments are as follows:
IWK — Work array of length 31 * NDATA + 2*(NXOUT * NYOUT).
WK — Work array of length 6 * NDATA.
2. Informational errors
|
Type |
Code |
Description |
|
4 |
5 |
The data point values must be distinct. |
|
4 |
6 |
The XOUT values must be strictly increasing. |
|
4 |
7 |
The YOUT values must be strictly increasing. |
3. This method of interpolation reproduces linear functions.
Example
In this example, the interpolant to the linear function 3 + 7x + 2y is computed from 20 data points equally spaced on the circle of radius 3. We then print the values on a 3 × 3 grid.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER LDSUR, NDATA, NXOUT, NYOUT
PARAMETER (NDATA=20, NXOUT=3, NYOUT=3, LDSUR=NXOUT)
!
INTEGER I, J, NOUT
REAL ABS, COS, F, FDATA(NDATA), FLOAT, PI,&
SIN, SUR(LDSUR,NYOUT), X, XOUT(NXOUT),&
XYDATA(2,NDATA), Y, YOUT(NYOUT)
INTRINSIC ABS, COS, FLOAT, SIN
! Define function
F(X,Y) = 3.0 + 7.0*X + 2.0*Y
! Get value for PI
PI = CONST('PI')
! Set up X, Y, and F data on a circle
DO 10 I=1, NDATA
XYDATA(1,I) = 3.0*SIN(2.0*PI*FLOAT(I-1)/FLOAT(NDATA))
XYDATA(2,I) = 3.0*COS(2.0*PI*FLOAT(I-1)/FLOAT(NDATA))
FDATA(I) = F(XYDATA(1,I),XYDATA(2,I))
10 CONTINUE
! Set up XOUT and YOUT data on [0,1] by
! [0,1] grid.
DO 20 I=1, NXOUT
XOUT(I) = FLOAT(I-1)/FLOAT(NXOUT-1)
20 CONTINUE
DO 30 I=1, NYOUT
YOUT(I) = FLOAT(I-1)/FLOAT(NYOUT-1)
30 CONTINUE
! Interpolate scattered data
CALL SURF (XYDATA, FDATA, XOUT, YOUT, SUR)
! Get output unit number
CALL UMACH (2, NOUT)
! Write heading
WRITE (NOUT,99998)
! Print results
DO 40 I=1, NYOUT
DO 40 J=1, NXOUT
WRITE (NOUT,99999) XOUT(J), YOUT(I), SUR(J,I),&
F(XOUT(J),YOUT(I)),&
ABS(SUR(J,I)-F(XOUT(J),YOUT(I)))
40 CONTINUE
99998 FORMAT (' ', 10X, 'X', 11X, 'Y', 9X, 'SURF', 6X, 'F(X,Y)', 7X,&
'ERROR', /)
99999 FORMAT (1X, 5F12.4)
END
X Y SURF F(X,Y) ERROR
0.0000 0.0000 3.0000 3.0000 0.0000
0.5000 0.0000 6.5000 6.5000 0.0000
1.0000 0.0000 10.0000 10.0000 0.0000
0.0000 0.5000 4.0000 4.0000 0.0000
0.5000 0.5000 7.5000 7.5000 0.0000
1.0000 0.5000 11.0000 11.0000 0.0000
0.0000 1.0000 5.0000 5.0000 0.0000
0.5000 1.0000 8.5000 8.5000 0.0000
1.0000 1.0000 12.0000 12.0000 0.0000
