Changeset 14049

Timestamp:

11/29/12 08:09:39 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: Added minimum curvature interpolation

File:

• issm/trunk-jpl/src/c/Container/Observations.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/Container/Observations.cpp

@@ -426 +426 @@
 /*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 +439 @@
         IssmPDouble *Green   = NULL;
         IssmPDouble *weights = NULL;
-        IssmPDouble *d       = NULL;
+        IssmPDouble *g       = NULL;
 
         /*Some checks*/
@@ -455 +459 @@
                 /*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 +493 @@
 #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 +504 @@
         xDelete<IssmPDouble>(y);
         xDelete<IssmPDouble>(obs);
+        xDelete<IssmPDouble>(Green);
+        xDelete<IssmPDouble>(g);
+        xDelete<IssmPDouble>(weights);
 }/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/