diff --git a/benchmark/leg_sims_flanagan_benchmark.cpp b/benchmark/leg_sims_flanagan_benchmark.cpp index d2d26f95..eaa1cc10 100644 --- a/benchmark/leg_sims_flanagan_benchmark.cpp +++ b/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); diff --git a/pykep/test_trajopt.py b/pykep/test_trajopt.py index 3dcc413d..69d6a151 100644 --- a/pykep/test_trajopt.py +++ b/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) \ No newline at end of file + self.assertTrue(float_rel_error(f,107787.63141728108) < 1e-14) \ No newline at end of file diff --git a/src/udpla/jpl_lp.cpp b/src/udpla/jpl_lp.cpp index 52ebffe4..e89e92d7 100644 --- a/src/udpla/jpl_lp.cpp +++ b/src/udpla/jpl_lp.cpp @@ -8,6 +8,7 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #include +#include #include #include @@ -47,18 +48,13 @@ static const std::array neptune_el = {30.06992276, 0.00859048, 1.7700 static const std::array 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 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 jpl_lp::_f_elements(double mjd2000) const // algorithm from https://ssd.jpl.nasa.gov/planets/approx_pos.html accessed // 2023. std::array 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; } diff --git a/test/udpla_jpl_lp_test.cpp b/test/udpla_jpl_lp_test.cpp index fe7b6862..8b2a9a22 100644 --- a/test/udpla_jpl_lp_test.cpp +++ b/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, 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);