minovitch flybys added

CMakeLists.txt

@@ -145,6 +145,7 @@ set(kep3_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/udpla/jpl_lp.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/udpla/vsop2013.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/leg/sims_flanagan.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/flyby.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"

doc/api.rst

@@ -14,6 +14,7 @@ pykep API is designed to maximize its usability. Let us know what you think abou
    elements
    epoch
    lambert
+   flyby
    leg
    planet
    plot

doc/flyby.rst

@@ -0,0 +1,20 @@
+.. _flyby:
+
+Fly-by routines
+########################
+The gravity assist technique, widely used in the Pioneer and Voyager missions, was developed by Michael Andrew Minovitch when he was a UCLA 
+graduate student working as intern at NASA's Jet Propulsion Laboratory. Many authors later refined the technique which, avoiding to solve
+directly the full three body problem, is based on patching the planetocentric hyperbola with incoming and outgoing (elliptic) trajectories.
+In `pykep` we offer the basic routines that allow to use this technqiue in the context of trajectory optimization and patched conics propagations.
+
+.. currentmodule:: pykep
+
+------------------------------------------------------
+
+Minovitch fly-by
+****************
+.. autofunction:: fb_con
+
+.. autofunction:: fb_dv
+
+.. autofunction:: fb_vout 

include/kep3/core_astro/eq2par2eq.hpp

@@ -7,12 +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/.
 
-/// From cartesian to osculating Keplerian
-/**
- * Transforms cartesian coordinates (r,v) to Keplerian elements (a,e,i,W,w,E).
- * Note that we use the eccentric anomaly (or Gudermannian if e > 1)
- */
-
 #ifndef kep3_EQ2PAR2EQ_H
 #define kep3_EQ2PAR2EQ_H
 

include/kep3/core_astro/flyby.hpp

@@ -7,12 +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/.
 
-/// From cartesian to osculating Keplerian
-/**
- * Transforms cartesian coordinates (r,v) to Keplerian elements (a,e,i,W,w,E).
- * Note that we use the eccentric anomaly (or Gudermannian if e > 1)
- */
-
 #ifndef kep3_FLYBY_H
 #define kep3_FLYBY_H
 
@@ -25,14 +19,24 @@
 namespace kep3
 {
 
+// Returns the constraints [v2, alpha] of a fly-by. The first is an equality constraint, the second an inequality
+// (negative if satisfied).
+kep3_DLL_PUBLIC std::pair<double, double> fb_con(const std::array<double, 3> &v_rel_in,
+                                                 const std::array<double, 3> &v_rel_out, double mu, double safe_radius);
+
 kep3_DLL_PUBLIC std::pair<double, double> fb_con(const std::array<double, 3> &v_rel_in,
-                                                 const std::array<double, 3> &v_rel_out, double, double);
+                                                 const std::array<double, 3> &v_rel_out, const kep3::planet &pl);
+
+// Returns the dv needed to make a fly-by feasible. (assuming one DV at the out conditions).
+kep3_DLL_PUBLIC double fb_dv(const std::array<double, 3> &v_rel_in, const std::array<double, 3> &v_rel_out, double mu,
+                             double safe_radius);
 
-kep3_DLL_PUBLIC double fb_dv(const std::array<double, 3> &v_rel_in, const std::array<double, 3> &v_rel_out, double,
-                             double);
+kep3_DLL_PUBLIC double fb_dv(const std::array<double, 3> &v_rel_in, const std::array<double, 3> &v_rel_out,
+                             const kep3::planet &pl);
 
-kep3_DLL_PUBLIC double fb_vout(const std::array<double, 3> &v_in, const std::array<double, 3> &v_pla, double rp,
-                               double beta, double mu);
+// Propagates an absolute incoming velocity through a fly-by.
+kep3_DLL_PUBLIC std::array<double, 3> fb_vout(const std::array<double, 3> &v_in, const std::array<double, 3> &v_pla,
+                                              double rp, double beta, double mu);
 
 } // namespace kep3
 #endif // kep3_FLYBY_H

include/kep3/core_astro/ic2eq2ic.hpp

@@ -7,12 +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/.
 
-/// From cartesian to osculating Keplerian
-/**
- * Transforms cartesian coordinates (r,v) to Keplerian elements (a,e,i,W,w,E).
- * Note that we use the eccentric anomaly (or Gudermannian if e > 1)
- */
-
 #ifndef kep3_IC2EQ2IC_H
 #define kep3_IC2EQ2IC_H
 

include/kep3/core_astro/ic2par2ic.hpp

