From 67c1fe38d523a35d30407795e1c7710721c8a49d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 1 Nov 2023 09:21:55 +0100 Subject: [PATCH 1/4] Initial files. --- CMakeLists.txt | 1 + include/kep3/detail/s11n.hpp | 2 - include/kep3/planets/jpl_lp.hpp | 5 +- include/kep3/planets/vsop2013.hpp | 55 ++++++++ src/planet.cpp | 15 +- src/planets/vsop2013.cpp | 223 ++++++++++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/planet_vsop2013_test.cpp | 35 +++++ 8 files changed, 325 insertions(+), 12 deletions(-) create mode 100644 include/kep3/planets/vsop2013.hpp create mode 100644 src/planets/vsop2013.cpp create mode 100644 test/planet_vsop2013_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 21c2f096..00a3fb05 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -134,6 +134,7 @@ set(kep3_SRC_FILES "${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/planets/vsop2013.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/detail/s11n.hpp b/include/kep3/detail/s11n.hpp index 23454d73..f318edd2 100644 --- a/include/kep3/detail/s11n.hpp +++ b/include/kep3/detail/s11n.hpp @@ -35,8 +35,6 @@ // Include the archives. #include #include -#include -#include namespace kep3::detail { diff --git a/include/kep3/planets/jpl_lp.hpp b/include/kep3/planets/jpl_lp.hpp index 8f17d3b0..eb31c9d9 100644 --- a/include/kep3/planets/jpl_lp.hpp +++ b/include/kep3/planets/jpl_lp.hpp @@ -68,8 +68,7 @@ class kep3_DLL_PUBLIC jpl_lp [[nodiscard]] double get_radius() const; [[nodiscard]] double get_safe_radius() const; [[nodiscard]] std::string get_extra_info() const; - [[nodiscard]] std::array elements(double = 0., - kep3::elements_type = kep3::elements_type::KEP_F) const; + [[nodiscard]] std::array elements(double = 0., kep3::elements_type = kep3::elements_type::KEP_F) const; private: [[nodiscard]] std::array _f_elements(double = 0.) const; @@ -85,4 +84,4 @@ struct fmt::formatter : ostream_formatter { // necessary for serialization TANUKI_S11N_WRAP_EXPORT_KEY(kep3::udpla::jpl_lp, kep3::detail::planet_iface) -#endif // KEP_TOOLBOX_PLANET_JPL_LP_H \ No newline at end of file +#endif // KEP_TOOLBOX_PLANET_JPL_LP_H diff --git a/include/kep3/planets/vsop2013.hpp b/include/kep3/planets/vsop2013.hpp new file mode 100644 index 00000000..8559f29b --- /dev/null +++ b/include/kep3/planets/vsop2013.hpp @@ -0,0 +1,55 @@ +// 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_VSOP2013_H +#define kep3_PLANET_VSOP2013_H + +#include +#include +#include + +#include +#include +#include + +namespace kep3::udpla +{ + +class kep3_DLL_PUBLIC vsop2013 +{ + struct impl; + + std::unique_ptr m_impl; + + friend class boost::serialization::access; + void save(boost::archive::binary_oarchive &, unsigned) const; + void load(boost::archive::binary_iarchive &, unsigned); + BOOST_SERIALIZATION_SPLIT_MEMBER() + +public: + vsop2013(); + explicit vsop2013(std::string, double = 1e-5); + vsop2013(vsop2013 &&) noexcept; + vsop2013(const vsop2013 &); + vsop2013 &operator=(vsop2013 &&) noexcept; + vsop2013 &operator=(const vsop2013 &); + ~vsop2013(); + + // Mandatory UDPLA methods. + [[nodiscard]] std::array, 2> eph(double) const; + + // Optional UDPLA methods + [[nodiscard]] std::string get_name() const; +}; + +} // namespace kep3::udpla + +TANUKI_S11N_WRAP_EXPORT_KEY(kep3::udpla::vsop2013, kep3::detail::planet_iface) + +#endif diff --git a/src/planet.cpp b/src/planet.cpp index da5d87f9..81e046e0 100644 --- a/src/planet.cpp +++ b/src/planet.cpp @@ -7,9 +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/convert_anomalies.hpp" -#include "kep3/core_astro/eq2par2eq.hpp" -#include "kep3/core_astro/ic2par2ic.hpp" #include #include #include @@ -17,6 +14,9 @@ #include #include +#include +#include +#include #include #include @@ -38,8 +38,8 @@ double period_from_energy(const std::array &r, const std::array elements_from_posvel(const std::array, 2> &pos_vel, - double mu, kep3::elements_type el_type) +std::array elements_from_posvel(const std::array, 2> &pos_vel, double mu, + kep3::elements_type el_type) { std::array retval = kep3::ic2par(pos_vel, mu); switch (el_type) { @@ -60,7 +60,7 @@ std::array elements_from_posvel(const std::array, 2> null_udpla::eph(double) std::array pos = {1., 0., 0.}; std::array vel = {0., 1., 0.}; return {pos, vel}; -}; +} + } // namespace kep3::detail namespace kep3::detail diff --git a/src/planets/vsop2013.cpp b/src/planets/vsop2013.cpp new file mode 100644 index 00000000..f93c3f9c --- /dev/null +++ b/src/planets/vsop2013.cpp @@ -0,0 +1,223 @@ +// 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 +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace kep3::udpla +{ + +namespace detail +{ + +namespace +{ + +// name -> planet index map. +// NOLINTNEXTLINE(cert-err58-cpp) +const std::unordered_map pl_idx_map + = {{"mercury", 1u}, {"venus", 2u}, {"earth_moon", 3u}, {"mars", 4u}, {"jupiter", 5u}, + {"saturn", 6u}, {"uranus", 7u}, {"neptune", 8u}, {"pluto", 9u}}; + +// Hasher for the JIT cache. +struct state_dict_hasher { + std::size_t operator()(const std::pair &p) const noexcept + { + std::size_t seed = std::hash{}(p.first); + boost::hash_combine(seed, std::hash{}(p.second)); + + return seed; + } +}; + +// The JIT cache. +// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) +std::unordered_map, heyoka::llvm_state, state_dict_hasher> state_dict; + +// Mutex for thread-safe access to the cache. +// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) +std::mutex state_dict_mutex; + +} // namespace + +} // namespace detail + +struct vsop2013::impl { + using fptr_t = void (*)(double *, const double *, const double *, const double *) noexcept; + + heyoka::llvm_state m_state; + fptr_t eval_f = nullptr; + double m_mu = 0; + std::string m_pl_name; + + impl() = default; + explicit impl(heyoka::llvm_state s, double mu, std::string pl_name) + : m_state(std::move(s)), m_mu(mu), m_pl_name(std::move(pl_name)) + { + eval_f = reinterpret_cast(m_state.jit_lookup("eval_f")); + } + impl(impl &&) noexcept = delete; + impl(const impl &other) : m_state(other.m_state), m_mu(other.m_mu), m_pl_name(other.m_pl_name) + { + eval_f = reinterpret_cast(m_state.jit_lookup("eval_f")); + } + impl &operator=(impl &&) noexcept = delete; + impl &operator=(const impl &) = delete; + ~impl() = default; + + void save(boost::archive::binary_oarchive &ar, unsigned) const + { + ar << m_state; + ar << m_mu; + ar << m_pl_name; + } + void load(boost::archive::binary_iarchive &ar, unsigned) + { + ar >> m_state; + ar >> m_mu; + ar >> m_pl_name; + + eval_f = reinterpret_cast(m_state.jit_lookup("eval_f")); + } + BOOST_SERIALIZATION_SPLIT_MEMBER() +}; + +vsop2013::vsop2013() : vsop2013("mercury") {} + +vsop2013::vsop2013(std::string pl_name, double thresh) +{ + if (!std::isfinite(thresh) || thresh < 0) { + throw std::invalid_argument( + fmt::format("The threshold value must be finite and non-negative, but it is {} instead", thresh)); + } + + boost::algorithm::to_lower(pl_name); + const auto pl_it = detail::pl_idx_map.find(pl_name); + if (pl_it == detail::pl_idx_map.end()) { + throw std::invalid_argument(fmt::format("The planet name '{}' is invalid", pl_name)); + } + + // Fetch the planet index. + const auto pl_index = pl_it->second; + + // Fetch the planet mu. + auto mu = heyoka::model::get_vsop2013_mus()[pl_index]; + // The original mu is in AU**3/day**2, convert it to SI units. + mu = mu * kep3::AU * kep3::AU * kep3::AU / (kep3::DAY2SEC * kep3::DAY2SEC); + + // Lock down for access to the cache. + std::lock_guard lock(detail::state_dict_mutex); + + if (auto it = detail::state_dict.find({pl_index, thresh}); it == detail::state_dict.end()) { + // Cache miss, we need to create a new compiled function. + + // Time variable. + auto tm = heyoka::expression("tm"); + + // Create the expression for the VSOP2013 solution. + // NOTE: we use (tm - 0.5) / 365250 as time coordinate because: + // - VSOP2013 counts time from J2000 and not MJD2000 (they differ by 12 hours), + // - VSOP2013 measures time in thousands of Julian years. + auto v_ex = heyoka::model::vsop2013_cartesian_icrf(pl_index, heyoka::kw::time = (tm - 0.5) / 365250., + heyoka::kw::thresh = thresh); + + // Create the llvm_state and add the compiled function. + heyoka::llvm_state s; + heyoka::add_cfunc(s, "eval_f", v_ex); + s.compile(); + + // Add the compiled state to the cache. + [[maybe_unused]] auto [_, flag] = detail::state_dict.insert({{pl_index, thresh}, s}); + assert(flag); + + // Build the impl. + m_impl = std::make_unique(std::move(s), mu, std::move(pl_name)); + } else { + // Cache hit, build the impl from the cache content. + m_impl = std::make_unique(it->second, mu, std::move(pl_name)); + } +} + +vsop2013::vsop2013(vsop2013 &&) noexcept = default; + +vsop2013::vsop2013(const vsop2013 &other) : m_impl(std::make_unique(*other.m_impl)) {} + +vsop2013 &vsop2013::operator=(vsop2013 &&) noexcept = default; + +vsop2013 &vsop2013::operator=(const vsop2013 &other) +{ + if (this != &other) { + *this = vsop2013(other); + } + + return *this; +} + +void vsop2013::save(boost::archive::binary_oarchive &ar, unsigned) const +{ + ar << m_impl; +} + +void vsop2013::load(boost::archive::binary_iarchive &ar, unsigned) +{ + try { + ar >> m_impl; + // LCOV_EXCL_START + } catch (...) { + *this = vsop2013{}; + throw; + } + // LCOV_EXCL_STOP +} + +vsop2013::~vsop2013() = default; + +std::array, 2> vsop2013::eph(double mjd2000) const +{ + double out[6]{}; + m_impl->eval_f(out, &mjd2000, nullptr, nullptr); + + return { + {{{out[0] * kep3::AU, out[1] * kep3::AU, out[2] * kep3::AU}}, + {{out[3] * kep3::AU / kep3::DAY2SEC, out[4] * kep3::AU / kep3::DAY2SEC, out[5] * kep3::AU / kep3::DAY2SEC}}}}; +} + +std::string vsop2013::get_name() const +{ + return fmt::format("vsop2013_{}", m_impl->m_pl_name); +} + +} // namespace kep3::udpla + +// NOLINTNEXTLINE +TANUKI_S11N_WRAP_EXPORT_IMPLEMENT(kep3::udpla::vsop2013, kep3::detail::planet_iface) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c648c462..678d22f9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -32,6 +32,7 @@ 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(planet_vsop2013_test) ADD_kep3_TESTCASE(stm_test) ADD_kep3_TESTCASE(ic2par2ic_test) ADD_kep3_TESTCASE(ic2eq2ic_test) diff --git a/test/planet_vsop2013_test.cpp b/test/planet_vsop2013_test.cpp new file mode 100644 index 00000000..9275abe3 --- /dev/null +++ b/test/planet_vsop2013_test.cpp @@ -0,0 +1,35 @@ +// 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 +#include + +#include "catch.hpp" +#include "test_helpers.hpp" + +using kep3::udpla::jpl_lp; +using kep3::udpla::vsop2013; + +TEST_CASE("constructor") +{ + vsop2013 pl1("mercury", 1e-5); + fmt::println("{}", pl1.eph(0)); + + jpl_lp pl2("mercury"); + fmt::println("{}", pl2.eph(0)); +} From 685f5f30b492b398047fa22a7fa47a3de8a2937b Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 23 Nov 2023 08:43:07 +0100 Subject: [PATCH 2/4] A first chunk of testing. --- src/planets/vsop2013.cpp | 16 ++-- test/planet_vsop2013_test.cpp | 175 ++++++++++++++++++++++++++++++++-- 2 files changed, 177 insertions(+), 14 deletions(-) diff --git a/src/planets/vsop2013.cpp b/src/planets/vsop2013.cpp index f93c3f9c..02980449 100644 --- a/src/planets/vsop2013.cpp +++ b/src/planets/vsop2013.cpp @@ -79,15 +79,17 @@ struct vsop2013::impl { fptr_t eval_f = nullptr; double m_mu = 0; std::string m_pl_name; + double m_thresh = 0; impl() = default; - explicit impl(heyoka::llvm_state s, double mu, std::string pl_name) - : m_state(std::move(s)), m_mu(mu), m_pl_name(std::move(pl_name)) + explicit impl(heyoka::llvm_state s, double mu, std::string pl_name, double thresh) + : m_state(std::move(s)), m_mu(mu), m_pl_name(std::move(pl_name)), m_thresh(thresh) { eval_f = reinterpret_cast(m_state.jit_lookup("eval_f")); } impl(impl &&) noexcept = delete; - impl(const impl &other) : m_state(other.m_state), m_mu(other.m_mu), m_pl_name(other.m_pl_name) + impl(const impl &other) + : m_state(other.m_state), m_mu(other.m_mu), m_pl_name(other.m_pl_name), m_thresh(other.m_thresh) { eval_f = reinterpret_cast(m_state.jit_lookup("eval_f")); } @@ -100,12 +102,14 @@ struct vsop2013::impl { ar << m_state; ar << m_mu; ar << m_pl_name; + ar << m_thresh; } void load(boost::archive::binary_iarchive &ar, unsigned) { ar >> m_state; ar >> m_mu; ar >> m_pl_name; + ar >> m_thresh; eval_f = reinterpret_cast(m_state.jit_lookup("eval_f")); } @@ -161,10 +165,10 @@ vsop2013::vsop2013(std::string pl_name, double thresh) assert(flag); // Build the impl. - m_impl = std::make_unique(std::move(s), mu, std::move(pl_name)); + m_impl = std::make_unique(std::move(s), mu, std::move(pl_name), thresh); } else { // Cache hit, build the impl from the cache content. - m_impl = std::make_unique(it->second, mu, std::move(pl_name)); + m_impl = std::make_unique(it->second, mu, std::move(pl_name), thresh); } } @@ -214,7 +218,7 @@ std::array, 2> vsop2013::eph(double mjd2000) const std::string vsop2013::get_name() const { - return fmt::format("vsop2013_{}", m_impl->m_pl_name); + return fmt::format("vsop2013 {}, threshold={}", m_impl->m_pl_name, m_impl->m_thresh); } } // namespace kep3::udpla diff --git a/test/planet_vsop2013_test.cpp b/test/planet_vsop2013_test.cpp index 9275abe3..9fa028bf 100644 --- a/test/planet_vsop2013_test.cpp +++ b/test/planet_vsop2013_test.cpp @@ -7,29 +7,188 @@ // 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 #include #include +#include #include #include #include #include "catch.hpp" -#include "test_helpers.hpp" -using kep3::udpla::jpl_lp; +using kep3::planet; using kep3::udpla::vsop2013; -TEST_CASE("constructor") +TEST_CASE("basic") +{ + using Catch::Matchers::Message; + + { + std::ostringstream oss; + + vsop2013 pl1; + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "mercury")); + REQUIRE(boost::contains(oss.str(), "thresh")); + } + + { + std::ostringstream oss; + + vsop2013 pl1("vEnUs"); + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + } + + { + std::ostringstream oss; + + vsop2013 pl1("venus", 0.5); + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + } + + // Run it again to hit the cache. + { + std::ostringstream oss; + + vsop2013 pl1("venus", 0.5); + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + } + + REQUIRE_THROWS_MATCHES(vsop2013("venus", -1), std::invalid_argument, + Message("The threshold value must be finite and non-negative, but it is -1 instead")); + + REQUIRE_THROWS_MATCHES(vsop2013("goofy"), std::invalid_argument, Message("The planet name 'goofy' is invalid")); + + // Move ctor. + { + std::ostringstream oss; + + vsop2013 pl0("venus", 0.5); + auto pl1 = std::move(pl0); + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + } + + // Copy ctor. + { + std::ostringstream oss; + + vsop2013 pl0("venus", 0.5); + auto pl1 = pl0; + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + + oss.str(""); + + oss << planet(pl0); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + } + + // Move assignment. + { + std::ostringstream oss; + + vsop2013 pl0("venus", 0.5), pl1; + pl1 = std::move(pl0); + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + } + + // Copy assignment. + { + std::ostringstream oss; + + vsop2013 pl0("venus", 0.5), pl1; + pl1 = pl0; + oss << planet(pl1); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + + oss.str(""); + + oss << planet(pl0); + + REQUIRE(boost::contains(oss.str(), "venus")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); + } +} + +TEST_CASE("s11n") +{ + std::stringstream ss; + + vsop2013 pl0("jupiter", 0.5); + + { + boost::archive::binary_oarchive oa(ss); + oa << pl0; + } + + pl0 = vsop2013{}; + + { + boost::archive::binary_iarchive ia(ss); + ia >> pl0; + } + + std::ostringstream oss; + oss << planet(pl0); + + REQUIRE(boost::contains(oss.str(), "jupiter")); + REQUIRE(boost::contains(oss.str(), "thresh")); + REQUIRE(boost::contains(oss.str(), "0.5")); +} + +TEST_CASE("eph") { - vsop2013 pl1("mercury", 1e-5); - fmt::println("{}", pl1.eph(0)); + planet p{vsop2013{"venus", 1e-9}}; + + auto eph = p.eph(123.); + + // NOTE: these values have been manually checked against DE440. + REQUIRE(eph[0][0] == Approx(103304986899.7981)); + REQUIRE(eph[0][1] == Approx(32220404104.1199)); + REQUIRE(eph[0][2] == Approx(7957719449.51538)); - jpl_lp pl2("mercury"); - fmt::println("{}", pl2.eph(0)); + REQUIRE(eph[1][0] == Approx(-10696.505905435035)); + REQUIRE(eph[1][1] == Approx(30061.035989651813)); + REQUIRE(eph[1][2] == Approx(14201.00090492195)); } From 9b29912bfa2b3c3e50c4b31423cb24bff4f22d72 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 23 Nov 2023 08:52:10 +0100 Subject: [PATCH 3/4] Initial Python exposition. --- pykep/docstrings.cpp | 19 ++++++++++++++++++- pykep/docstrings.hpp | 3 +-- pykep/expose_udplas.cpp | 11 +++++++++-- pykep/python_udpla.cpp | 2 ++ 4 files changed, 30 insertions(+), 5 deletions(-) diff --git a/pykep/docstrings.cpp b/pykep/docstrings.cpp index 7901f92c..a8df819a 100644 --- a/pykep/docstrings.cpp +++ b/pykep/docstrings.cpp @@ -1103,6 +1103,24 @@ model from JPL (https://ssd.jpl.nasa.gov/planets/approx_pos.html). )"; } +std::string udpla_vsop2013_docstring() +{ + return R"(__init__(body = "mercury", thresh = 1e-5) + +Constructs a solar system planet with ephemerides computed using the VSOP2013 analytical +theory (https://en.wikipedia.org/wiki/VSOP_model). + +Args: + *body* (:class:`str`): the name of the solar system planet. + *thresh* (:class:`float`): the truncation threshold for the theory's coefficients. + +Examples: + >>> import pykep as pk + >>> udpla = pk.udpla.vsop2013(body="venus") + >>> pla = pk.planet(udpla) +)"; +} + std::string lambert_problem_docstring() { return R"(__init__(r0 = [1,0,0], r1 = [0,1,0], tof = pi/2, mu = 1., cw = False, max_revs = 0) @@ -1206,5 +1224,4 @@ std::string propagate_lagrangian_v_docstring() )"; } - } // namespace pykep \ No newline at end of file diff --git a/pykep/docstrings.hpp b/pykep/docstrings.hpp index c748c28d..9d28b1e2 100644 --- a/pykep/docstrings.hpp +++ b/pykep/docstrings.hpp @@ -76,7 +76,7 @@ std::string planet_elements_docstring(); std::string udpla_keplerian_from_elem_docstring(); std::string udpla_keplerian_from_posvel_docstring(); std::string udpla_jpl_lp_docstring(); - +std::string udpla_vsop2013_docstring(); // Lambert Problem std::string lambert_problem_docstring(); @@ -85,7 +85,6 @@ std::string lambert_problem_docstring(); std::string propagate_lagrangian_docstring(); std::string propagate_lagrangian_v_docstring(); - } // namespace pykep #endif \ No newline at end of file diff --git a/pykep/expose_udplas.cpp b/pykep/expose_udplas.cpp index d0b8c5a0..1d5c130d 100644 --- a/pykep/expose_udplas.cpp +++ b/pykep/expose_udplas.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include "common_utils.hpp" #include "docstrings.hpp" @@ -57,10 +58,16 @@ void expose_all_udplas(py::module &udpla_module, py::class_ &plane udpla_module, planet_class, "_jpl_lp", "Low precision ephemerides from JPL for the solar system planets"); // Constructors. jpl_lp_udpla - .def(py::init(), py::arg("body") = "earth", - pykep::udpla_jpl_lp_docstring().c_str()) + .def(py::init(), py::arg("body") = "earth", pykep::udpla_jpl_lp_docstring().c_str()) // repr(). .def("__repr__", &pykep::ostream_repr); + + // vsop2013 udpla. + auto vsop2013_udpla = pykep::expose_one_udpla( + udpla_module, planet_class, "_vsop2013", "Analytical planetary ephemerides from the VSOP2013 theory"); + // Constructors. + vsop2013_udpla.def(py::init(), py::arg("body") = "mercury", py::arg("thresh") = 1e-5, + pykep::udpla_vsop2013_docstring().c_str()); } } // namespace pykep diff --git a/pykep/python_udpla.cpp b/pykep/python_udpla.cpp index 739c7715..a1e05789 100644 --- a/pykep/python_udpla.cpp +++ b/pykep/python_udpla.cpp @@ -13,8 +13,10 @@ #include #include + #include #include + #include #include From cd7f9d4574ec1c58e3108511af956de38463d4d4 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 23 Nov 2023 09:15:41 +0100 Subject: [PATCH 4/4] Some small testing for the vsop2013 udpla. --- pykep/__init__.py | 41 +++++++++++++++++++++++++++++++++++++++-- pykep/test.py | 12 ++++++++++++ 2 files changed, 51 insertions(+), 2 deletions(-) diff --git a/pykep/__init__.py b/pykep/__init__.py index 179ca323..d5f2ca36 100644 --- a/pykep/__init__.py +++ b/pykep/__init__.py @@ -10,9 +10,43 @@ from ._version import __version__ del _version +import os as _os -# Importing cpp functionalities -from .core import * +if _os.name == "posix": + # NOTE: on some platforms Python by default opens extensions + # with the RTLD_LOCAL flag, which creates problems because + # public symbols used by heyoka (e.g., sleef functions, quad + # precision math) are then not found by the LLVM jit machinery. + # Thus, before importing core, we temporarily flip on the + # RTLD_GLOBAL flag, which makes the symbols visible and + # solves these issues. Another possible approach suggested + # in the llvm discord is to manually and explicitly add + # libheyoka.so to the DL search path: + # DynamicLibrarySearchGenerator::Load(“/path/to/libheyoka.so”) + # See: + # https://docs.python.org/3/library/ctypes.html + import ctypes as _ctypes + import sys as _sys + + _orig_dlopen_flags = _sys.getdlopenflags() + _sys.setdlopenflags(_orig_dlopen_flags | _ctypes.RTLD_GLOBAL) + + try: + # Importing cpp functionalities + from .core import * + finally: + # Restore the original dlopen flags whatever + # happens. + _sys.setdlopenflags(_orig_dlopen_flags) + + del _ctypes + del _sys + del _orig_dlopen_flags +else: + # Importing cpp functionalities + from .core import * + +del _os # Importing python udplas from . import udpla @@ -28,6 +62,9 @@ udpla.jpl_lp = core._jpl_lp udpla.jpl_lp.__name__ = "jpl_lp" udpla.jpl_lp.__module__ = "udpla" +udpla.vsop2013 = core._vsop2013 +udpla.vsop2013.__name__ = "vsop2013" +udpla.vsop2013.__module__ = "udpla" # Importing the python utils from .utils import * diff --git a/pykep/test.py b/pykep/test.py index 5ca0327b..41257b57 100644 --- a/pykep/test.py +++ b/pykep/test.py @@ -258,7 +258,18 @@ def test_spice(self): self.assertTrue(np.allclose(rvpk, np.array(rv)*1000, atol=1e-13)) +class vsop2013_test(_ut.TestCase): + def test_basic(self): + import pykep as pk + + p = pk.planet(pk.udpla.vsop2013()) + self.assertTrue("mercury" in str(p)) + self.assertTrue("1e-05" in str(p)) + + def run_test_suite(): + tl = _ut.TestLoader() + suite = _ut.TestSuite() suite.addTest(anomaly_conversions_tests("test_m2e")) suite.addTest(anomaly_conversions_tests("test_m2f")) @@ -275,6 +286,7 @@ def run_test_suite(): suite.addTest(epoch_test("test_epoch_operators")) suite.addTest(py_udplas_test("test_tle")) suite.addTest(py_udplas_test("test_spice")) + suite.addTest(tl.loadTestsFromTestCase(vsop2013_test))