From dd1503cdee0fe9d5942c5ca8873b0778eab633f9 Mon Sep 17 00:00:00 2001
From: Willem Deconinck <willem.deconinck@ecmwf.int>
Date: Wed, 26 Apr 2023 13:51:27 +0200
Subject: [PATCH] Configurable eps in compare_pointxyz via
 ATLAS_COMPAREPOINTXYZ_EPS_FACTOR env var

New default 1e8*eps instead of 1e4*eps.
For previous behaviour:
 ATLAS_COMPAREPOINTXYZ_EPS_FACTOR=1e4

For high resolution we want to have it to 1e4, for low resolution 1e8
The detection mechanism of source points needs to be replaced with something deterministic that
does not rely on comparing points.
---
 ...nservativeSphericalPolygonInterpolation.cc | 52 ++++++++++---------
 1 file changed, 27 insertions(+), 25 deletions(-)

diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc
index 515868f75..0a39b95e4 100644
--- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc
+++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc
@@ -896,30 +896,6 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr
     }
 }
 
-namespace {
-// needed for intersect_polygons only, merely for detecting duplicate points
-struct ComparePointXYZ {
-    bool operator()(const PointXYZ& f, const PointXYZ& s) const {
-        // eps = ConvexSphericalPolygon::EPS which is the threshold when two points are "same"
-        double eps = 1e4 * std::numeric_limits<double>::epsilon();
-        if (f[0] < s[0] - eps) {
-            return true;
-        }
-        else if (std::abs(f[0] - s[0]) < eps) {
-            if (f[1] < s[1] - eps) {
-                return true;
-            }
-            else if (std::abs(f[1] - s[1]) < eps) {
-                if (f[2] < s[2] - eps) {
-                    return true;
-                }
-            }
-        }
-        return false;
-    }
-};
-}  // namespace
-
 void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolygonArray& src_csp,
                                                                    const CSPolygonArray& tgt_csp) {
     ATLAS_TRACE();
@@ -947,7 +923,33 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg
     StopWatch stopwatch_polygon_intersections;
 
     stopwatch_src_already_in.start();
-    std::set<PointXYZ, ComparePointXYZ> src_cent;
+
+
+    // needed for intersect_polygons only, merely for detecting duplicate points
+    // Treshold at which points are considered same
+    double compare_pointxyz_eps = 1.e8 * std::numeric_limits<double>::epsilon();
+    if (::getenv("ATLAS_COMPAREPOINTXYZ_EPS_FACTOR")) {
+        compare_pointxyz_eps = std::atof(::getenv("ATLAS_COMPAREPOINTXYZ_EPS_FACTOR")) * std::numeric_limits<double>::epsilon();
+    }
+
+    auto compare_pointxyz = [eps=compare_pointxyz_eps] (const PointXYZ& f, const PointXYZ& s) -> bool {
+        if (f[0] < s[0] - eps) {
+            return true;
+        }
+        else if (std::abs(f[0] - s[0]) < eps) {
+            if (f[1] < s[1] - eps) {
+                return true;
+            }
+            else if (std::abs(f[1] - s[1]) < eps) {
+                if (f[2] < s[2] - eps) {
+                    return true;
+                }
+            }
+        }
+        return false;
+    };
+
+    std::set<PointXYZ, decltype(compare_pointxyz)> src_cent(compare_pointxyz);
     auto src_already_in = [&](const PointXYZ& halo_cent) {
         if (src_cent.find(halo_cent) == src_cent.end()) {
             atlas_omp_critical{