refactor and expose ta

CMakeLists.txt

@@ -145,7 +145,7 @@ set(kep3_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/eq2par2eq.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/stm.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/propagate_lagrangian.cpp"
-    "${CMAKE_CURRENT_SOURCE_DIR}/src/ta/lt_kepler.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/ta/stark.cpp"
 )
 
 # Setup of the kep3 shared library.

include/kep3/ta/lt_kepler.hpp → include/kep3/ta/stark.hpp

@@ -21,16 +21,16 @@
 namespace kep3::ta
 {
 // Returns the low-thrust dynamics (heyoka API) in a Keplerian context and Cartesian throttles. 7 states, 5 parameters: mu, veff, ux, uy, uz.
-std::vector<std::pair<heyoka::expression, heyoka::expression>> lt_kepler_dyn();
+kep3_DLL_PUBLIC std::vector<std::pair<heyoka::expression, heyoka::expression>> stark_dyn();
 
 // These return const references to function level static variables of type heyoka::taylor_adaptive<double>.
 // NOTE: The object retruned are expected to be copied to then be modified.
-kep3_DLL_PUBLIC const heyoka::taylor_adaptive<double> &get_ta_lt_kepler(double tol);
-kep3_DLL_PUBLIC const heyoka::taylor_adaptive<double> &get_ta_lt_kepler_var(double tol); // variational (x,y,z,vx,vy,vz,ux,uy,uz) first order
+kep3_DLL_PUBLIC const heyoka::taylor_adaptive<double> &get_ta_stark(double tol);
+kep3_DLL_PUBLIC const heyoka::taylor_adaptive<double> &get_ta_stark_var(double tol); // variational (x,y,z,vx,vy,vz,ux,uy,uz) first order
 
 // Methods to access the cache dimensions.
-kep3_DLL_PUBLIC size_t get_ta_lt_kepler_cache_dim();
-kep3_DLL_PUBLIC size_t get_ta_lt_kepler_var_cache_dim();
+kep3_DLL_PUBLIC size_t get_ta_stark_cache_dim();
+kep3_DLL_PUBLIC size_t get_ta_stark_var_cache_dim();
 } // namespace kep3::ta
 
 #endif

pykep/CMakeLists.txt

@@ -53,6 +53,7 @@ install(TARGETS core
 install(FILES ${PYKEP_PYTHON_FILES} DESTINATION ${_PYKEP_INSTALL_DIR})
 ADD_SUBDIRECTORY(udpla)
 ADD_SUBDIRECTORY(leg)
+ADD_SUBDIRECTORY(ta)
 ADD_SUBDIRECTORY(trajopt)
 ADD_SUBDIRECTORY(utils)
 ADD_SUBDIRECTORY(data)

pykep/__init__.py

@@ -55,6 +55,9 @@
 # Importing trajectory legs
 from . import leg
 
+# Importing Taylor adaptive integrators
+from . import ta
+
 # Importing the python utils
 from .utils import *
 

pykep/core.cpp

@@ -22,6 +22,7 @@
 #include <kep3/leg/sims_flanagan.hpp>
 #include <kep3/planet.hpp>
 #include <kep3/stark_problem.hpp>
+#include <kep3/ta/stark.hpp>
 #include <kep3/udpla/keplerian.hpp>
 
 #include <pybind11/chrono.h>
@@ -288,6 +289,11 @@ PYBIND11_MODULE(core, m)
         .def_property_readonly("iters", &kep3::lambert_problem::get_iters, "The number of iterations made.")
         .def_property_readonly("Nmax", &kep3::lambert_problem::get_Nmax, "The maximum number of iterations allowed.");
 
+    // Exposing taylor adaptive propagators
+    m.def("_stark", &kep3::ta::get_ta_stark, py::arg("tol") = 1e-16, pykep::ta_stark_docstring().c_str());
+    m.def("_stark_var", &kep3::ta::get_ta_stark_var, py::arg("tol") = 1e-16, pykep::ta_stark_var_docstring().c_str());
+    m.def("_stark_dyn", &kep3::ta::stark_dyn, pykep::ta_stark_dyn_docstring().c_str());
+
     // Exposing propagators
     m.def(
         "propagate_lagrangian",
@@ -349,8 +355,8 @@ PYBIND11_MODULE(core, m)
                double tof) {
                 auto sp_retval = sp.propagate_var(rvm_state, thrust, tof);
                 // Lets transfer ownership of dxdx to python (not sure this is actually needed to
-                // get an efficient return value ... maybe its overkill here). It surely avoid one more copy / allocation of 49+21
-                // values, but in the overall algorithm maybe irrelevant. 
+                // get an efficient return value ... maybe its overkill here). It surely avoid one more copy /
+                // allocation of 49+21 values, but in the overall algorithm maybe irrelevant.
                 std::array<double, 49> &dxdx = std::get<1>(sp_retval);
 
                 // We create a capsule for the py::array_t to manage ownership change.

pykep/docstrings.cpp

@@ -1271,6 +1271,14 @@ as well as the gradients of the final states with respect to the thrust directio
 )";
 }
 
