where R is a user-specified tension parameter and s_{k} is
the angular distance
between the interpolation point and the kth observation value.
The weights A_{k} are computed so that the estimated function
agrees with the
observations at the observation points.
For n observation points, this requires the solution of n simultaneous
equations and the inversion
of an nxn matrix.
As a result, a multiquadric fit is limited to a reasonable number of points.
In Spherekit, if more than 500 points are entered, the domain is broken
up into overlapping
regions, and a surface fit is performed on each region. Points lying in
more than one region are
assigned a value that is a weighted average of the multiple estimates.

An option is available to use reciprocal multiquadric interpolation
(RMQ). In this case, the square root in the above equation is replaced
by the reciprocal of the square root.

### References

Hardy, R. L. and W. M. Gofert (1975), Least squares prediction of
gravity anomalies, geoidal
undulations, and deflections of the vertical multiquadric harmonic
functions, Geophysical
Research Letters, 2, 423-426.

Pottmann, H. and M. Eck (1990), Modified multiquadric methods for
scattered data
interpolation over a sphere, Computer Aided Design, 7, 313-321.