Index: /issm/trunk-jpl/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 14048)
+++ /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 14049)
@@ -426,4 +426,8 @@
 /*FUNCTION Observations::InterpolationV4{{{*/
 void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
+	/* Reference:  David T. Sandwell, Biharmonic spline interpolation of GEOS-3
+	 * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
+	 * 1987.  Describes interpolation using value or gradient of value in any
+	 * dimension.*/
 
 	/*Intermediaries*/
@@ -435,5 +439,5 @@
 	IssmPDouble *Green   = NULL;
 	IssmPDouble *weights = NULL;
-	IssmPDouble *d       = NULL;
+	IssmPDouble *g       = NULL;
 
 	/*Some checks*/
@@ -455,9 +459,15 @@
 		/*Allocate intermediary matrix and vectors*/
 		Green = xNew<IssmPDouble>(n_obs*n_obs);
-		d     = xNew<IssmPDouble>(n_obs);
+		g     = xNew<IssmPDouble>(n_obs);
 
 		/*First: distance vector*/
 		for(i=0;i<n_obs;i++){
-			d[i] = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
+			h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
+			if(h>0){
+				g[i] = h*h*(log(h)-1.);
+			}
+			else{
+				g[i] = 0.;
+			}
 		}
 
@@ -483,7 +493,8 @@
 #endif
 
-		/*ok now we can interpolate with these weights*/
-
-		_error_("Not implemented yet");
+		/*Interpolate*/
+		prediction = 0;
+		for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
+
 	}
 
@@ -493,4 +504,7 @@
 	xDelete<IssmPDouble>(y);
 	xDelete<IssmPDouble>(obs);
+	xDelete<IssmPDouble>(Green);
+	xDelete<IssmPDouble>(g);
+	xDelete<IssmPDouble>(weights);
 }/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/