@@ -7,12 +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/.
 
-/// From cartesian to osculating Keplerian
-/**
- * Transforms cartesian coordinates (r,v) to Keplerian elements (a,e,i,W,w,E).
- * Note that we use the eccentric anomaly (or Gudermannian if e > 1)
- */
-
 #ifndef kep3_IC2PAR2IC_H
 #define kep3_IC2PAR2IC_H
 

include/kep3/linalg.hpp

@@ -45,11 +45,14 @@ xt::xtensor_fixed<double, xt::xshape<m, n>> _dot(const xt::xtensor_fixed<double,
     }
     return C;
 }
-
 mat33 _skew(const mat31 &v);
 mat31 _cross(const mat31 &v1, const mat31 &v2);
 // ---------------------------------------------------------------------------------------
 
+// Linear algebra helpers for std::array<double, 3> types.
+std::array<double, 3> _cross(const std::array<double, 3> &v1, const std::array<double, 3> &v2);
+void _normalize(std::array<double, 3> &v1);
+
 } // namespace kep3::linalg
 
 #endif

pykep/core.cpp

@@ -14,6 +14,7 @@
 #include <kep3/core_astro/constants.hpp>
 #include <kep3/core_astro/convert_anomalies.hpp>
 #include <kep3/core_astro/eq2par2eq.hpp>
+#include <kep3/core_astro/flyby.hpp>
 #include <kep3/core_astro/ic2eq2ic.hpp>
 #include <kep3/core_astro/ic2par2ic.hpp>
 #include <kep3/core_astro/propagate_lagrangian.hpp>
@@ -438,6 +439,26 @@ PYBIND11_MODULE(core, m)
             py::arg("rvm_state"), py::arg("thrust"), py::arg("tof"),
             pykep::stark_problem_propagate_var_docstring().c_str());
 
+    // Exposing fly-by routines
+    m.def("fb_con",
+          py::overload_cast<const std::array<double, 3> &, const std::array<double, 3> &, const kep3::planet &>(
+              &kep3::fb_con),
+          py::arg("v_rel_in"), py::arg("v_rel_out"), py::arg("planet"), pykep::fb_con_docstring().c_str());
+    m.def(
+        "fb_con",
+        py::overload_cast<const std::array<double, 3> &, const std::array<double, 3> &, double, double>(&kep3::fb_con),
+        py::arg("v_rel_in"), py::arg("v_rel_out"), py::arg("mu"), py::arg("safe_radius"));
+
+    m.def("fb_dv",
+          py::overload_cast<const std::array<double, 3> &, const std::array<double, 3> &, const kep3::planet &>(
+              &kep3::fb_dv),
+          py::arg("v_rel_in"), py::arg("v_rel_out"), py::arg("planet"), pykep::fb_dv_docstring().c_str());
+    m.def("fb_dv",
+          py::overload_cast<const std::array<double, 3> &, const std::array<double, 3> &, double, double>(&kep3::fb_dv),
+          py::arg("v_rel_in"), py::arg("v_rel_out"), py::arg("mu"), py::arg("safe_radius"));
+    m.def("fb_vout", &kep3::fb_vout, py::arg("v_in"), py::arg("v_pla"), py::arg("rp"), py::arg("beta"), py::arg("mu"),
+          pykep::fb_vout_docstring().c_str());
+
     // Exposing the sims_flanagan leg
     py::class_<kep3::leg::sims_flanagan> sims_flanagan(m, "_sims_flanagan", pykep::leg_sf_docstring().c_str());
     sims_flanagan

pykep/docstrings.cpp