+std::string ta_stark_docstring()
+{return "";}
+std::string ta_stark_var_docstring()
+{return "";}
+std::string ta_stark_dyn_docstring()
+{return "";}
+
+
 std::string propagate_lagrangian_docstring()
 {
     return R"(propagate_lagrangian(rv = [[1,0,0], [0,1,0]], tof = pi/2, mu = 1, stm = False)

pykep/docstrings.hpp

@@ -78,6 +78,12 @@ std::string udpla_keplerian_from_posvel_docstring();
 std::string udpla_jpl_lp_docstring();
 std::string udpla_vsop2013_docstring();
 
+// Taylor Adaptive propagators
+std::string ta_stark_docstring();
+std::string ta_stark_var_docstring();
+std::string ta_stark_dyn_docstring();
+
+
 // Lambert Problem
 std::string lambert_problem_docstring();
 

pykep/ta/CMakeLists.txt

@@ -0,0 +1,2 @@
+set(PYKEP_TA_PYTHON_FILES __init__.py)
+install(FILES ${PYKEP_TA_PYTHON_FILES} DESTINATION ${_PYKEP_INSTALL_DIR}/ta)

pykep/ta/__init__.py

@@ -0,0 +1,27 @@
+## 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/.
+
+"""
+Various types of interplanetary trajectory legs
+"""
+from .. import core as _core
+
+# We make sure to register the various conversions from/to heyoka objects
+import heyoka as _hy 
+
+# Renaming cpp taylor adaptive integrators (we need to create an alias first and then 
+# to fool sphinx into thinking these are not aliases, else the sphinx built docs
+# would report them as aliases and fail to document these classes)
+stark = _core._stark
+stark.__module__ = "ta"
+stark_var = _core._stark_var
+stark_var.__module__ = "ta"
+stark_dyn = _core._stark_dyn
+stark_dyn.__module__ = "ta"
+# Removing core from the list of imported symbols.
+del _core
+del _hy

src/stark_problem.cpp

@@ -14,18 +14,18 @@
 
 #include <kep3/core_astro/constants.hpp>
 #include <kep3/stark_problem.hpp>
-#include <kep3/ta/lt_kepler.hpp>
+#include <kep3/ta/stark.hpp>
 
 namespace kep3
 {
 
 using heyoka::taylor_outcome;
-using ta::get_ta_lt_kepler;
-using ta::get_ta_lt_kepler_var;
+using ta::get_ta_stark;
+using ta::get_ta_stark_var;
 
 // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
 stark_problem::stark_problem(double mu, double veff, double tol)
-    : m_mu(mu), m_veff(veff), m_tol(tol), m_ta(get_ta_lt_kepler(tol)), m_ta_var(get_ta_lt_kepler_var(tol)), m_var_ic(70)
+    : m_mu(mu), m_veff(veff), m_tol(tol), m_ta(get_ta_stark(tol)), m_ta_var(get_ta_stark_var(tol)), m_var_ic(70)
 {
     // We set mu and veff for the non variational
     *(m_ta.get_pars_data()) = m_mu;

src/ta/lt_kepler.cpp → src/ta/stark.cpp

@@ -24,7 +24,7 @@
 #include <heyoka/taylor.hpp>
 
 #include <kep3/core_astro/constants.hpp>
-#include <kep3/ta/lt_kepler.hpp>
+#include <kep3/ta/stark.hpp>
 
 using heyoka::eq;
 using heyoka::expression;
@@ -40,7 +40,7 @@ using heyoka::var_ode_sys;
 
 namespace kep3::ta
 {
-std::vector<std::pair<expression, expression>> lt_kepler_dyn()
+std::vector<std::pair<expression, expression>> stark_dyn()
 {
     // The symbolic variables
     auto [x, y, z, vx, vy, vz, m] = make_vars("x", "y", "z", "vx", "vy", "vz", "m");
@@ -71,63 +71,63 @@ std::vector<std::pair<expression, expression>> lt_kepler_dyn()
 };
 
 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
-std::mutex ta_lt_kepler_mutex;
+std::mutex ta_stark_mutex;
 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
-std::unordered_map<double, taylor_adaptive<double>> ta_lt_kepler_cache;
+std::unordered_map<double, taylor_adaptive<double>> ta_stark_cache;
 
-const heyoka::taylor_adaptive<double> &get_ta_lt_kepler(double tol)
+const heyoka::taylor_adaptive<double> &get_ta_stark(double tol)
 {
     // Lock down for access to cache.
-    std::lock_guard const lock(ta_lt_kepler_mutex);
+    std::lock_guard const lock(ta_stark_mutex);
 
     // Lookup.
-    if (auto it = ta_lt_kepler_cache.find(tol); it == ta_lt_kepler_cache.end()) {
+    if (auto it = ta_stark_cache.find(tol); it == ta_stark_cache.end()) {
         // Cache miss, create new one.
         const std::vector init_state = {1., 1., 1., 1., 1., 1., 1.};
-        auto new_ta = taylor_adaptive<double>{lt_kepler_dyn(), init_state, heyoka::kw::tol = tol};
-        return ta_lt_kepler_cache.insert(std::make_pair(tol, std::move(new_ta))).first->second;
+        auto new_ta = taylor_adaptive<double>{stark_dyn(), init_state, heyoka::kw::tol = tol};
+        return ta_stark_cache.insert(std::make_pair(tol, std::move(new_ta))).first->second;
     } else {
         // Cache hit, return existing.
         return it->second;
     }
 }
 
 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
-std::mutex ta_lt_kepler_var_mutex;
+std::mutex ta_stark_var_mutex;
 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
-std::unordered_map<double, taylor_adaptive<double>> ta_lt_kepler_var_cache;
+std::unordered_map<double, taylor_adaptive<double>> ta_stark_var_cache;
 
-const heyoka::taylor_adaptive<double> &get_ta_lt_kepler_var(double tol)
+const heyoka::taylor_adaptive<double> &get_ta_stark_var(double tol)
 {
     // Lock down for access to cache.
-    std::lock_guard const lock(ta_lt_kepler_var_mutex);
+    std::lock_guard const lock(ta_stark_var_mutex);
 
     // Lookup.
-    if (auto it = ta_lt_kepler_var_cache.find(tol); it == ta_lt_kepler_var_cache.end()) {
+    if (auto it = ta_stark_var_cache.find(tol); it == ta_stark_var_cache.end()) {
         auto [x, y, z, vx, vy, vz, m] = make_vars("x", "y", "z", "vx", "vy", "vz", "m");
-        auto vsys = var_ode_sys(lt_kepler_dyn(), {x, y, z, vx, vy, vz, m, par[2], par[3], par[4]}, 1);
+        auto vsys = var_ode_sys(stark_dyn(), {x, y, z, vx, vy, vz, m, par[2], par[3], par[4]}, 1);
         // Cache miss, create new one.
         const std::vector init_state = {1., 1., 1., 1., 1., 1., 1.};
         auto new_ta = taylor_adaptive<double>{vsys, init_state, heyoka::kw::tol = tol, heyoka::kw::compact_mode = true};
-        return ta_lt_kepler_var_cache.insert(std::make_pair(tol, std::move(new_ta))).first->second;
+        return ta_stark_var_cache.insert(std::make_pair(tol, std::move(new_ta))).first->second;
     } else {
         // Cache hit, return existing.
         return it->second;
     }
 }
 
-size_t get_ta_lt_kepler_cache_dim()
+size_t get_ta_stark_cache_dim()
 {
     // Lock down for access to cache.
-    std::lock_guard const lock(ta_lt_kepler_mutex);
-    return ta_lt_kepler_cache.size();
+    std::lock_guard const lock(ta_stark_mutex);
+    return ta_stark_cache.size();
 }
 
-size_t get_ta_lt_kepler_var_cache_dim()
+size_t get_ta_stark_var_cache_dim()
 {
     // Lock down for access to cache.
-    std::lock_guard const lock(ta_lt_kepler_mutex);
-    return ta_lt_kepler_var_cache.size();
+    std::lock_guard const lock(ta_stark_mutex);
+    return ta_stark_var_cache.size();
 }
 
 } // namespace kep3::ta

test/CMakeLists.txt

@@ -41,4 +41,4 @@ ADD_kep3_TESTCASE(propagate_keplerian_test)
 ADD_kep3_TESTCASE(lambert_problem_test)
 ADD_kep3_TESTCASE(stark_problem_test)
 ADD_kep3_TESTCASE(leg_sims_flanagan_test)
-ADD_kep3_TESTCASE(lt_kepler_test)
+ADD_kep3_TESTCASE(ta_stark_test)