Changeset 12282

Timestamp:

05/24/12 09:31:58 (13 years ago)

Author:

Mathieu Morlighem

Message:

minor: better deal with nobs=0

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (1 diff)
• objects/Kriging/GaussianVariogram.cpp (modified) (1 diff)

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

@@ -110 +110 @@
 
         /*Find all observations that are in range*/
-        this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range);
-
-        if(nobs==0){
-                /*No observation found, double range*/
-                //printf("No observation found within range, doubling range\n");
+        nobs=0;
+        while(nobs==0){
                 xfree((void**)&indices);
-                this->ObservationList(&x,&y,&obs,&nobs,x_interp,y_interp,range*2);
-        }
-        else{
-                //if(nobs>1000) printf("Taking more than 1000 observations\n");
-                /*Allocate vectors*/
-                x   = (double*)xmalloc(nobs*sizeof(double));
-                y   = (double*)xmalloc(nobs*sizeof(double));
-                obs = (double*)xmalloc(nobs*sizeof(double));
-
-                /*Loop over all observations and fill in x, y and obs*/
-                for (i=0;i<nobs;i++){
-                        observation=(Observation*)this->GetObjectByOffset(indices[i]);
-                        observation->WriteXYObs(&x[i],&y[i],&obs[i]);
-                }
+                this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range);
+                if(!nobs){
+                        /*No observation found, double range*/
+                        range=2*range;
+                }
+        }
+        /*Check that we do not have more than 50 observations*/
+        //while(nobs>50){
+        //      range=range/2;
+        //      xfree((void**)&indices);
+        //      this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range);
+        //}
+
+        if(!nobs) _error_("Unexpected error: no observation within current range and more than 50 observations if doubling range");
+
+        /*Allocate vectors*/
+        x   = (double*)xmalloc(nobs*sizeof(double));
+        y   = (double*)xmalloc(nobs*sizeof(double));
+        obs = (double*)xmalloc(nobs*sizeof(double));
+
+        /*Loop over all observations and fill in x, y and obs*/
+        for (i=0;i<nobs;i++){
+                observation=(Observation*)this->GetObjectByOffset(indices[i]);
+                observation->WriteXYObs(&x[i],&y[i],&obs[i]);
         }
 

issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp

@@ -69 +69 @@
         /*return semi-variogram*/
         a     = 1./3.;
-        gamma = (sill-nugget)*(1-exp(-h2/(a*pow(range,2.)))) + nugget;
+        gamma = (sill-nugget)*(1.-exp(-h2/(a*range*range))) + nugget;
         return gamma;
 }