cr3bp added

doc/taylor_adaptive.rst

@@ -22,3 +22,15 @@ Stark dynamics
 
 .. autofunction:: stark_dyn 
 
+
+Circular Restricted Three Body Problem
+=======================================
+
+.. currentmodule:: pykep.ta
+
+.. autofunction:: get_cr3bp
+
+.. autofunction:: get_cr3bp_var
+
+.. autofunction:: cr3bp_dyn 
+

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/cr3bp.hpp>
 #include <kep3/ta/stark.hpp>
 #include <kep3/udpla/keplerian.hpp>
 
@@ -290,9 +291,42 @@ PYBIND11_MODULE(core, m)
         .def_property_readonly("Nmax", &kep3::lambert_problem::get_Nmax, "The maximum number of iterations allowed.");
 
     // Exposing taylor adaptive propagators
-    m.def("_get_stark", &kep3::ta::get_ta_stark, py::arg("tol") = 1e-16, pykep::get_stark_docstring().c_str());
-    m.def("_get_stark_var", &kep3::ta::get_ta_stark_var, py::arg("tol") = 1e-16, pykep::get_stark_var_docstring().c_str());
+    // Stark
+    m.def(
+        "_get_stark",
+        [](double tol) {
+            auto ta_cache = kep3::ta::get_ta_stark(tol);
+            heyoka::taylor_adaptive<double> ta(ta_cache);
+            return ta;
+        },
+        py::arg("tol") = 1e-16, pykep::get_stark_docstring().c_str());
+    m.def(
+        "_get_stark_var",
+        [](double tol) {
+            auto ta_cache = kep3::ta::get_ta_stark_var(tol);
+            heyoka::taylor_adaptive<double> ta(ta_cache);
+            return ta;
+        },
+        py::arg("tol") = 1e-16, pykep::get_stark_var_docstring().c_str());
     m.def("_stark_dyn", &kep3::ta::stark_dyn, pykep::stark_dyn_docstring().c_str());
