+ int i;
+ gint16 elevs[4], dists[4];
+ gint16 max_dist;
+ gdouble t = 0.0;
+ gdouble b = 0.0;
+
+ if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
+ return VIK_DEM_INVALID_ELEVATION;
+
+ max_dist = 0;
+ for (i = 0; i < 4; i++) {
+ if (dists[i] < 1) {
+ return(elevs[i]);
+ }
+ if (dists[i] > max_dist)
+ max_dist = dists[i];
+ }
+
+ gdouble tmp;
+#if 0 /* derived method by Franke & Nielson. Does not seem to work too well here */
+ for (i = 0; i < 4; i++) {
+ tmp = pow((1.0*(max_dist - dists[i])/max_dist*dists[i]), 2);
+ t += tmp*elevs[i];
+ b += tmp;
+ }
+#endif
+
+ for (i = 0; i < 4; i++) {
+ tmp = pow((1.0/dists[i]), 2);
+ t += tmp*elevs[i];
+ b += tmp;
+ }
+
+ // fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);
+
+ return(t/b);
+