sanitize eph of jpl_lp

benchmark/leg_sims_flanagan_benchmark.cpp

@@ -54,7 +54,7 @@ void perform_test_speed(unsigned N, unsigned nseg, unsigned pop_size)
     uda.set_ftol_abs(1e-8);
     uda.set_maxeval(1000);
     pagmo::algorithm algo{uda};
-    //algo.set_verbosity(0u);
+    algo.set_verbosity(0u);
 
     // The initial positions
     kep3::udpla::vsop2013 udpla_earth("earth_moon", 1e-2);

pykep/test_trajopt.py

@@ -56,15 +56,16 @@ class gym_cassini1_tests(_ut.TestCase):
     def test_fitness(self):
         import pykep as pk
         udp = pk.trajopt.gym.cassini1
-        # Three random values. Ground truth provided by the old pykep code
+        # Three random values. Ground truth provided by a snapshot of the kep3 code validated agains the old pykep code 
+        # (not identical, just validated ... as a the old pykep had a buggy jpl_lp eph method nhere fixed)
         x = [-554.5189290104555, 103.27184879471751, 335.41655259663474, 80.50258543604521, 862.0950563689543, 2865.018040480413  ]
         f = udp.fitness(x)[0]
-        self.assertTrue(float_rel_error(f,80400.08897584089) < 1e-14)
+        self.assertTrue(float_rel_error(f,80747.26667221037) < 1e-14)
 
         x = [-932.0532394941108, 37.534681289972674, 162.5093144821548, 336.970139545233, 1743.2915882004586, 2527.8785180526465  ]
         f = udp.fitness(x)[0]
-        self.assertTrue(float_rel_error(f,217105.87501757513) < 1e-14)
+        self.assertTrue(float_rel_error(f,216694.84791232392) < 1e-14)
 
         x = [-583.0776694390058, 388.65047998036107, 334.9959782156864, 65.57508619540917, 1520.2982946551908, 2132.7771932619144 ]
         f = udp.fitness(x)[0]
-        self.assertTrue(float_rel_error(f,107218.0849582104) < 1e-14)
+        self.assertTrue(float_rel_error(f,107787.63141728108) < 1e-14)

src/udpla/jpl_lp.cpp

@@ -8,6 +8,7 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <cmath>
+#include <random>
 #include <stdexcept>
 #include <unordered_map>
 
@@ -47,18 +48,13 @@ static const std::array<double, 6> neptune_el = {30.06992276, 0.00859048, 1.7700
 static const std::array<double, 6> neptune_el_dot = {0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664};
 // clang-format on
 
-// The old AU values used in pykep2 are here used to match legacy code.
-double const MU_SUN = 1.32712440018e20;
-double const AU = 149597870691.0;
-
-
 // NOLINTNEXTLINE(cert-err58-cpp)
 const std::unordered_map<std::string, int> mapped_planets
     = {{"mercury", 1}, {"venus", 2},  {"earth", 3},   {"mars", 4}, {"jupiter", 5},
        {"saturn", 6},  {"uranus", 7}, {"neptune", 8}, {"pluto", 9}};
 
 jpl_lp::jpl_lp(std::string name)
-    : m_elements(), m_elements_dot(), m_name(std::move(name)), m_mu_central_body(MU_SUN)
+    : m_elements(), m_elements_dot(), m_name(std::move(name)), m_mu_central_body(kep3::MU_SUN)
 {
 
     boost::algorithm::to_lower(m_name);
@@ -138,11 +134,11 @@ std::array<double, 6> jpl_lp::_f_elements(double mjd2000) const
     // algorithm from https://ssd.jpl.nasa.gov/planets/approx_pos.html accessed
     // 2023.
     std::array<double, 6> elements_updated{}, elements_f{};
-    double dt = (mjd2000 - 0.5) / 36525.; // Number of centuries passed since J2000.0
+    double dt = (mjd2000) / 36525.; // Number of centuries passed since J2000.0
     for (unsigned int i = 0; i < 6; ++i) {
         elements_updated[i] = (m_elements[i] + m_elements_dot[i] * dt);
     }
-    elements_f[0] = elements_updated[0] * AU;
+    elements_f[0] = elements_updated[0] * kep3::AU;
     elements_f[1] = elements_updated[1];
     elements_f[2] = elements_updated[2] * kep3::DEG2RAD;
     elements_f[3] = elements_updated[5] * kep3::DEG2RAD;
@@ -210,7 +206,7 @@ std::string jpl_lp::get_extra_info() const
     auto pos_vel = eph(0.);
 
     std::string retval = fmt::format("\nLow-precision planet elements: \n")
-                         + fmt::format("Semi major axis (AU): {}\n", par[0] / AU)
+                         + fmt::format("Semi major axis (AU): {}\n", par[0] / kep3::AU)
                          + fmt::format("Eccentricity: {}\n", par[1])
                          + fmt::format("Inclination (deg.): {}\n", par[2] * kep3::RAD2DEG)
                          + fmt::format("Big Omega (deg.): {}\n", par[3] * kep3::RAD2DEG)
@@ -225,7 +221,7 @@ std::string jpl_lp::get_extra_info() const
 
 std::ostream &operator<<(std::ostream &os, const kep3::udpla::jpl_lp &udpla)
 {
-    os << udpla.get_extra_info() << std::endl;
+    os << udpla.get_extra_info() << "/n";
     return os;
 }
 

test/udpla_jpl_lp_test.cpp

@@ -40,14 +40,17 @@ TEST_CASE("constructor")
 
 TEST_CASE("eph")
 {
-    // We use 2030-01-01 as a reference epoch for all these tests
+    // We use 2020-01-01 as a reference epoch for all these tests. To compute the ground truth 
+    // we queried JPL Horizon. Since the queries are made at different times the eph used are also
+    // different. As a consequence if you repeat the query today you may get a different result.
+    // Note also the low error tolerance requested as low precision eph are indeed low precision.
     double ref_epoch = kep3::epoch(2458849.5, kep3::epoch::julian_type::JD).mjd2000();
     {
-        // This is Mercury w.r.t. the Sun queried from JPL Horizon at
+        // This is Mercury w.r.t. the Sun queried from JPL Horizon (DE441) at
         // 2020-01-01
         std::array<std::array<double, 3>, 2> pos_vel_0{
-            {{-9.474762662376745E+09, -6.894147965135109E+10, -4.764334842347469E+09},
-             {3.848711305256677E+04, -4.155242103836629E+03, -3.870162659830893E+03}}};
+            {{-1.004313454665499E+10, -6.782852744259514E+10, -4.760875668975301E+09},
+             {3.847265152100766E+04, -4.158689617302953E+03, -3.869763814829368E+03}}};
         // Mercury in jpl_lp mode
         jpl_lp udpla{"mercury"};
         auto [r, v] = udpla.eph(ref_epoch);