+    // CR3BP
+    m.def(
+        "_get_cr3bp",
+        [](double tol) {
+            auto ta_cache = kep3::ta::get_ta_cr3bp(tol);
+            heyoka::taylor_adaptive<double> ta(ta_cache);
+            return ta;
+        },
+        py::arg("tol") = 1e-16, pykep::get_cr3bp_docstring().c_str());
+    m.def(
+        "_get_cr3bp_var",
+        [](double tol) {
+            auto ta_cache = kep3::ta::get_ta_cr3bp_var(tol);
+            heyoka::taylor_adaptive<double> ta(ta_cache);
+            return ta;
+        },
+        py::arg("tol") = 1e-16, pykep::get_cr3bp_var_docstring().c_str());
+    m.def("_cr3bp_dyn", &kep3::ta::cr3bp_dyn, pykep::cr3bp_dyn_docstring().c_str());
 
     // Exposing propagators
     m.def(

pykep/docstrings.cpp

@@ -1286,7 +1286,7 @@ std::string get_stark_docstring()
 {
     return R"(get_stark(tol)
 
-Gets the Taylor adaptive propagator (Heyoka) for the Stark problem from the global cache.
+Gets the Taylor adaptive propagator (Heyoka) for the Stark problem from the global cache and returns a copy.
 
 In `pykep`, abusing a term well established in electrodynamics, 
 this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.
@@ -1320,7 +1320,7 @@ std::string get_stark_var_docstring()
 {
     return R"(get_stark_var(tol)
 
-Gets the variational (order 1) Taylor adaptive propagator (Heyoka) for the Stark problem from the global cache.
+Gets the variational (order 1) Taylor adaptive propagator (Heyoka) for the Stark problem from the global cache and returns a copy.
 
 .. note:
    Variations are only considered with repsect to initial conditions and the fixed inertial thurst.
@@ -1373,6 +1373,98 @@ where :math:`\mu, v_{eff} = I_{sp}g_0` and :math:`\mathbf T = [T_x, T_y, T_z]` a
 )";
 }
 
+std::string get_cr3bp_docstring()
+{
+    return R"(get_cr3bp(tol)
+
+Gets the Taylor adaptive propagator (Heyoka) for the CR3BP problem from the global cache and returns a copy. 
+If the requested propagator was never created this will create it, else it will
+return the one from the global cache, thus avoiding jitting.
+
+In `pykep`, the CR3BP is defined in Cartesian coordinates (thus non symplectic as not in a Hamiltonian form). 
+
+The dynamics is that returned by :func:`~pykep.ta.cr3bp_dyn`.
+
+Args:
+    *tol* (:class:`float`): the tolerance of the Taylor adaptive propagator. 
+
+Returns:
+    :class:`hy::taylor_adaptive`: The Taylor adaptive propagator.
+
+Examples:
+  >>> import pykep as pk
+  >>> ta = pk.ta.get_cr3bp(tol = 1e-16)
+  >>> ta.time = 0.
+  >>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
+  >>> mu = 0.01215058560962404
+  >>> tof = 5.7856656782589234
+  >>> ta.pars[0] = mu
+  >>> ta.propagate_until(tof)
+)";
+}
+
+std::string get_cr3bp_var_docstring()
+{
+    return R"(get_cr3bp_var(tol)
+
+Gets the variational (order 1) Taylor adaptive propagator (Heyoka) for the CR3BP problem from the global cache and returns a copy.
+If the requested propagator was never created this will create it, else it will
+return the one from the global cache, thus avoiding jitting.
+
+.. note:
+   Variations are only considered with respect to initial conditions.
+
+In `pykep`, the CR3BP is defined in Cartesian coordinates (thus non symplectic as not in a Hamiltonian form). 
+
+The dynamics is that returned by :func:`~pykep.ta.cr3bp_dyn`: and also used in :func:`~pykep.ta.get_cr3bp`
+
+Args:
+    *tol* (:class:`float`): the tolerance of the Taylor adaptive propagator. 
+
+Returns:
+    :class:`hy::taylor_adaptive`: The Taylor adaptive propagator.
+
+Examples:
+  >>> import pykep as pk
+  >>> ta = pk.ta.get_cr3bp_var(tol = 1e-16)
+  >>> ta.time = 0.
+  >>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
+  >>> mu = 0.01215058560962404
+  >>> tof = 5.7856656782589234
+  >>> ta.pars[0] = mu
+  >>> ta.propagate_until(tof)
+)";
+}
+
+std::string cr3bp_dyn_docstring()
+{return R"(cr3bp_dyn()
+
+The dynamics of the Circular Restricted Body Problem. 
+
+In `pykep`, the CR3BP is defined in Cartesian coordinates (thus non symplectic as not in a Hamiltonian form). 
+
+The parameter :math:`\mu` is defined as :math:`\frac{m_2}{m_1+m_2}` where :math:`m_2` is the mass of the
+secondary body (i.e. placed on the positive x axis). 
+
+The equations are non-dimensional with units :math:`L = r_{12}` (distance between the primaries), :math:`M = m_1 + m_2` (total system mass) and
+:math:`T = \sqrt(\frac{r_{12}^3}{m_1+m_2})` (period of rotation of the primaries).
+
+.. math::
+   \left\{
+   \begin{array}{l}
+       \dot{\mathbf r} = \mathbf v \\
+       \dot v_x = 2v_y + x - (1 - \mu) \frac{x + \mu}{r_1^3} - \mu \frac{x + \mu - 1}{r_2^3} \\
+       \dot v_y = -2 v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3} \\
+       \dot v_z = -(1 - \mu) \frac{z}{r_1^3} - \mu \frac{z}{r_2^3}
+   \end{array}\right.
+
+where :math:`\mu` is the only parameter.
+
+Returns:
+    :class:`list` [ :class:`tuple` (:class:`hy::expression`, :class:`hy::expression` )]: The dynamics in the form [(x, dx), ...]
+)";
+}
+
 
 std::string propagate_lagrangian_docstring()
 {

pykep/docstrings.hpp

@@ -82,7 +82,9 @@ std::string udpla_vsop2013_docstring();
 std::string get_stark_docstring();
 std::string get_stark_var_docstring();
 std::string stark_dyn_docstring();
-
+std::string get_cr3bp_docstring();
+std::string get_cr3bp_var_docstring();
+std::string cr3bp_dyn_docstring();
 
 // Lambert Problem
 std::string lambert_problem_docstring();

pykep/ta/__init__.py

@@ -23,6 +23,13 @@
 get_stark_var.__module__ = "ta"
 stark_dyn = _core._stark_dyn
 stark_dyn.__module__ = "ta"
+
+get_cr3bp = _core._get_cr3bp
+get_cr3bp.__module__ = "ta"
+get_cr3bp_var = _core._get_cr3bp_var
+get_cr3bp_var.__module__ = "ta"
+cr3bp_dyn = _core._cr3bp_dyn
+cr3bp_dyn.__module__ = "ta"
 # Removing core from the list of imported symbols.
 del _core
 del _hy

src/ta/stark.cpp

@@ -8,6 +8,7 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <array>
+#include <heyoka/kw.hpp>
 #include <mutex>
 #include <tuple>
 #include <unordered_map>
@@ -59,12 +60,12 @@ std::vector<std::pair<expression, expression>> stark_dyn()
     const auto xdot = vx;
     const auto ydot = vy;
     const auto zdot = vz;
-    const auto vxdot = -mu * pow(r2, -3. / 2) * x +  ux / m;
-    const auto vydot = -mu * pow(r2, -3. / 2) * y +  uy / m;
-    const auto vzdot = -mu * pow(r2, -3. / 2) * z +  uz / m;
+    const auto vxdot = -mu * pow(r2, -3. / 2) * x + ux / m;
+    const auto vydot = -mu * pow(r2, -3. / 2) * y + uy / m;
+    const auto vzdot = -mu * pow(r2, -3. / 2) * z + uz / m;
     // To avoid singularities in the corner case u_norm=0. we use a select here. Implications on performances should be
     // studied.
-    const auto mdot =  select(eq(u_norm, 0.), 0., - u_norm / veff);
+    const auto mdot = select(eq(u_norm, 0.), 0., -u_norm / veff);
     return {prime(x) = xdot,   prime(y) = ydot,   prime(z) = zdot, prime(vx) = vxdot,
             prime(vy) = vydot, prime(vz) = vzdot, prime(m) = mdot};
 };
@@ -83,7 +84,8 @@ const heyoka::taylor_adaptive<double> &get_ta_stark(double tol)
     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>{stark_dyn(), init_state, heyoka::kw::tol = tol};
+        auto new_ta = taylor_adaptive<double>{stark_dyn(), init_state, heyoka::kw::tol = tol,
+                                              heyoka::kw::pars = {1., 1., 0., 0., 0.}};
         return ta_stark_cache.insert(std::make_pair(tol, std::move(new_ta))).first->second;
     } else {
         // Cache hit, return existing.
@@ -107,7 +109,8 @@ const heyoka::taylor_adaptive<double> &get_ta_stark_var(double tol)
         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};
+        auto new_ta = taylor_adaptive<double>{vsys, init_state, heyoka::kw::tol = tol, heyoka::kw::compact_mode = true,
+                                              heyoka::kw::pars = {1., 1., 0., 0., 1e-32}};
         return ta_stark_var_cache.insert(std::make_pair(tol, std::move(new_ta))).first->second;
     } else {
         // Cache hit, return existing.