Kepler's equation precision fix

benchmark/convert_anomalies_benchmark.cpp

@@ -34,7 +34,7 @@ void perform_test_speed(double min_ecc, double max_ecc, unsigned N) {
   // Distributions
   //
   std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
-  std::uniform_real_distribution<double> M_d(-100, 100.);
+  std::uniform_real_distribution<double> M_d(-1e8, 1e8);
 
   // We generate the random dataset
   std::vector<double> eccenricities(N);
@@ -68,7 +68,7 @@ void perform_test_accuracy(double min_ecc, double max_ecc, unsigned N) {
   // Distributions
   //
   std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
-  std::uniform_real_distribution<double> M_d(-100, 100.);
+  std::uniform_real_distribution<double> M_d(-1e8, 1e8);
 
   // We generate the random dataset
   std::vector<double> eccenricities(N);
@@ -85,7 +85,10 @@ void perform_test_accuracy(double min_ecc, double max_ecc, unsigned N) {
   std::vector<double> err(N);
   for (auto i = 0u; i < N; ++i) {
     auto res = e2m(m2e(mean_anomalies[i], eccenricities[i]), eccenricities[i]);
-    err[i] = std::abs(res - mean_anomalies[i]);
+    // error is arbitrarily: (|sinM-sinMtrue| +|cosM-cosMtrue|)/2
+    err[i] = (std::abs(std::sin(res) - std::sin(mean_anomalies[i])) +
+              std::abs(std::cos(res) - std::cos(mean_anomalies[i]))) /
+             2.;
   }
   auto max_it = max_element(std::begin(err), std::end(err));
   auto min_it = min_element(std::begin(err), std::end(err));

include/kep3/core_astro/convert_anomalies.hpp

@@ -10,6 +10,7 @@
 #ifndef kep3_CONVERT_ANOMALIES_H
 #define kep3_CONVERT_ANOMALIES_H
 
+#include "kep3/core_astro/constants.hpp"
 #include <cmath>
 
 #include <boost/math/constants/constants.hpp>
@@ -21,19 +22,32 @@ namespace kep3 {
 
 // mean to eccentric (only ellipses) e<1. Preserves the sign and integer number
 // of revolutions.
+// NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
 inline double m2e(double M, double ecc) {
-  // We use as intial guess the Mean Anomaly
-  // (tests indicated that any higher order expansion does not really improve)
-  double IG = M;
+  // We compute the sin and cos of the mean anomaly which are used also later
+  // for the initial guess (3rd order expansion of Kepler's equation).
+  double sinM = std::sin(M), cosM = std::cos(M);
+  // Here we use the atan2 to recover the mean anomaly in the [0,2pi] range.
+  // This makes sure that for high value of M no catastrophic cancellation
+  // occurs, as would be the case using std::fmod(M, 2pi)
+  double M_cropped = std::atan2(sinM, cosM);
+  if (M_cropped < 0) {
+    M_cropped += 2 * kep3::pi;
+  }
+  // The Initial guess follows from a third order expansion of Kepler's equation.
+  double IG = M_cropped + ecc * sinM + ecc * ecc * sinM * cosM +
+              ecc * ecc * ecc * sinM * (1.5 * cosM * cosM - 0.5);
+
   const int digits = std::numeric_limits<double>::digits;
   std::uintmax_t max_iter = 100u;
-  double sol = boost::math::tools::halley_iterate(
-      [M, ecc](double E) {
-        return std::make_tuple(kepE(E, M, ecc), d_kepE(E, ecc),
-                               dd_kepE(E, ecc));
+
+  // Newton-raphson or Halley iterates can be used here. Similar performances, thus
+  // we choose the simplest algorithm.
+  double sol = boost::math::tools::newton_raphson_iterate(
+      [M_cropped, ecc](double E) {
+        return std::make_tuple(kepE(E, M_cropped, ecc), d_kepE(E, ecc));
       },
-      IG, IG - boost::math::constants::pi<double>(),
-      IG + boost::math::constants::pi<double>(), digits, max_iter);
+      IG, IG - kep3::pi, IG + kep3::pi, digits, max_iter);
   if (max_iter == 100u) {
     throw std::domain_error(
         "Maximum number of iterations exceeded when solving Kepler's "
@@ -54,13 +68,9 @@ inline double f2e(double f, double ecc) {
 }
 
 // mean to true (only ellipses) e<1 (returns in range [-pi,pi])
-inline double m2f(double M, double ecc) {
-  return e2f(m2e(M, ecc), ecc);
-}
+inline double m2f(double M, double ecc) { return e2f(m2e(M, ecc), ecc); }
 // true to mean (only ellipses) e<1 (returns in range [-pi,pi])
-inline double f2m(double f, double ecc) {
-  return e2m(f2e(f, ecc), ecc);
-}
+inline double f2m(double f, double ecc) { return e2m(f2e(f, ecc), ecc); }
 
 // gudermannian to true (only hyperbolas) e>1 (returns in range [-pi,pi])
 inline double zeta2f(double f, double ecc) {

include/kep3/planets/keplerian.hpp

@@ -47,11 +47,12 @@ class kep3_DLL_PUBLIC keplerian {
   }
 
 public:
-  // NOTE: in here elem is a,e,i,W,w,M (Mean anomaly, not true anomaly)
+  // NOTE: in here par is a,e,i,W,w,M (Mean anomaly, not true anomaly)
   // NOTE: added_param contains mu_self, radius and safe_radius
   explicit keplerian(const epoch &ref_epoch, const std::array<double, 6> &par,
                      double mu_central_body = 1., std::string name = "Unknown",
                      std::array<double, 3> added_params = {-1., -1., -1.});
+  // Constructor from pos_vel
   explicit keplerian(const epoch &ref_epoch = kep3::epoch(0),
                      const std::array<std::array<double, 3>, 2> &pos_vel =
                          {{{1.0, 0.0, 0.0}, {0., 1.0, 0.0}}},

test/convert_anomalies_test.cpp

@@ -14,13 +14,12 @@
 
 #include "catch.hpp"
 
-using kep3::e2m;
-using kep3::m2e;
 using kep3::e2f;
+using kep3::e2m;
 using kep3::f2e;
-using kep3::zeta2f;
 using kep3::f2zeta;
-
+using kep3::m2e;
+using kep3::zeta2f;
 
 TEST_CASE("m2e") {
   using Catch::Detail::Approx;
@@ -34,31 +33,30 @@ TEST_CASE("m2e") {
   //
   std::uniform_real_distribution<double> ecc_difficult_d(0.9, 0.99);
   std::uniform_real_distribution<double> ecc_easy_d(0., 0.9);
-  std::uniform_real_distribution<double> M_d(-100, 100.);
+  std::uniform_real_distribution<double> M_d(-1e6, 1e6);
 
   // Testing on N random calls (easy)
   unsigned N = 10000;
   for (auto i = 0u; i < N; ++i) {
     double mean_anom = M_d(rng_engine);
     double ecc = ecc_easy_d(rng_engine);
     double res = e2m(m2e(mean_anom, ecc), ecc);
-    // We require 1e-13 and not 1e-14 as the range for m is [-100,100]. For high
-    // values of M some precision loss is to be expected as to maintain digits
-    // of the higher number.
+
     REQUIRE(std::sin(res) ==
-            Approx(std::sin(mean_anom)).epsilon(0.).margin(1e-13));
+            Approx(std::sin(mean_anom)).epsilon(0.).margin(1e-14));
     REQUIRE(std::cos(res) ==
-            Approx(std::cos(mean_anom)).epsilon(0.).margin(1e-13));
+            Approx(std::cos(mean_anom)).epsilon(0.).margin(1e-14));
   }
-  // Testing on N random calls (difficult)
+
+  // Testing on N random calls (difficult e> 0.9 e<1)
   for (auto i = 0u; i < N; ++i) {
     double mean_anom = M_d(rng_engine);
     double ecc = ecc_difficult_d(rng_engine);
     double res = e2m(m2e(mean_anom, ecc), ecc);
     REQUIRE(std::sin(res) ==
-            Approx(std::sin(mean_anom)).epsilon(0.).margin(1e-13));
+            Approx(std::sin(mean_anom)).epsilon(0.).margin(1e-14));
     REQUIRE(std::cos(res) ==
-            Approx(std::cos(mean_anom)).epsilon(0.).margin(1e-13));
+            Approx(std::cos(mean_anom)).epsilon(0.).margin(1e-14));
   }
 }
 

test/planet_keplerian_test.cpp

@@ -8,12 +8,57 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <fmt/core.h>
+#include <fmt/ranges.h>
 
+#include <kep3/core_astro/constants.hpp>
 #include <kep3/planets/keplerian.hpp>
 
 #include "catch.hpp"
 
 TEST_CASE("constructor") {
-  kep3::udpla::keplerian udpla{};
-  fmt::print("{}", udpla);
+  REQUIRE_NOTHROW(kep3::udpla::keplerian{});
+  kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000};
+  // From posvel
+  REQUIRE_NOTHROW(kep3::udpla::keplerian{
+      ref_epoch, {{{0.3, 1., 0.2}, {0.0, 1.12, 0.}}}, 1., "unknown"});
+  REQUIRE_NOTHROW(kep3::udpla::keplerian{ref_epoch,
+                                         {{{0.3, 1., 0.2}, {0.0, 1.12, 0.}}},
+                                         1.,
+                                         "unknown",
+                                         {-1, -1, -1}});
+  // From parameters
+  std::array<double, 6> par{{1., 0., 0., 0., 0., 0.}};
+  REQUIRE_NOTHROW(kep3::udpla::keplerian{ref_epoch, par, 1., "unknown"});
+  REQUIRE_NOTHROW(
+      kep3::udpla::keplerian{ref_epoch, par, 1., "unknown", {-1, -1, -1}});
+  // Checking the data members initializations:
+  kep3::udpla::keplerian udpla{ref_epoch, par, 1.1, "unknown", {1.2, 2.2, 1.9}};
+  REQUIRE(udpla.get_ref_epoch() == ref_epoch);
+  REQUIRE(udpla.get_name() == "unknown");
+  REQUIRE(udpla.get_mu_central_body() == 1.1);
+  REQUIRE(udpla.get_mu_self() == 1.2);
+  REQUIRE(udpla.get_radius() == 2.2);
+  REQUIRE(udpla.get_safe_radius() == 1.9);
+}
+
+TEST_CASE("eph") {
+  kep3::epoch ref_epoch{0., kep3::epoch::MJD2000};
+  kep3::udpla::keplerian udpla1{
+      ref_epoch,
+      {{{kep3::AU, 0., 0.}, {0., kep3::EARTH_VELOCITY, 0.}}},
+      kep3::MU_SUN};
+  kep3::udpla::keplerian udpla2{
+      ref_epoch,
+      {{{kep3::AU, 0., 0.}, {0., -kep3::EARTH_VELOCITY, 0.}}},
+      kep3::MU_SUN};
+  std::array<double, 6> par3{{kep3::AU, 0., 0., 0., 0., 0.}};
+  kep3::udpla::keplerian udpla3{ref_epoch, par3};
+  std::array<double, 6> par4{{kep3::AU, 0., 0., 0., 0., 0.}};
+  kep3::udpla::keplerian udpla4{ref_epoch, par4};
+  double period_in_days =
+      (2 * kep3::pi * kep3::AU / kep3::EARTH_VELOCITY) * kep3::SEC2DAY;
+  auto [r, v] = udpla1.eph(ref_epoch + 200* period_in_days);
+  fmt::print("r: {},{},{}", r[0] / kep3::AU, r[1] / kep3::AU, r[2] / kep3::AU);
+  fmt::print("r: {},{},{}", v[0] / kep3::EARTH_VELOCITY,
+             v[1] / kep3::EARTH_VELOCITY, v[2] / kep3::EARTH_VELOCITY);
 }