From 724878b97a2472c971a3e433918d7beb1261ee80 Mon Sep 17 00:00:00 2001
From: Jerome St-Louis <jerome@ecere.com>
Date: Wed, 14 Aug 2024 06:34:42 -0400
Subject: [PATCH] ISEA: Fix bug introduced in
 5f0aedf1e80a1a979649f9a24b809773513b9b36 - Also added tests for
 `+orient=pole` - step 1: double z = fabs(cosZ) < 1E-15 ? 0 : acos(cosZ); was
 wrong - acos(0) is Pi/2, not 0

---
 src/projections/isea.cpp | 30 ++++++++++++++++++++++--------
 test/gie/builtins.gie    | 37 +++++++++++++++++++++++++++++++++++++
 2 files changed, 59 insertions(+), 8 deletions(-)

diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp
index fb5bb696c5..144bedd53c 100644
--- a/src/projections/isea.cpp
+++ b/src/projections/isea.cpp
@@ -342,6 +342,24 @@ struct pj_isea_data {
 #pragma warning(disable : 4702)
 #endif
 
+#define SAFE_ARC_EPSILON 1E-15
+
+static inline double safeArcSin(double t) {
+    return fabs(t) < SAFE_ARC_EPSILON         ? 0
+           : fabs(t - 1.0) < SAFE_ARC_EPSILON ? M_PI / 2
+           : fabs(t + 1.0) < SAFE_ARC_EPSILON ? -M_PI / 2
+                                              : asin(t);
+}
+
+static inline double safeArcCos(double t) {
+    return fabs(t) < SAFE_ARC_EPSILON       ? M_PI / 2
+           : fabs(t + 1) < SAFE_ARC_EPSILON ? M_PI
+           : fabs(t - 1) < SAFE_ARC_EPSILON ? 0
+                                            : acos(t);
+}
+
+#undef SAFE_ARC_EPSILON
+
 /* coord needs to be in radians */
 static int isea_snyder_forward(const struct pj_isea_data *data,
                                const struct GeoPoint *ll, struct isea_pt *out) {
@@ -371,7 +389,8 @@ static int isea_snyder_forward(const struct pj_isea_data *data,
         double sinAz, cosAz;
 
         /* step 1 */
-        double z = fabs(cosZ) < 1E-15 ? 0 : acos(cosZ);
+        double z = safeArcCos(cosZ);
+
         /* not on this triangle */
         if (z > sdc2vos /*g*/ + 0.000005) { /* TODO DBL_EPSILON */
             continue;
@@ -533,9 +552,7 @@ static struct GeoPoint snyder_ctran(const struct GeoPoint &np,
     while (lambdap < -M_PI)
         lambdap += 2 * M_PI;
 
-    result.lat = fabs(sin_phip - 1.0) < 1E-15   ? M_PI / 2
-                 : fabs(sin_phip + 1.0) < 1E-15 ? -M_PI / 2
-                                                : asin(sin_phip);
+    result.lat = safeArcSin(sin_phip);
     result.lon = lambdap;
     return result;
 }
@@ -1242,10 +1259,7 @@ class ISEAPlanarProjection {
                         double cosLat0SinZ = cosLat0 * sinZ;
                         double latSin =
                             sinLat0 * cosZ + cosLat0SinZ * cos(Az_earth);
-                        double lat = fabs(latSin - 1.0) < 1E-15 ? M_PI / 2
-                                     : fabs(latSin + 1.0) < 1E-15
-                                         ? -M_PI / 2
-                                         : asin(latSin);
+                        double lat = safeArcSin(latSin);
                         double lon =
                             facesCenterDodecahedronVertices[c.face].lon +
                             atan2(sin(Az_earth) * cosLat0SinZ,
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index c500c3fb04..639f5eccff 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -2820,6 +2820,43 @@ accept 0 45
 expect -837334.699958428042009 8323409.759132191538811
 roundtrip 1
 
+-------------------------------------------------------------------------------
+operation +proj=isea +R=6371007.18091875 +orient=pole
+-------------------------------------------------------------------------------
+tolerance 0.2 mm
+
+accept -168.75 58.282525588539
+expect -16702163.549901897087693 6386395.630649688653648
+roundtrip 1
+
+accept 11.25 58.282525588539
+expect 619648.646531744743697 6212947.536539182066917
+roundtrip 1
+
+accept -110 54
+expect -13285649.857057726010680 6149501.348902118392289
+roundtrip 1
+
+accept -75 45
+expect -7921366.529368571005762 4728387.055336073972285
+roundtrip 1
+
+accept 2 49
+expect 152616.434999307675753 5152048.791301283054054
+roundtrip 1
+
+accept 0 0
+expect 0 -195097.133640714135254
+roundtrip 1
+
+accept 90 0
+expect 9593072.435467451811 0
+roundtrip 1
+
+accept 0 45
+expect 0 4726854.770339427515864
+roundtrip 1
+
 ===============================================================================
 # Kavrayskiy V
 # 	PCyl., Sph.