Index: /issm/trunk-jpl/src/c/Container/Elements.h
===================================================================
--- /issm/trunk-jpl/src/c/Container/Elements.h	(revision 12418)
+++ /issm/trunk-jpl/src/c/Container/Elements.h	(revision 12419)
@@ -15,5 +15,4 @@
 class DataSet;
 class Inputs;
-
 
 class Elements: public DataSet{
Index: /issm/trunk-jpl/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 12418)
+++ /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 12419)
@@ -20,4 +20,5 @@
 #include "../shared/shared.h"
 #include "../include/include.h"
+#include "../modules/modules.h"
 #include "../EnumDefinitions/EnumDefinitions.h"
 #include "../io/io.h"
@@ -115,5 +116,54 @@
 
 /*Methods*/
-/*FUNCTION Observations::ObservationList{{{*/
+/*FUNCTION Observations::ClosestObservation{{{*/
+void Observations::ClosestObservation(double *px,double *py,double *pobs,double x_interp,double y_interp,double radius){
+
+	/*Output and Intermediaries*/
+	bool         stop;
+	int          nobs,i,index;
+	double       h2,hmin2,radius2;
+	int         *indices      = NULL;
+	Observation *observation  = NULL;
+
+	/*If radius is not provided or is 0, return all observations*/
+	if(radius==0) radius=this->quadtree->root->length;
+
+	/*Compute radius square*/
+	radius2 = radius*radius;
+
+	/*Find all observations that are in radius*/
+	this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,radius);
+	for (i=0;i<nobs;i++){
+		observation=(Observation*)this->GetObjectByOffset(indices[i]);
+		h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
+
+		if(i==0){
+			hmin2 = h2;
+			index = i;
+		}
+		else{
+			if(h2<hmin2){
+				hmin2 = h2;
+				index = i;
+			}
+		}
+	}  
+
+	/*Assign output pointer*/
+	if(!nobs){
+		*px=UNDEF;
+		*py=UNDEF;
+		*pobs=UNDEF;
+	}
+	else{
+		observation=(Observation*)this->GetObjectByOffset(indices[index]);
+		*px=observation->x;
+		*py=observation->y;
+		*pobs=observation->value;
+	}
+	xfree((void**)&indices);
+
+}/*}}}*/
+/*FUNCTION Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){{{*/
 void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){
 
@@ -130,6 +180,8 @@
 	Observation *observation  = NULL;
 
+	/*If radius is not provided or is 0, return all observations*/
+	if(radius==0) radius=this->quadtree->root->length;
+
 	/*Compute radius square*/
-	if(radius==0) radius=this->quadtree->root->length;
 	radius2 = radius*radius;
 
@@ -196,4 +248,119 @@
 	*pnobs=nobs;
 }/*}}}*/
+/*FUNCTION Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs){{{*/
+void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs){
+
+	/*Output and Intermediaries*/
+	int          nobs;
+	double      *x            = NULL;
+	double      *y            = NULL;
+	double      *obs          = NULL;
+	Observation *observation  = NULL;
+
+	nobs = this->Size();
+
+	if(nobs){
+		x   = (double*)xmalloc(nobs*sizeof(double));
+		y   = (double*)xmalloc(nobs*sizeof(double));
+		obs = (double*)xmalloc(nobs*sizeof(double));
+		for(int i=0;i<this->Size();i++){
+			observation=(Observation*)this->GetObjectByOffset(i);
+			observation->WriteXYObs(&x[i],&y[i],&obs[i]);
+		}
+	}
+
+	/*Assign output pointer*/
+	*px=x;
+	*py=y;
+	*pobs=obs;
+	*pnobs=nobs;
+}/*}}}*/
+/*FUNCTION Observations::InterpolationKriging{{{*/
+void Observations::InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram){
+
+	/*Intermediaries*/
+	int           i,j,n_obs;
+	double        prediction,error;
+	double        numerator,denominator,ratio;
+	double       *x            = NULL;
+	double       *y            = NULL;
+	double       *obs          = NULL;
+	double       *Gamma        = NULL;
+	double       *GinvG0       = NULL;
+	double       *Ginv1        = NULL;
+	double       *GinvZ        = NULL;
+	double       *gamma0       = NULL;
+	double       *ones         = NULL;
+
+	/*Some checks*/
+	_assert_(mindata>0 && maxdata>0);
+	_assert_(pprediction && perror);
+
+	/*If radius is not provided or is 0, return all observations*/
+	if(radius==0) radius=this->quadtree->root->length;
+
+	/*Get list of observations for current point*/
+	this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+	/*If we have less observations than mindata, return UNDEF*/
+	if(n_obs<mindata){
+		*pprediction = -999.0; 
+		*perror      = -999.0; 
+		return;
+	}
+
+	/*Allocate intermediary matrix and vectors*/
+	Gamma  = (double*)xmalloc(n_obs*n_obs*sizeof(double));
+	gamma0 = (double*)xmalloc(n_obs*sizeof(double));
+	ones   = (double*)xmalloc(n_obs*sizeof(double));
+
+	/*First: Create semivariogram matrix for observations*/
+	for(i=0;i<n_obs;i++){
+		for(j=0;j<=i;j++){
+			//Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
+			Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
+			Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
+		}
+	}
+	for(i=0;i<n_obs;i++) ones[i]=1;
+
+	/*Get semivariogram vector associated to this location*/
+	//for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp,y[i]-y_interp);
+	for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
+
+	/*Solve the three linear systems*/
+#if _HAVE_GSL_
+	SolverxGsl(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
+	SolverxGsl(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
+	SolverxGsl(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
+#else
+	_error_("GSL is required");
+#endif
+
+	/*Prepare predictor*/
+	numerator=-1.; denominator=0.;
+	for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
+	for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
+	ratio=numerator/denominator;
+
+	prediction = 0.;
+	error      = - numerator*numerator/denominator;
+	for(i=0;i<n_obs;i++) prediction += (gamma0[i]-ratio)*GinvZ[i];
+	for(i=0;i<n_obs;i++) error += gamma0[i]*GinvG0[i];
+
+	/*clean-up*/
+	*pprediction = prediction;
+	*perror = error;
+	xfree((void**)&x);
+	xfree((void**)&y);
+	xfree((void**)&obs);
+	xfree((void**)&Gamma);
+	xfree((void**)&gamma0);
+	xfree((void**)&ones);
+	xfree((void**)&GinvG0);
+	xfree((void**)&Ginv1);
+	xfree((void**)&GinvZ);
+
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){
Index: /issm/trunk-jpl/src/c/Container/Observations.h
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.h	(revision 12418)
+++ /issm/trunk-jpl/src/c/Container/Observations.h	(revision 12419)
@@ -8,4 +8,5 @@
 class Obsevration;
 class Quadtree;
+class Variogram;
 class Options;
 
@@ -23,4 +24,7 @@
 
 		/*Methods*/
+		void ClosestObservation(double *px,double *py,double *pobs,double x_interp,double y_interp,double radius);
+		void InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram);
+		void ObservationList(double **px,double **py,double **pobs,int* pnobs);
 		void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata);
 		void QuadtreeColoring(double* A,double *x,double *y,int n);