uses new max_search_distance of nnn_search in interpolation method spmap

src/core/interp_method_spmap.c

@@ -237,31 +237,6 @@ static size_t check_tgt_result_points(
   return new_offset;
 }
 
-static void check_max_search_distance(
-  struct yac_interp_grid * interp_grid, double max_search_distance,
-  size_t * num_src_points_, size_t * src_points,
-  yac_coordinate_pointer src_coords, size_t * tgt_result_points) {
-
-  if (max_search_distance > 0.0) {
-    yac_const_coordinate_pointer tgt_field_coords =
-      yac_interp_grid_get_tgt_field_coords(interp_grid);
-    size_t const num_src_points = *num_src_points_;
-    size_t new_num_src_points = 0;
-    for (size_t i = 0; i < num_src_points; ++i) {
-      if (get_vector_angle(
-            src_coords[i], tgt_field_coords[tgt_result_points[i]]) <=
-          max_search_distance) {
-        if (i != new_num_src_points) {
-          src_points[new_num_src_points] = src_points[i];
-          tgt_result_points[new_num_src_points] = tgt_result_points[i];
-        }
-        new_num_src_points++;
-      }
-    }
-    *num_src_points_ = new_num_src_points;
-  }
-}
-
 static double * compute_weights(
   struct yac_interp_grid * interp_grid,
   enum yac_interp_spmap_weight_type weight_type, size_t num_src_points,
@@ -522,21 +497,27 @@ static size_t do_search_spmap (struct interp_method * method,
   size_t * tgt_result_points =
     xmalloc(num_src_points * sizeof(*tgt_result_points));
   yac_interp_grid_do_nnn_search_tgt(
-    interp_grid, src_coords, num_src_points, 1, tgt_result_points, M_PI);
-
-  // check that we found a target for each source point
-  for (size_t i = 0; i < num_src_points; ++i)
-    YAC_ASSERT(
-      tgt_result_points[i] != SIZE_MAX,
-      "ERROR(do_search_spmap): could not find matching target point")
-
-  // remove results that exceed the maximum search distance
-  check_max_search_distance(
-    interp_grid, ((struct interp_method_spmap*)method)->max_search_distance,
-    &num_src_points, src_points, src_coords, tgt_result_points);
+    interp_grid, src_coords, num_src_points, 1, tgt_result_points,
+    ((struct interp_method_spmap*)method)->max_search_distance);
 
   free(src_coords);
 
+  // remove source points for which matching target point was found
+  {
+    size_t new_num_src_points = 0;
+    for (size_t i = 0; i < num_src_points; ++i) {
+      if (tgt_result_points[i] != SIZE_MAX) {
+        if (i != new_num_src_points) {
+          src_points[new_num_src_points] = src_points[i];
+          tgt_result_points[new_num_src_points] =
+            tgt_result_points[i];
+        }
+        ++new_num_src_points;
+      }
+    }
+    num_src_points = new_num_src_points;
+  }
+
   // spread the data from each source point to multiple target points
   // (weight_data is set to NULL if no spreading was applied)
   double * weight_data;
@@ -638,7 +619,8 @@ struct interp_method * yac_interp_method_spmap_new(
 
   method->vtable = &interp_method_spmap_vtable;
   method->spread_distance = spread_distance;
-  method->max_search_distance = max_search_distance;
+  method->max_search_distance =
+    (max_search_distance == 0.0)?M_PI:max_search_distance;
   method->weight_type = weight_type;
   method->scale_type = scale_type;
   method->src_area_scale = src_sphere_radius * src_sphere_radius;