@@ -1059,12 +1059,12 @@ Constructs a Keplerian udpla from its position and velocity at epoch.
 
 std::string udpla_keplerian_from_elem_docstring()
 {
-    return R"(__init__(ep, elem, mu_central_body, name = "unkown", added_params = [-1,-1,-1], elem_type = KEP_F)
+    return R"(__init__(when, elem, mu_central_body, name = "unkown", added_params = [-1,-1,-1], elem_type = KEP_F)
 
 Constructs a Keplerian udpla from its orbital elements at epoch.
 
 Args:
-    *ep* (:class:`~pykep.epoch`): the epoch at which the orbital elements are provided.
+    *when* (:class:`~pykep.epoch`): the epoch at which the orbital elements are provided.
 
     *elem* (:class:`list` or :class:`numpy.ndarray`): the orbital elements. by default.
 
@@ -1074,13 +1074,13 @@ Constructs a Keplerian udpla from its orbital elements at epoch.
 
     *added_params* (:class:`list`): the body gravitational parameter, its radius and its safe radius. (if -1 they are assumed unkown)
 
-    *el_type* (:class:`~pykep.el_type`): the elements type. Deafulets to osculating Keplrian (a ,e ,i, W, w, f) with true anomaly.
+    *el_type* (:class:`~pykep.el_type`): the elements type. Deafulets to osculating Keplerian (a ,e ,i, W, w, f) with true anomaly.
 
 Examples:
     >>> import pykep as pk
     >>> elem = [1, 0, 0, 0, 0, 0]
-    >>> ep = pk.epoch("2025-03-22")
-    >>> udpla = pk.udpla.keplerian(ep = ep, elem = elem, mu_central_body =1, name = "my_pla")
+    >>> when = pk.epoch("2025-03-22")
+    >>> udpla = pk.udpla.keplerian(when = when, elem = elem, mu_central_body =1, name = "my_pla")
     >>> pla = pk.planet(udpla)
 )";
 }
@@ -1349,9 +1349,10 @@ The dynamics is that returned by :func:`~pykep.ta.stark_dyn`: and also used in :
   >>> ta.propagate_until(tof)
 )";
 }
-   
+
 std::string stark_dyn_docstring()
-{return R"(stark_dyn()
+{
+    return R"(stark_dyn()
 
 The dynamics of the Stark problem. 
 
@@ -1437,9 +1438,10 @@ The specific dynamics used is that returned by :func:`~pykep.ta.cr3bp_dyn`.
   >>> ta.propagate_until(tof)
 )";
 }
-   
+
 std::string cr3bp_dyn_docstring()
-{return R"(cr3bp_dyn()
+{
+    return R"(cr3bp_dyn()
 
 The dynamics of the Circular Restricted Three Body Problem (CR3BP).
 
@@ -1467,7 +1469,6 @@ where :math:`\mu` is the only parameter.
 )";
 }
 
-
 std::string propagate_lagrangian_docstring()
 {
     return R"(propagate_lagrangian(rv = [[1,0,0], [0,1,0]], tof = pi/2, mu = 1, stm = False)
@@ -1720,4 +1721,130 @@ Computes the gradients of the throttles constraints. Indicating the total time o
 )";
 };
 
