From a92ce6c630fa9113cc31196a094e87de36aa1b42 Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Thu, 28 Sep 2023 16:41:18 +0200
Subject: [PATCH] Fix for orbital elements singularity when using vsop2013 at
 low precision.

---
 src/model/vsop2013.cpp | 22 +++++++++++++++++++---
 test/vsop2013.cpp      | 24 ++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 3 deletions(-)

diff --git a/src/model/vsop2013.cpp b/src/model/vsop2013.cpp
index 5ec31564a..491dd42da 100644
--- a/src/model/vsop2013.cpp
+++ b/src/model/vsop2013.cpp
@@ -324,6 +324,12 @@ std::vector<expression> vsop2013_cartesian_impl(std::uint32_t pl_idx, expression
                              qp = sqrt(qp_2);
                          });
 
+    // Check if qp is zero. This indicates a zero inclination orbit.
+    bool zero_inc = false;
+    if (const auto *num_ptr = std::get_if<number>(&qp.value())) {
+        zero_inc = is_zero(*num_ptr);
+    }
+
     // E, e, sqrt(1 - e**2), cos(i), sin(i), cos(Om), sin(Om), sin(E), cos(E)
     expression E, e, sqrt_1me2, ci, si, cOm, sOm, sin_E, cos_E;
     tbb::parallel_invoke(
@@ -337,12 +343,22 @@ std::vector<expression> vsop2013_cartesian_impl(std::uint32_t pl_idx, expression
             ci = 1_dbl - 2_dbl * qp_2;
             si = sqrt(1_dbl - ci * ci);
         },
-        [&]() { cOm = q / qp; }, [&]() { sOm = p / qp; });
+        // NOTE: for zero inclination, set conventionally Om to zero.
+        [&]() { cOm = zero_inc ? 1_dbl : q / qp; }, [&]() { sOm = zero_inc ? 0_dbl : p / qp; });
+
+    // Check for zero eccentricity.
+    bool zero_ecc = false;
+    if (const auto *num_ptr = std::get_if<number>(&e.value())) {
+        zero_ecc = is_zero(*num_ptr);
+    }
 
     // cos(om), sin(om), q1/a, q2/a.
     expression com, som, q1_a, q2_a;
-    tbb::parallel_invoke([&]() { com = (k * cOm + h * sOm) / e; }, [&]() { som = (h * cOm - k * sOm) / e; },
-                         [&]() { q1_a = cos_E - e; }, [&]() { q2_a = sqrt_1me2 * sin_E; });
+    tbb::parallel_invoke(
+        // NOTE: for zero eccentricity, set conventionally om to zero.
+        [&]() { com = zero_ecc ? 1_dbl : (k * cOm + h * sOm) / e; },
+        [&]() { som = zero_ecc ? 0_dbl : (h * cOm - k * sOm) / e; }, [&]() { q1_a = cos_E - e; },
+        [&]() { q2_a = sqrt_1me2 * sin_E; });
 
     // Prepare the entries of the rotation matrix, and a few auxiliary quantities.
     expression R00, R01, R10, R11, R20, R21, v_num, v_den;
diff --git a/test/vsop2013.cpp b/test/vsop2013.cpp
index 7d9fdae9e..8a00ce100 100644
--- a/test/vsop2013.cpp
+++ b/test/vsop2013.cpp
@@ -14,6 +14,7 @@
 #include <vector>
 
 #include <heyoka/expression.hpp>
+#include <heyoka/llvm_state.hpp>
 #include <heyoka/model/vsop2013.hpp>
 #include <heyoka/taylor.hpp>
 
@@ -780,3 +781,26 @@ TEST_CASE("vsop2013 mus")
     REQUIRE(mus[3] == 8.9970116036316091182e-10);
     REQUIRE(mus[9] == 2.1886997654259696800e-12);
 }
+
+// Bug: low-precision ephmeris may have zero inclination,
+// which leads to singularities when converting to cartesian.
+TEST_CASE("vsop2013 low prec zero inc")
+{
+    auto ex = vsop2013_cartesian(1, kw::time = "tm"_var, kw::thresh = 1e-1)[0];
+
+    llvm_state s;
+
+    add_cfunc<double>(s, "f", {ex});
+
+    s.compile();
+
+    auto *cf_ptr
+        = reinterpret_cast<void (*)(double *, const double *, const double *, const double *)>(s.jit_lookup("f"));
+
+    double out = 0, in = 0;
+
+    cf_ptr(&out, &in, nullptr, nullptr);
+
+    REQUIRE(!std::isnan(out));
+    REQUIRE(out != 0);
+}