This is an implementation of a bilinear interpolating function.
For a point (x0,y0) contained in a rectangle (x1,y1),(x2,y1),
(x2,y2),(x1,y2) and x1<x2, y1<y2, the first step is to get z()
at locations (x0,y1) and (x0,y2) as convex linear combinations
z(x0,y*)=a*z(x1,y*)+(1-a)*z(x2,y*) where a=(x2-x1)/(x0-x1) for
y*=y1,y2. In a second step z(x0,y0) is calculated as convex linear
combination between z(x0,y1) and z(x0,y2) as
z(x0,y1)=b*z(x0,y1)+(1-b)*z(x0,y2) where b=(y2-y1)/(y0-y1).
Finally, z(x0,y0) is a convex linear combination of the z values at
the corners of the containing rectangle with weights according to
the distance from (x0,y0) to these corners.
The grid lines can be unevenly spaced.