27. 
where C is a prefixed constant and p is either % or 1. To avoid sharp peaks and 
dips at the reference points, one must not use the value C-0 with p-X. Chiles and 
Delfiner [23], who call their method kriging, do not specify the form of the 
distance function. 
It is important to note that, once the heights have been reduced to a 
proper trend surface, there will be a minimum distance for which they tend to have 
opposite signs. In the case of regularly undulating terrain that distance will be 
half the wavelength of the undulations. At that distance, the correlation function 
has a negative value. As a result, this function has the appearance of a rather 
rapidly damped oscillation. None of the above distance functions has this feature. 
Because of the above limitations, one cannot justifiably claim that 
according to theory one of the above described methods is the most accurate one. 
Still, in one experiment [9], the best results obtained with summation methods 
proved to be slightly better than the best results obtained with moving surface 
methods. In that experiment, a smooth undulating reference surface was used. 
It is noteworthy also that in that experiment the best results were 
obtained when the distance function and the surfaces of revolution were so flat 
that even with double-precision arithmetic the matrix B was nearly singular. 
Similarly, Lauer [52] obtained optimum results with a correlation function which, 
when written in the form of eq. (4c), has a very small value of the coefficient a. 
In some programs, the inversion of a very large matrix B in eq. (3) is 
avoided by restricting the number of reference points that contributes to this 
matrix. For the Stuttgart Contour Program, the interpolation area which may 
contain several thousand reference heights is subdivided into overlapping units 
with an average of 70 reference points each. Koch gives examples in which the 
summation method is used with from 12 to 16 and in one case 20 points. (The moving 
surface method is used with from 8 to 12 points) A matrix B must here be inverted 
for every combination of reference points. Lauer [53] subdivides the interpolation 
area into triangles with vertices at the reference points. The interpolation of 
all points in a triangle makes use of the same matrix B, computed from the heights 
at the vertices of the triangle and at the immediately surrounding ones. These 
programs are suitable especially to compute heights at the nodes of a grid, in 
which case a supplementary interpolation will prevent discontinuities. 
5. Simultaneous patchwise polynomials 
  
A third group of interpolation methods starts by dividing the region of 
interest into square or rectangular elements by means of a regular grid and 
represents the terrain surface in each element by a low-degree polynomial in such 
à way that the total surface is continuous and possibly also smooth. In other 
words, at their common boundaries the local polynomial surfaces agree in height and 
possibly in tilt. The common characteristic of these methods is the simultaneous 
computation of all the local surfaces. Three such methods have been described in 
some detail. 
Bosman, Eckhart, and Kubik [50], [22] specify that the polynomials must 
satisfy the condition that the weighted sum of squares of the residuals at the 
reference points plus a weighted sum which expresses the flatness of the total 
surface must be a minimum. This condition leads to the formation and solution of 
normal equations by the method of least squares. The second one of the two sums 
is introduced to ensure non-singularity of the matrix of the normal equations and 
to flatten out the surfaces in poorly controlled areas. The two weights and the 
size of the grid can be varied to reduce the residuals to the size of the measur- 
ing errors. Each local polynomial is either a full 16-term bicubic polynomial or 
it is composed of two linear parts, each serving one of two triangles into which 
a grid element is subdivided.