Merge pull request #29 from bluescarni/pr/vsop2013_planet

VSOP2013 planet

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"

include/kep3/detail/s11n.hpp

@@ -35,8 +35,6 @@
 // Include the archives.
 #include <boost/archive/binary_iarchive.hpp>
 #include <boost/archive/binary_oarchive.hpp>
-#include <boost/archive/text_iarchive.hpp>
-#include <boost/archive/text_oarchive.hpp>
 
 namespace kep3::detail
 {

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<double, 6> elements(double = 0.,
-                                                 kep3::elements_type = kep3::elements_type::KEP_F) const;
+    [[nodiscard]] std::array<double, 6> elements(double = 0., kep3::elements_type = kep3::elements_type::KEP_F) const;
 
 private:
     [[nodiscard]] std::array<double, 6> _f_elements(double = 0.) const;
@@ -85,4 +84,4 @@ struct fmt::formatter<kep3::udpla::jpl_lp> : 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
+#endif // KEP_TOOLBOX_PLANET_JPL_LP_H

include/kep3/planets/vsop2013.hpp

@@ -0,0 +1,55 @@
+// Copyright 2023, 2024 Dario Izzo ([email protected]), Francesco Biscani
+// ([email protected])
+//
+// 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 <array>
+#include <memory>
+#include <string>
+
+#include <kep3/detail/s11n.hpp>
+#include <kep3/detail/visibility.hpp>
+#include <kep3/planet.hpp>
+
+namespace kep3::udpla
+{
+
+class kep3_DLL_PUBLIC vsop2013
+{
+    struct impl;
+
+    std::unique_ptr<impl> 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<std::array<double, 3>, 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

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 *

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

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

pykep/expose_udplas.cpp

@@ -14,6 +14,7 @@
 #include <kep3/planet.hpp>
 #include <kep3/planets/jpl_lp.hpp>
 #include <kep3/planets/keplerian.hpp>
+#include <kep3/planets/vsop2013.hpp>
 
 #include "common_utils.hpp"
 #include "docstrings.hpp"
@@ -57,10 +58,16 @@ void expose_all_udplas(py::module &udpla_module, py::class_<kep3::planet> &plane
         udpla_module, planet_class, "_jpl_lp", "Low precision ephemerides from JPL for the solar system planets");
     // Constructors.
     jpl_lp_udpla
-        .def(py::init<std::string>(), py::arg("body") = "earth",
-             pykep::udpla_jpl_lp_docstring().c_str())
+        .def(py::init<std::string>(), py::arg("body") = "earth", pykep::udpla_jpl_lp_docstring().c_str())
         // repr().
         .def("__repr__", &pykep::ostream_repr<kep3::udpla::jpl_lp>);
+
+    // vsop2013 udpla.
+    auto vsop2013_udpla = pykep::expose_one_udpla<kep3::udpla::vsop2013>(
+        udpla_module, planet_class, "_vsop2013", "Analytical planetary ephemerides from the VSOP2013 theory");
+    // Constructors.
+    vsop2013_udpla.def(py::init<std::string, double>(), py::arg("body") = "mercury", py::arg("thresh") = 1e-5,
+                       pykep::udpla_vsop2013_docstring().c_str());
 }
 
 } // namespace pykep

pykep/python_udpla.cpp

@@ -13,8 +13,10 @@
 
 #include <fmt/core.h>
 #include <fmt/std.h>
+
 #include <kep3/core_astro/constants.hpp>
 #include <kep3/planet.hpp>
+
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 

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))
 
 
 

src/planet.cpp

@@ -7,16 +7,16 @@
 // 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 <array>
 #include <cmath>
 #include <limits>
 
 #include <boost/core/demangle.hpp>
 
 #include <kep3/core_astro/constants.hpp>
+#include <kep3/core_astro/convert_anomalies.hpp>
+#include <kep3/core_astro/eq2par2eq.hpp>
+#include <kep3/core_astro/ic2par2ic.hpp>
 #include <kep3/exceptions.hpp>
 #include <kep3/planet.hpp>
 
@@ -38,8 +38,8 @@ double period_from_energy(const std::array<double, 3> &r, const std::array<doubl
     }
 }
 
-std::array<double, 6> elements_from_posvel(const std::array<std::array<double, 3>, 2> &pos_vel,
-                                                           double mu, kep3::elements_type el_type)
+std::array<double, 6> elements_from_posvel(const std::array<std::array<double, 3>, 2> &pos_vel, double mu,
+                                           kep3::elements_type el_type)
 {
     std::array<double, 6> retval = kep3::ic2par(pos_vel, mu);
     switch (el_type) {
@@ -60,7 +60,7 @@ std::array<double, 6> elements_from_posvel(const std::array<std::array<double, 3
         // LCOV_EXCL_START
         default:
             throw std::logic_error("You should not go here!");
-        // LCOV_EXCL_END
+            // LCOV_EXCL_END
     }
     return retval;
 }
@@ -70,7 +70,8 @@ std::array<std::array<double, 3>, 2> null_udpla::eph(double)
     std::array<double, 3> pos = {1., 0., 0.};
     std::array<double, 3> vel = {0., 1., 0.};
     return {pos, vel};
-};
+}
+
 } // namespace kep3::detail
 
 namespace kep3::detail