From fb51afa3c5a553a61bf495b85bea3d57a7ddeeb2 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sun, 24 Sep 2023 22:27:42 +0200 Subject: [PATCH 01/10] constructor added --- include/kep3/planets/jpl_lp.hpp | 83 +++++++++++++++++ include/kep3/planets/keplerian.hpp | 5 +- src/planets/jpl_lp.cpp | 141 +++++++++++++++++++++++++++++ src/planets/keplerian.cpp | 2 +- 4 files changed, 228 insertions(+), 3 deletions(-) create mode 100644 include/kep3/planets/jpl_lp.hpp create mode 100644 src/planets/jpl_lp.cpp diff --git a/include/kep3/planets/jpl_lp.hpp b/include/kep3/planets/jpl_lp.hpp new file mode 100644 index 00000000..7b9ab57d --- /dev/null +++ b/include/kep3/planets/jpl_lp.hpp @@ -0,0 +1,83 @@ +// Copyright 2023, 2024 Dario Izzo (dario.izzo@gmail.com), Francesco Biscani +// (bluescarni@gmail.com) +// +// This file is part of the kep3 library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef kep3_PLANET_JPL_LP_H +#define kep3_PLANET_JPL_LP_H + +#include + +#include + +#include +#include +#include +#include +#include + +namespace kep3::udpla { + +/// Solar System Planet (jpl simplified ephemerides) +/** + * This class allows to instantiate planets of + * the solar system by referring to their common names. The ephemeris used + * are low_precision ephemeris defined in the JPL pages: + * https://ssd.jpl.nasa.gov/planets/approx_pos.html and valid for the timeframe 1800AD - + * 2050 AD + */ + +class kep3_DLL_PUBLIC jpl_lp { + + std::array m_elements; + std::array m_elements_dot; + std::string m_name; + double m_mu_central_body; + double m_mu_self; + double m_radius; + double m_safe_radius; + + friend class boost::serialization::access; + template void serialize(Archive &ar, unsigned) { + ar &m_elements; + ar &m_elements_dot; + ar &m_name; + ar &m_mu_central_body; + ar &m_mu_self; + ar &m_radius; + ar &m_safe_radius; + } + +public: + // Constructor + explicit jpl_lp(const std::string & = "earth"); + // Mandatory UDPLA methods + [[nodiscard]] std::array, 2> eph(const epoch &) const; + + // Optional UDPLA methods + [[nodiscard]] std::string get_name() const; + [[nodiscard]] double get_mu_central_body() const; + [[nodiscard]] double get_mu_self() const; + [[nodiscard]] double get_radius() const; + [[nodiscard]] double get_safe_radius() const; + [[nodiscard]] std::string get_extra_info() const; + + // Other methods + [[nodiscard]] std::array + elements(const kep3::epoch & = kep3::epoch(), + kep3::elements_type = kep3::elements_type::KEP_F) const; +}; +kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &, + const kep3::udpla::jpl_lp &); +} // namespace kep3::udpla + +// fmt formatter redirecting to the stream operator +template <> struct fmt::formatter : ostream_formatter {}; +// necessary for serialization +kep3_S11N_PLANET_EXPORT_KEY(kep3::udpla::jpl_lp); + +#endif // KEP_TOOLBOX_PLANET_JPL_LP_H \ No newline at end of file diff --git a/include/kep3/planets/keplerian.hpp b/include/kep3/planets/keplerian.hpp index 7c9998be..66cf120c 100644 --- a/include/kep3/planets/keplerian.hpp +++ b/include/kep3/planets/keplerian.hpp @@ -1,4 +1,3 @@ - // Copyright 2023, 2024 Dario Izzo (dario.izzo@gmail.com), Francesco Biscani // (bluescarni@gmail.com) // @@ -80,8 +79,10 @@ kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &, const kep3::udpla::keplerian &); } // namespace kep3::udpla +// fmt formatter redirecting to the stream operator template <> struct fmt::formatter : ostream_formatter {}; +// necessary for serialization kep3_S11N_PLANET_EXPORT_KEY(kep3::udpla::keplerian); -#endif // kep3_EPOCH_H \ No newline at end of file +#endif // kep3_UDPLA_KEPLERIAN_H \ No newline at end of file diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp new file mode 100644 index 00000000..a5d830b6 --- /dev/null +++ b/src/planets/jpl_lp.cpp @@ -0,0 +1,141 @@ +// Copyright 2023, 2024 Dario Izzo (dario.izzo@gmail.com), Francesco Biscani +// (bluescarni@gmail.com) +// +// This file is part of the kep3 library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +namespace kep3::udpla { + +// Data from: https://ssd.jpl.nasa.gov/planets/approx_pos.html +// clang-format off +static const std::array mercury_el = {0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593}; +static const std::array mercury_el_dot = {0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081}; +static const std::array venus_el = {0.72333566, 0.00677672, 3.39467605, 181.97909950, 131.60246718, 76.67984255}; +static const std::array venus_el_dot = {0.00000390, -0.00004107, -0.00078890, 58517.81538729, 0.00268329, -0.27769418}; +static const std::array earth_moon_el = {1.00000261, 0.01671123, -0.00001531, 100.46457166, 102.93768193, 0.0}; +static const std::array earth_moon_el_dot = {0.00000562, -0.00004392, -0.01294668, 35999.37244981, 0.32327364, 0.0}; +static const std::array mars_el = {1.52371034, 0.09339410, 1.84969142, -4.55343205, -23.94362959, 49.55953891}; +static const std::array mars_el_dot = {0.00001847, 0.00007882, -0.00813131, 19140.30268499, 0.44441088, -0.29257343}; +static const std::array jupiter_el = {5.20288700, 0.04838624, 1.30439695, 34.39644051, 14.72847983, 100.47390909}; +static const std::array jupiter_el_dot = {-0.00011607, -0.00013253, -0.00183714, 3034.74612775, 0.21252668, 0.20469106}; +static const std::array saturn_el = {9.53667594, 0.05386179, 2.48599187, 49.95424423, 92.59887831, 113.66242448}; +static const std::array saturn_el_dot = {-0.00125060, -0.00050991, 0.00193609, 1222.49362201, -0.41897216, -0.28867794}; +static const std::array uranus_el = {19.18916464, 0.04725744, 0.77263783, 313.23810451, 170.95427630, 74.01692503}; +static const std::array uranus_el_dot = {-0.00196176, -0.00004397, -0.00242939, 428.48202785, 0.40805281, 0.04240589}; +static const std::array neptune_el = {30.06992276, 0.00859048, 1.77004347, -55.12002969, 44.96476227, 131.78422574}; +static const std::array neptune_el_dot = {0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664}; +// clang-format on + +// Cannot be marked const as operator[] is not and we use it. +static 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(const std::string &name) + : m_elements(), m_elements_dot(), m_name(name), + m_mu_central_body(kep3::MU_SUN) { + + std::string lower_case_name = name; + boost::algorithm::to_lower(lower_case_name); + switch (mapped_planets[lower_case_name]) { + case (1): { + m_elements = mercury_el; + m_elements_dot = mercury_el_dot; + m_radius = 2440000.; + m_safe_radius = 1.1; + m_mu_self = 22032e9; + } break; + case (2): { + m_elements = venus_el; + m_elements_dot = venus_el_dot; + m_radius = 6052000.; + m_safe_radius = 1.1; + m_mu_self = 324859e9; + } break; + case (3): { + m_elements = earth_moon_el; + m_elements_dot = earth_moon_el_dot; + m_radius = 6378000.; + m_safe_radius = 1.1; + m_mu_self = 398600.4418e9; + } break; + case (4): { + m_elements = mars_el; + m_elements_dot = mars_el_dot; + m_radius = 3397000.; + m_safe_radius = 1.1; + m_mu_self = 42828e9; + } break; + case (5): { + m_elements = jupiter_el; + m_elements_dot = jupiter_el_dot; + m_radius = 71492000.; + m_safe_radius = 9.; + m_mu_self = 126686534e9; + } break; + case (6): { + m_elements = saturn_el; + m_elements_dot = saturn_el_dot; + m_radius = 60330000.; + m_safe_radius = 1.1; + m_mu_self = 37931187e9; + } break; + case (7): { + m_elements = uranus_el; + m_elements_dot = uranus_el_dot; + m_radius = 25362000.; + m_safe_radius = 1.1; + m_mu_self = 5793939e9; + } break; + case (8): { + m_elements = neptune_el; + m_elements_dot = neptune_el_dot; + m_radius = 24622000.; + m_safe_radius = 1.1; + m_mu_self = 6836529e9; + } break; + default: { + throw std::logic_error("unknown planet name: "); + } + } +} + +std::string jpl_lp::get_name() const { return m_name; } + +double jpl_lp::get_mu_central_body() const { return m_mu_central_body; } + +double jpl_lp::get_mu_self() const { return m_mu_self; } + +double jpl_lp::get_radius() const { return m_radius; } + +double jpl_lp::get_safe_radius() const { return m_safe_radius; } + +std::string jpl_lp::get_extra_info() const { + auto par = elements(); + std::string retval = fmt::format("Keplerian planet elements: \n"); + return retval; +} + +std::ostream &operator<<(std::ostream &os, const kep3::udpla::jpl_lp &udpla) { + os << udpla.get_extra_info() << std::endl; + return os; +} + +} // namespace kep3::udpla + +kep3_S11N_PLANET_IMPLEMENT(kep3::udpla::jpl_lp) diff --git a/src/planets/keplerian.cpp b/src/planets/keplerian.cpp index 02125aa5..f926ae1b 100644 --- a/src/planets/keplerian.cpp +++ b/src/planets/keplerian.cpp @@ -7,7 +7,6 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -#include "kep3/core_astro/eq2par2eq.hpp" #include #include #include @@ -18,6 +17,7 @@ #include #include +#include #include #include #include From 863b610c912a719a954e6b9bc3d6a7fe08ff811f Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sun, 24 Sep 2023 22:30:06 +0200 Subject: [PATCH 02/10] small change --- src/planets/jpl_lp.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index a5d830b6..92dab75b 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -41,15 +41,16 @@ 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 -// Cannot be marked const as operator[] is not and we use it. -static 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(const std::string &name) : m_elements(), m_elements_dot(), m_name(name), m_mu_central_body(kep3::MU_SUN) { + // Cannot be marked const as operator[] is not and we use it. + static std::unordered_map mapped_planets = { + {"mercury", 1}, {"venus", 2}, {"earth", 3}, + {"mars", 4}, {"jupiter", 5}, {"saturn", 6}, + {"uranus", 7}, {"neptune", 8}, {"pluto", 9}}; + std::string lower_case_name = name; boost::algorithm::to_lower(lower_case_name); switch (mapped_planets[lower_case_name]) { From 5d9f618c90324bf52fd8428da4f1138d2ff93c48 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Wed, 27 Sep 2023 19:42:29 +0200 Subject: [PATCH 03/10] jpl low-precision ephs --- CMakeLists.txt | 1 + include/kep3/planets/jpl_lp.hpp | 8 +++- src/planets/jpl_lp.cpp | 74 ++++++++++++++++++++++++++++++++- 3 files changed, 80 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 752df921..384663cc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -133,6 +133,7 @@ set(kep3_SRC_FILES "${CMAKE_CURRENT_SOURCE_DIR}/src/planet.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/lambert_problem.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/planets/keplerian.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/src/planets/jpl_lp.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/ic2par2ic.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/ic2eq2ic.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/eq2par2eq.cpp" diff --git a/include/kep3/planets/jpl_lp.hpp b/include/kep3/planets/jpl_lp.hpp index 7b9ab57d..78e92e38 100644 --- a/include/kep3/planets/jpl_lp.hpp +++ b/include/kep3/planets/jpl_lp.hpp @@ -27,8 +27,8 @@ namespace kep3::udpla { * This class allows to instantiate planets of * the solar system by referring to their common names. The ephemeris used * are low_precision ephemeris defined in the JPL pages: - * https://ssd.jpl.nasa.gov/planets/approx_pos.html and valid for the timeframe 1800AD - - * 2050 AD + * https://ssd.jpl.nasa.gov/planets/approx_pos.html and valid for the timeframe + * 1800AD - 2050 AD */ class kep3_DLL_PUBLIC jpl_lp { @@ -70,6 +70,10 @@ class kep3_DLL_PUBLIC jpl_lp { [[nodiscard]] std::array elements(const kep3::epoch & = kep3::epoch(), kep3::elements_type = kep3::elements_type::KEP_F) const; + +private: + [[nodiscard]] std::array + _f_elements(const kep3::epoch & = kep3::epoch()) const; }; kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &, const kep3::udpla::jpl_lp &); diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index 92dab75b..bf619f8f 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -7,6 +7,8 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "kep3/core_astro/eq2par2eq.hpp" +#include "kep3/epoch.hpp" #include #include #include @@ -16,12 +18,14 @@ #include #include +#include #include #include namespace kep3::udpla { // Data from: https://ssd.jpl.nasa.gov/planets/approx_pos.html +// a,e,i,L,W,w // clang-format off static const std::array mercury_el = {0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593}; static const std::array mercury_el_dot = {0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081}; @@ -116,6 +120,57 @@ jpl_lp::jpl_lp(const std::string &name) } } +// Computes the kep3::KEP_F elements (osculating with true anomaly) at epoch. +std::array jpl_lp::_f_elements(const kep3::epoch &ep) const { + double mjd2000 = ep.mjd2000(); + if (mjd2000 <= -73048.0 || mjd2000 >= 18263.0) { + throw std::domain_error("Low precision Ephemeris are only valid in the " + "range range [1800-2050]"); + } + // algorithm from https://ssd.jpl.nasa.gov/planets/approx_pos.html accessed + // 2023. + std::array elements_updated{}, elements_f{}; + double dt = (mjd2000 - 2451545.0) / + 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] * get_mu_central_body(); + elements_f[1] = elements_updated[1]; + elements_f[2] = elements_updated[2] * kep3::DEG2RAD; + elements_f[3] = elements_updated[5] * kep3::DEG2RAD; + elements_f[4] = (elements_updated[4] - elements_updated[5]) * kep3::DEG2RAD; + elements_f[5] = (elements_updated[3] - elements_updated[4]) * kep3::DEG2RAD; + elements_f[5] = kep3::m2e(elements_updated[5], elements_updated[1]); + return elements_f; +} + +std::array, 2> jpl_lp::eph(const epoch &ep) const { + auto elements_f = _f_elements(ep); + return par2ic(elements_f, get_mu_central_body()); +} + +std::array jpl_lp::elements(const kep3::epoch &ep, + kep3::elements_type el_type) const { + auto elements = _f_elements(ep); + switch (el_type) { + case kep3::elements_type::KEP_F: + break; + case kep3::elements_type::KEP_M: + elements[5] = kep3::f2m(elements[5], elements[1]); + break; + case kep3::elements_type::MEQ: + elements = kep3::par2eq(elements, false); + break; + case kep3::elements_type::MEQ_R: + elements = kep3::par2eq(elements, true); + break; + default: + throw std::logic_error("You should not go here!"); + } + return elements; +} + std::string jpl_lp::get_name() const { return m_name; } double jpl_lp::get_mu_central_body() const { return m_mu_central_body; } @@ -128,7 +183,24 @@ double jpl_lp::get_safe_radius() const { return m_safe_radius; } std::string jpl_lp::get_extra_info() const { auto par = elements(); - std::string retval = fmt::format("Keplerian planet elements: \n"); + kep3::epoch ep{0., kep3::epoch::MJD2000}; + auto pos_vel = eph(ep); + + std::string retval = + fmt::format("Low-precision planet elements: \n") + + 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) + + fmt::format("Small omega (deg.): {}\n", par[4] * kep3::RAD2DEG) + + fmt::format("True anomly (deg.): {}\n", par[5] * kep3::RAD2DEG); + retval += fmt::format("Mean anomly (deg.): {}\n", + kep3::f2m(par[5], par[1]) * kep3::RAD2DEG); + retval += + fmt::format("Elements reference epoch (MJD2000): {}\n", ep.mjd2000()) + + fmt::format("Elements reference epoch (date): {}\n", ep) + + fmt::format("r at ref. = {}\n", pos_vel[0]) + + fmt::format("v at ref. = {}\n", pos_vel[1]); return retval; } From 9e1c718a248f5bfe424c3a49c1f783d0f44e26ef Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Wed, 27 Sep 2023 22:22:45 +0200 Subject: [PATCH 04/10] test added wip --- test/planet_jpl_lp_test.cpp | 66 +++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 test/planet_jpl_lp_test.cpp diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp new file mode 100644 index 00000000..9fa4efb5 --- /dev/null +++ b/test/planet_jpl_lp_test.cpp @@ -0,0 +1,66 @@ +// Copyright 2023, 2024 Dario Izzo (dario.izzo@gmail.com), Francesco Biscani +// (bluescarni@gmail.com) +// +// This file is part of the kep3 library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "catch.hpp" +#include "test_helpers.hpp" + +using kep3::udpla::jpl_lp; + +TEST_CASE("constructor") { + REQUIRE_NOTHROW(jpl_lp{}); + kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000}; + // From name + REQUIRE_NOTHROW(jpl_lp{"Mars"}); + REQUIRE_NOTHROW(jpl_lp{"mars"}); + REQUIRE_NOTHROW(kep3::planet{jpl_lp{"neptune"}}); + + REQUIRE_THROWS_AS(jpl_lp{"gigi"}, std::logic_error); + jpl_lp udpla{"earth"}; + kep3::planet earth{udpla}; + REQUIRE(kep3_tests::floating_point_error(earth.period() * kep3::SEC2DAY, 365.25) < 0.01); +} + +TEST_CASE("stream_operator") { + REQUIRE_NOTHROW((std::cout << jpl_lp{} << '\n')); +} + +TEST_CASE("serialization_test") { + // Instantiate a generic jpl_lp + kep3::epoch ref_epoch{2423.4343, kep3::epoch::MJD2000}; + jpl_lp udpla{"neptune"}; + + // Store the string representation. + std::stringstream ss; + auto before = boost::lexical_cast(udpla); + // Now serialize + { + boost::archive::binary_oarchive oarchive(ss); + oarchive << udpla; + } + // Deserialize + // Create a new udpla object + jpl_lp udpla2{}; + { + boost::archive::binary_iarchive iarchive(ss); + iarchive >> udpla2; + } + auto after = boost::lexical_cast(udpla2); + // Compare the string represetation + REQUIRE(before == after); +} \ No newline at end of file From 1513ebc8441011ed823fa2ed718a3a62d94eca22 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Fri, 29 Sep 2023 22:23:48 +0200 Subject: [PATCH 05/10] wip --- include/kep3/planet.hpp | 2 +- src/planets/jpl_lp.cpp | 11 +++++------ test/CMakeLists.txt | 1 + 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/include/kep3/planet.hpp b/include/kep3/planet.hpp index 308f7990..3fdf3038 100644 --- a/include/kep3/planet.hpp +++ b/include/kep3/planet.hpp @@ -228,7 +228,7 @@ struct kep3_DLL_PUBLIC_INLINE_CLASS planet_inner final : planet_inner_base { // period from the energy at epoch auto [r, v] = eph(ep); double mu = get_mu_central_body(); - double R = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + double R = std::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]); double v2 = v[0] * v[0] + v[1] * v[1] + v[2] * v[2]; double en = v2 / 2. - mu / R; if (en > 0) { diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index bf619f8f..33b00f03 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -130,18 +130,17 @@ std::array jpl_lp::_f_elements(const kep3::epoch &ep) const { // algorithm from https://ssd.jpl.nasa.gov/planets/approx_pos.html accessed // 2023. std::array elements_updated{}, elements_f{}; - double dt = (mjd2000 - 2451545.0) / - 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] * get_mu_central_body(); + 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; elements_f[4] = (elements_updated[4] - elements_updated[5]) * kep3::DEG2RAD; elements_f[5] = (elements_updated[3] - elements_updated[4]) * kep3::DEG2RAD; - elements_f[5] = kep3::m2e(elements_updated[5], elements_updated[1]); + elements_f[5] = kep3::m2f(elements_f[5], elements_f[1]); return elements_f; } @@ -182,12 +181,12 @@ double jpl_lp::get_radius() const { return m_radius; } double jpl_lp::get_safe_radius() const { return m_safe_radius; } std::string jpl_lp::get_extra_info() const { - auto par = elements(); kep3::epoch ep{0., kep3::epoch::MJD2000}; + auto par = elements(ep); auto pos_vel = eph(ep); std::string retval = - fmt::format("Low-precision planet elements: \n") + + fmt::format("\nLow-precision planet elements: \n") + 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) + diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0bbc24fb..e077a613 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -30,6 +30,7 @@ ADD_kep3_TESTCASE(convert_anomalies_test) ADD_kep3_TESTCASE(epoch_test) ADD_kep3_TESTCASE(planet_test) ADD_kep3_TESTCASE(planet_keplerian_test) +ADD_kep3_TESTCASE(planet_jpl_lp_test) ADD_kep3_TESTCASE(ic2par2ic_test) ADD_kep3_TESTCASE(ic2eq2ic_test) ADD_kep3_TESTCASE(eq2par2eq_test) From 3f9879201ace70f71931db4e95255341c5cfe0ae Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Fri, 29 Sep 2023 23:52:28 +0200 Subject: [PATCH 06/10] test round finished --- test/planet_jpl_lp_test.cpp | 86 ++++++++++++++++++++++++++++++++++++- 1 file changed, 85 insertions(+), 1 deletion(-) diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp index 9fa4efb5..f92811bb 100644 --- a/test/planet_jpl_lp_test.cpp +++ b/test/planet_jpl_lp_test.cpp @@ -33,7 +33,91 @@ TEST_CASE("constructor") { REQUIRE_THROWS_AS(jpl_lp{"gigi"}, std::logic_error); jpl_lp udpla{"earth"}; kep3::planet earth{udpla}; - REQUIRE(kep3_tests::floating_point_error(earth.period() * kep3::SEC2DAY, 365.25) < 0.01); + REQUIRE(kep3_tests::floating_point_error(earth.period() * kep3::SEC2DAY, + 365.25) < 0.01); +} + +TEST_CASE("eph") { + // We use 2030-01-01 as a reference epoch for all these tests + kep3::epoch ref_epoch{2458849.5, kep3::epoch::JD}; + { + // This is the Earth-Moon w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{-2.488023054631234E+10, 1.449771522542222E+11, + -6.590293144971132E+02}, + {-2.984589828430694E+04, -5.151004951052294E+03, + 3.108878527788850E-01}}}; + // The Earth in jpl_lp mode + jpl_lp udpla{"earth"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); + } + { + // This is Mercury w.r.t. the Sun queried from JPL Horizon 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}}}; + // Mercury in jpl_lp mode + jpl_lp udpla{"mercury"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.05); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.05); + } + { + // This is Venus w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{1.081892249749067E+11, 7.861125522626230E+09, -6.135421905733132E+09}, + {-2.679023504807336E+03, 3.476995213020635E+04, + 6.316923820826013E+02}}}; + // Venus in jpl_lp mode + jpl_lp udpla{"venus"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.02); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.02); + } +} + +TEST_CASE("elements") { + kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000}; + // We use Neptune + jpl_lp udpla{"nePTUne"}; // casing is not important + auto pos_vel = udpla.eph(ref_epoch); + // Test on various element types + { + auto par = udpla.elements(ref_epoch, kep3::elements_type::KEP_F); + auto [r, v] = kep3::par2ic(par, kep3::MU_SUN); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel[0]) < 1e-13); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel[1]) < 1e-13); + } + { + auto par = udpla.elements(ref_epoch, kep3::elements_type::KEP_M); + par[5] = kep3::m2f(par[5], par[1]); + auto [r, v] = kep3::par2ic(par, kep3::MU_SUN); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel[0]) < 1e-13); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel[1]) < 1e-13); + } + { + auto par = udpla.elements(ref_epoch, kep3::elements_type::MEQ); + auto [r, v] = kep3::eq2ic(par, kep3::MU_SUN); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel[0]) < 1e-13); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel[1]) < 1e-13); + } + { + auto par = udpla.elements(ref_epoch, kep3::elements_type::MEQ_R); + auto [r, v] = kep3::eq2ic(par, kep3::MU_SUN, true); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel[0]) < 1e-13); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel[1]) < 1e-13); + } + { + REQUIRE_THROWS_AS(udpla.elements(ref_epoch, kep3::elements_type::POSVEL), + std::logic_error); + } } TEST_CASE("stream_operator") { From ec1dd6fe5b67fdb0725c43a7f169658b8f7dcfa1 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sat, 30 Sep 2023 00:02:50 +0200 Subject: [PATCH 07/10] test round finished --- src/planets/jpl_lp.cpp | 16 ++++++++-------- test/planet_jpl_lp_test.cpp | 9 +++++++++ 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index 33b00f03..d4487c55 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -62,56 +62,56 @@ jpl_lp::jpl_lp(const std::string &name) m_elements = mercury_el; m_elements_dot = mercury_el_dot; m_radius = 2440000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 22032e9; } break; case (2): { m_elements = venus_el; m_elements_dot = venus_el_dot; m_radius = 6052000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 324859e9; } break; case (3): { m_elements = earth_moon_el; m_elements_dot = earth_moon_el_dot; m_radius = 6378000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 398600.4418e9; } break; case (4): { m_elements = mars_el; m_elements_dot = mars_el_dot; m_radius = 3397000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 42828e9; } break; case (5): { m_elements = jupiter_el; m_elements_dot = jupiter_el_dot; m_radius = 71492000.; - m_safe_radius = 9.; + m_safe_radius = 9. * m_radius; m_mu_self = 126686534e9; } break; case (6): { m_elements = saturn_el; m_elements_dot = saturn_el_dot; m_radius = 60330000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 37931187e9; } break; case (7): { m_elements = uranus_el; m_elements_dot = uranus_el_dot; m_radius = 25362000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 5793939e9; } break; case (8): { m_elements = neptune_el; m_elements_dot = neptune_el_dot; m_radius = 24622000.; - m_safe_radius = 1.1; + m_safe_radius = 1.1 * m_radius; m_mu_self = 6836529e9; } break; default: { diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp index f92811bb..64a4e3eb 100644 --- a/test/planet_jpl_lp_test.cpp +++ b/test/planet_jpl_lp_test.cpp @@ -120,6 +120,15 @@ TEST_CASE("elements") { } } +TEST_CASE("getters") { + jpl_lp udpla{"nePTUne"}; // casing is not important + REQUIRE(udpla.get_name() == "neptune"); + REQUIRE(udpla.get_mu_central_body() == kep3::MU_SUN); + REQUIRE(udpla.get_mu_self() == 6836529e9); + REQUIRE(udpla.get_radius() == 24622000.); + REQUIRE(udpla.get_safe_radius() == 1.1 * 24622000.); +} + TEST_CASE("stream_operator") { REQUIRE_NOTHROW((std::cout << jpl_lp{} << '\n')); } From 8786c1909bb310851c54757c3775e52e0844e013 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sat, 30 Sep 2023 00:11:10 +0200 Subject: [PATCH 08/10] test coverage improved --- src/planets/jpl_lp.cpp | 5 ++--- test/planet_jpl_lp_test.cpp | 40 +++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 3 deletions(-) diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index d4487c55..93c04f77 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -55,9 +55,8 @@ jpl_lp::jpl_lp(const std::string &name) {"mars", 4}, {"jupiter", 5}, {"saturn", 6}, {"uranus", 7}, {"neptune", 8}, {"pluto", 9}}; - std::string lower_case_name = name; - boost::algorithm::to_lower(lower_case_name); - switch (mapped_planets[lower_case_name]) { + boost::algorithm::to_lower(m_name); + switch (mapped_planets[m_name]) { case (1): { m_elements = mercury_el; m_elements_dot = mercury_el_dot; diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp index 64a4e3eb..e784a978 100644 --- a/test/planet_jpl_lp_test.cpp +++ b/test/planet_jpl_lp_test.cpp @@ -81,6 +81,46 @@ TEST_CASE("eph") { REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.02); REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.02); } + { + // This is Mars w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{-1.974852868472516E+11, -1.325074305699912E+11, + 2.068800235454373E+09}, + {1.440720082303430E+04, -1.804659323991406E+04, + -7.316474757792575E+02}}}; + // Mars in jpl_lp mode + jpl_lp udpla{"mars"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); + } + { + // This is Jupiter w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{7.871048884706007E+10, -7.780620023532844E+11, 1.470618693758428E+09}, + {1.285491315086331E+04, 1.933721291664580E+03, + -2.956500740610059E+02}}}; + // Jupiter in jpl_lp mode + jpl_lp udpla{"jupiter"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); + } + { + // This is Uranus w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{2.427299831783689E+12, 1.702279661193254E+12, -2.511573638659978E+09}, + {-3.947889815690411E+03, 5.259970919483185E+03, + 7.055304073010626E+01}}}; + // Uranus in jpl_lp mode + jpl_lp udpla{"uranus"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); + } } TEST_CASE("elements") { From 206fa703a6e0272a483646e0885f314c11cbb1db Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sat, 30 Sep 2023 00:21:01 +0200 Subject: [PATCH 09/10] full coverage --- test/planet_jpl_lp_test.cpp | 45 +++++++++++++++++++++++++------------ 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp index e784a978..d1bd4e2d 100644 --- a/test/planet_jpl_lp_test.cpp +++ b/test/planet_jpl_lp_test.cpp @@ -40,20 +40,6 @@ TEST_CASE("constructor") { TEST_CASE("eph") { // We use 2030-01-01 as a reference epoch for all these tests kep3::epoch ref_epoch{2458849.5, kep3::epoch::JD}; - { - // This is the Earth-Moon w.r.t. the Sun queried from JPL Horizon at - // 2020-01-01 - std::array, 2> pos_vel_0{ - {{-2.488023054631234E+10, 1.449771522542222E+11, - -6.590293144971132E+02}, - {-2.984589828430694E+04, -5.151004951052294E+03, - 3.108878527788850E-01}}}; - // The Earth in jpl_lp mode - jpl_lp udpla{"earth"}; - auto [r, v] = udpla.eph(ref_epoch); - REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); - REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); - } { // This is Mercury w.r.t. the Sun queried from JPL Horizon at // 2020-01-01 @@ -81,6 +67,21 @@ TEST_CASE("eph") { REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.02); REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.02); } + { + // This is the Earth-Moon w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{-2.488023054631234E+10, 1.449771522542222E+11, + -6.590293144971132E+02}, + {-2.984589828430694E+04, -5.151004951052294E+03, + 3.108878527788850E-01}}}; + // The Earth in jpl_lp mode + jpl_lp udpla{"earth"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); + } + { // This is Mars w.r.t. the Sun queried from JPL Horizon at // 2020-01-01 @@ -108,6 +109,19 @@ TEST_CASE("eph") { REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); } + { + // This is Saturn w.r.t. the Sun queried from JPL Horizon at + // 2020-01-01 + std::array, 2> pos_vel_0{ + {{5.680597453102431E+11, -1.389479460523918E+12, 1.545819892540634E+09}, + {8.420955066542843E+03, 3.631222339233865E+03, + -3.987639953503348E+02}}}; + // Uranus in jpl_lp mode + jpl_lp udpla{"sAtURN"}; + auto [r, v] = udpla.eph(ref_epoch); + REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); + REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); + } { // This is Uranus w.r.t. the Sun queried from JPL Horizon at // 2020-01-01 @@ -121,6 +135,9 @@ TEST_CASE("eph") { REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01); REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); } + jpl_lp udpla{"uranus"}; + REQUIRE_THROWS_AS(udpla.eph(kep3::epoch(5347534, kep3::epoch::MJD2000)), + std::domain_error); } TEST_CASE("elements") { From 20cbb52128ceca39b1eaa04151c22fb1fef027ad Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sat, 30 Sep 2023 00:22:56 +0200 Subject: [PATCH 10/10] full coverage --- test/planet_jpl_lp_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp index d1bd4e2d..fc2b829d 100644 --- a/test/planet_jpl_lp_test.cpp +++ b/test/planet_jpl_lp_test.cpp @@ -75,7 +75,7 @@ TEST_CASE("eph") { -6.590293144971132E+02}, {-2.984589828430694E+04, -5.151004951052294E+03, 3.108878527788850E-01}}}; - // The Earth in jpl_lp mode + // The Earth-Moon in jpl_lp mode jpl_lp udpla{"earth"}; auto [r, v] = udpla.eph(ref_epoch); REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 0.01);