+std::string fb_con_docstring()
+{
+    return R"(fb_con(v_rel_in, v_rel_out, mu, safe_radius)
+Alternative signature: fb_con(v_rel_in, v_rel_out, planet)
+
+Computes the constraint violation during a fly-by modelled as an instantaneous rotation of the incoming and outgoing relative velocities (Mivovitch).
+The two must be identical in magnitude (equality constraint) and the angle :math:`\alpha` between them must be less than the 
+:math:`\alpha_{max}`: the maximum value allowed for that particular *planet* (inequality constraint), as computed from its gravitational
+parameter and safe radius using the formula:
+
+.. math::
+  \alpha_{max} = - 2 \arcsin\left(\frac{1}{e_{min}}\right)
+
+where:
+
+.. math::
+  e_{min} = 1 + V_{\infty}^2 \frac{R_{safe}}{\mu}
+
+.. note::
+  This function is often used in the multiple gravity assist low-thrust (MGA-LT) encoding of an interplanetary trajectory where multiple
+  low-thrust legs are patched at the fly-by planets by forcing satisfaction for these constraints.
+
+Args:
+  *v_rel_in* (:class:`list` (3,)): Cartesian components of the incoming relative velocity.
+
+  *v_rel_out* (:class:`list` (3,)): Cartesian components of the outgoing relative velocity.
+
+  *mu* (:class:`float`): planet gravitational parameter
+
+  *safe_radius* (:class:`float`): planet safe radius
+
+  *planet* (:class:`~pykep.planet`): planet (in which case *mu* and *safe_radius* will be extracted from this object). Note: this signature is slower and to be avoided.
+
+Returns:
+  :class:`tuple` [:class:`float`, :class:`float`]: The equality constraint violation (defined as the difference between the squared velocities) and the inequality constraint violation 
+  (negative if satisfied).
+
+Examples:
+  >>> import pykep as pk
+  >>> eq, ineq = pk.fb_con(v_rel_in = [10.,1.,-4.], v_rel_out = [10.,1.,-4.], mu = 1., safe_radius = 1.)
+
+)";
+};
+
+std::string fb_dv_docstring()
+{
+    return R"(fb_dv(v_rel_in, v_rel_out, mu, safe_radius)
+Alternative signature: fb_dv(v_rel_in, v_rel_out, planet)
+
+Computes the :math:`\Delta V` necessary to perform a fly-by modelled as an instantaneous rotation of the incoming and outgoing
+relative velocities (Mivovitch). If planetary gravity is not enough to patch the incoming and outcoming conditions
+(i.e. the :func:`~pykep.fb_con()` returns some constraint violation) a :math:`\Delta V` is assumed
+at the end of the planetocentric hyperbola.
+
+.. note::
+  This function is often used in the multiple gravity assist (MGA) encoding of an interplanetary trajectory where multiple
+  Lambert arcs are patched at the fly-by planets by applying a hopefully vanishing :math:`\Delta V`.
+
+Args:
+  *v_rel_in* (:class:`list` (3,)): Cartesian components of the incoming relative velocity.
+
+  *v_rel_out* (:class:`list` (3,)): Cartesian components of the outgoing relative velocity.
+
+  *mu* (:class:`float`): planet gravitational parameter
+
+  *safe_radius* (:class:`float`): planet safe radius
+
+  *planet* (:class:`~pykep.planet`): planet (in which case *mu* and *safe_radius* will be extracted from this object). Note: this signature is slower and to be avoided.
+
+Returns:
+  :class:`float`: The magnitude of the required :math:`\Delta V`
+
+Examples:
+  >>> import pykep as pk
+  >>> DV = pk.fb_dv(v_rel_in = [10.,1.,-4.], v_rel_out = [10.,1.,-4.], mu = 1., safe_radius = 1.)
+)";
+};
+
+std::string fb_vout_docstring()
+{
+    return R"(fb_vout(v_in, v_pla, rp, beta, mu)
+
+Propagates incoming conditions through a planetary encounter (fly-by) assuming an instantaneous rotation of magnitude :math:`\delta` of 
+the incoming and outgoing relative velocities (Mivovitch). The planetocentric hyperbola (hence the angle :math:`\delta`) is fully
+determined by the incoming condition :math:`v_{\infty} = \mathbf v_{in} - \mathbf v_{pla}` as well as by its pericenter :math:`r_p`
+and an angle :math:`\beta` defining the orientation of the orbital plane. Eventually the outgoing conditions are computed as:
+
+.. math::
+  \mathbf v_{out} = \mathbf v_{in} + \mathbf v_{\infty}^{out}
+
+.. math::
+  \mathbf v_{\infty}^{out} = |\mathbf v_{\infty}^{in}|
+  \left(
+  \cos\delta\hat{\mathbf b}_1
+  +\sin\delta\cos\beta\hat{\mathbf b}_2
+  +\sin\delta\sin\beta\hat{\mathbf b}_3
+  \right)
+
+where the :math:`[\hat{\mathbf b}_1, \hat{\mathbf b}_2, \hat{\mathbf b}_3]` frame is defined by the incoming 
+relative velocity (normalized), :math:`\hat{\mathbf b}_1 \times \mathbf v_{pla}` (normalized) and completing the right handed frame.
+
+.. note::
+  This function is often used in the multiple gravity assist with DSM (MGA-DSM) encoding of an interplanetary trajectory where multiple
+  impulsive manouvres are allowed betwwen planetary fly-bys.
+
+Args:
+  *v_in* (:class:`list` (3,)): Cartesian components of the incoming velocity 
+  (note: this is NOT the relative velocity, rather the absolute and in the same frame as *v_pla*).
+
+  *v_pla* (:class:`list` (3,)): Cartesian components of the planet velocity.
+
+  *rp* (:class:`float`): planetocentric hyperbola pericenter radius.
+
+  *beta* (:class:`float`): planetocentric hyperbola plane angle.
+
+  *mu* (:class:`float`): planet gravitational parameter.
+
+Returns:
+  :class:`list` (3,): The outgoing velocity in the frame of *v_pla* (inertial).
+
+Examples:
+  >>> import pykep as pk
+  >>> DV = pk.fb_vout(v_in = [1.,1.,1], v_pla = [10.,1.,-4.], rp = 1., mu = 1., beta=1.2)
+)";
+};
+
 } // namespace pykep

pykep/docstrings.hpp

@@ -89,6 +89,11 @@ std::string cr3bp_dyn_docstring();
 // Lambert Problem
 std::string lambert_problem_docstring();
 
+// Flybys
+std::string fb_con_docstring();
+std::string fb_dv_docstring();
+std::string fb_vout_docstring();
+
 // Stark problem
 std::string stark_problem_docstring();
 std::string stark_problem_propagate_docstring();