From 61940f994e3fb2ea172931c3c44fffb4c769f105 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Thu, 26 Sep 2024 11:24:42 +0200 Subject: [PATCH] minovitch flybys added --- CMakeLists.txt | 1 + doc/api.rst | 1 + doc/flyby.rst | 20 ++++ include/kep3/core_astro/eq2par2eq.hpp | 6 -- include/kep3/core_astro/flyby.hpp | 26 +++-- include/kep3/core_astro/ic2eq2ic.hpp | 6 -- include/kep3/core_astro/ic2par2ic.hpp | 6 -- include/kep3/linalg.hpp | 5 +- pykep/core.cpp | 21 ++++ pykep/docstrings.cpp | 147 ++++++++++++++++++++++++-- pykep/docstrings.hpp | 5 + src/core_astro/flyby.cpp | 102 ++++++++++++++++++ src/linalg.cpp | 13 +++ src/udpla/keplerian.cpp | 6 +- test/CMakeLists.txt | 3 +- test/flyby_test.cpp | 125 ++++++++++++++++++++++ test/ic2eq2ic_test.cpp | 8 +- test/test_helpers.hpp | 2 +- 18 files changed, 453 insertions(+), 50 deletions(-) create mode 100644 doc/flyby.rst create mode 100644 src/core_astro/flyby.cpp create mode 100644 test/flyby_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 8f769f6..515e8a8 100644 --- a/CMakeLists.txt +++ b/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" diff --git a/doc/api.rst b/doc/api.rst index 19b3da1..9e1cb3f 100644 --- a/doc/api.rst +++ b/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 diff --git a/doc/flyby.rst b/doc/flyby.rst new file mode 100644 index 0000000..be42517 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/include/kep3/core_astro/eq2par2eq.hpp b/include/kep3/core_astro/eq2par2eq.hpp index 9315e27..43adc7d 100644 --- a/include/kep3/core_astro/eq2par2eq.hpp +++ b/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 diff --git a/include/kep3/core_astro/flyby.hpp b/include/kep3/core_astro/flyby.hpp index eb57136..7a4604a 100644 --- a/include/kep3/core_astro/flyby.hpp +++ b/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 fb_con(const std::array &v_rel_in, + const std::array &v_rel_out, double mu, double safe_radius); + kep3_DLL_PUBLIC std::pair fb_con(const std::array &v_rel_in, - const std::array &v_rel_out, double, double); + const std::array &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 &v_rel_in, const std::array &v_rel_out, double mu, + double safe_radius); -kep3_DLL_PUBLIC double fb_dv(const std::array &v_rel_in, const std::array &v_rel_out, double, - double); +kep3_DLL_PUBLIC double fb_dv(const std::array &v_rel_in, const std::array &v_rel_out, + const kep3::planet &pl); -kep3_DLL_PUBLIC double fb_vout(const std::array &v_in, const std::array &v_pla, double rp, - double beta, double mu); +// Propagates an absolute incoming velocity through a fly-by. +kep3_DLL_PUBLIC std::array fb_vout(const std::array &v_in, const std::array &v_pla, + double rp, double beta, double mu); } // namespace kep3 #endif // kep3_FLYBY_H \ No newline at end of file diff --git a/include/kep3/core_astro/ic2eq2ic.hpp b/include/kep3/core_astro/ic2eq2ic.hpp index a056b61..4dc6afb 100644 --- a/include/kep3/core_astro/ic2eq2ic.hpp +++ b/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 diff --git a/include/kep3/core_astro/ic2par2ic.hpp b/include/kep3/core_astro/ic2par2ic.hpp index 885e22a..0ec8816 100644 --- a/include/kep3/core_astro/ic2par2ic.hpp +++ b/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 diff --git a/include/kep3/linalg.hpp b/include/kep3/linalg.hpp index a0175a5..dae1223 100644 --- a/include/kep3/linalg.hpp +++ b/include/kep3/linalg.hpp @@ -45,11 +45,14 @@ xt::xtensor_fixed> _dot(const xt::xtensor_fixed types. +std::array _cross(const std::array &v1, const std::array &v2); +void _normalize(std::array &v1); + } // namespace kep3::linalg #endif diff --git a/pykep/core.cpp b/pykep/core.cpp index fd5af84..5d80f88 100644 --- a/pykep/core.cpp +++ b/pykep/core.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -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 &, 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, 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 &, 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, 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_ sims_flanagan(m, "_sims_flanagan", pykep::leg_sf_docstring().c_str()); sims_flanagan diff --git a/pykep/docstrings.cpp b/pykep/docstrings.cpp index a18cc3b..9f1cf99 100644 --- a/pykep/docstrings.cpp +++ b/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 \ No newline at end of file diff --git a/pykep/docstrings.hpp b/pykep/docstrings.hpp index 945604c..a3163cc 100644 --- a/pykep/docstrings.hpp +++ b/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(); diff --git a/src/core_astro/flyby.cpp b/src/core_astro/flyby.cpp new file mode 100644 index 0000000..4bef7f9 --- /dev/null +++ b/src/core_astro/flyby.cpp @@ -0,0 +1,102 @@ +// 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 + +namespace kep3 +{ + +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +std::pair fb_con(const std::array &v_rel_in, const std::array &v_rel_out, + double mu, double safe_radius) +{ + const double Vin2 = v_rel_in[0] * v_rel_in[0] + v_rel_in[1] * v_rel_in[1] + v_rel_in[2] * v_rel_in[2]; + const double Vout2 = v_rel_out[0] * v_rel_out[0] + v_rel_out[1] * v_rel_out[1] + v_rel_out[2] * v_rel_out[2]; + const double eq_V2 = Vin2 - Vout2; + + const double e_min = 1 + safe_radius / mu * Vin2; + const double alpha + = std::acos((v_rel_in[0] * v_rel_out[0] + v_rel_in[1] * v_rel_out[1] + v_rel_in[2] * v_rel_out[2]) + / std::sqrt(Vin2 * Vout2)); + const double ineq_delta = alpha - 2 * std::asin(1 / e_min); + return {eq_V2, ineq_delta}; +} + +std::pair fb_con(const std::array &v_rel_in, const std::array &v_rel_out, + const kep3::planet &pl) +{ + return std::move(fb_con(v_rel_in, v_rel_out, pl.get_mu_self(), pl.get_safe_radius())); +} + +double fb_dv(const std::array &v_rel_in, const std::array &v_rel_out, double mu, + double safe_radius) +{ + const double Vin2 = v_rel_in[0] * v_rel_in[0] + v_rel_in[1] * v_rel_in[1] + v_rel_in[2] * v_rel_in[2]; + const double Vout2 = v_rel_out[0] * v_rel_out[0] + v_rel_out[1] * v_rel_out[1] + v_rel_out[2] * v_rel_out[2]; + // eq_V2 = Vin2 - Vout2; + + const double e_min = 1 + safe_radius / mu * Vin2; + const double alpha + = std::acos((v_rel_in[0] * v_rel_out[0] + v_rel_in[1] * v_rel_out[1] + v_rel_in[2] * v_rel_out[2]) + / std::sqrt(Vin2 * Vout2)); + const double ineq_delta = alpha - 2 * std::asin(1 / e_min); + + double dv = 0.; + if (ineq_delta > 0.0) { + dv = std::sqrt(Vout2 + Vin2 - 2.0 * std::sqrt(Vout2 * Vin2) * std::cos(ineq_delta)); + } else { + dv = std::abs(sqrt(Vout2) - std::sqrt(Vin2)); + } + return dv; +} + +double fb_dv(const std::array &v_rel_in, const std::array &v_rel_out, const kep3::planet &pl) +{ + return fb_dv(v_rel_in, v_rel_out, pl.get_mu_self(), pl.get_safe_radius()); +} + +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +std::array fb_vout(const std::array &v_in, const std::array &v_pla, double rp, + double beta, double mu) +{ + std::array v_rel_in = {v_in[0] - v_pla[0], v_in[1] - v_pla[1], v_in[2] - v_pla[2]}; + std::array v_out{}; + + const double v_rel_in2 = v_rel_in[0] * v_rel_in[0] + v_rel_in[1] * v_rel_in[1] + v_rel_in[2] * v_rel_in[2]; + const double v_rel_in_norm = std::sqrt(v_rel_in2); + const double ecc = 1 + rp / mu * v_rel_in2; + const double delta = 2 * std::asin(1.0 / ecc); + const std::array i_hat + = {{v_rel_in[0] / v_rel_in_norm, v_rel_in[1] / v_rel_in_norm, v_rel_in[2] / v_rel_in_norm}}; + + std::array j_hat = linalg::_cross(i_hat, v_pla); + linalg::_normalize(j_hat); + std::array k_hat = linalg::_cross(i_hat, j_hat); + + v_out[0] = v_pla[0] + v_rel_in_norm * std::cos(delta) * i_hat[0] + + v_rel_in_norm * std::cos(beta) * std::sin(delta) * j_hat[0] + + v_rel_in_norm * std::sin(beta) * std::sin(delta) * k_hat[0]; + v_out[1] = v_pla[1] + v_rel_in_norm * std::cos(delta) * i_hat[1] + + v_rel_in_norm * std::cos(beta) * std::sin(delta) * j_hat[1] + + v_rel_in_norm * std::sin(beta) * std::sin(delta) * k_hat[1]; + v_out[2] = v_pla[2] + v_rel_in_norm * std::cos(delta) * i_hat[2] + + v_rel_in_norm * std::cos(beta) * std::sin(delta) * j_hat[2] + + v_rel_in_norm * std::sin(beta) * std::sin(delta) * k_hat[2]; + return v_out; +} + +} // namespace kep3 diff --git a/src/linalg.cpp b/src/linalg.cpp index 73ae88b..3ba3873 100644 --- a/src/linalg.cpp +++ b/src/linalg.cpp @@ -21,4 +21,17 @@ mat31 _cross(const mat31 &v1, const mat31 &v2) return {{v1(1, 0) * v2(2, 0) - v2(1, 0) * v1(2, 0), v2(0, 0) * v1(2, 0) - v1(0, 0) * v2(2, 0), v1(0, 0) * v2(1, 0) - v2(0, 0) * v1(1, 0)}}; } + +std::array _cross(const std::array &v1, const std::array &v2) +{ + return {{v1[1] * v2[2] - v2[1] * v1[2], v2[0] * v1[2] - v1[0] * v2[2], v1[0] * v2[1] - v2[0] * v1[1]}}; +} + +void _normalize(std::array &v1) +{ + const double norm = std::sqrt(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]); + v1[0] /= norm; + v1[1] /= norm; + v1[2] /= norm; +} } // namespace kep3::linalg diff --git a/src/udpla/keplerian.cpp b/src/udpla/keplerian.cpp index a26304d..7707cec 100644 --- a/src/udpla/keplerian.cpp +++ b/src/udpla/keplerian.cpp @@ -36,11 +36,11 @@ keplerian::keplerian(const epoch &ref_epoch, const std::array 0) ? m_ellipse = false : m_ellipse = true; - double a = -m_mu_central_body / 2. / en; + const double a = -m_mu_central_body / 2. / en; if (m_ellipse) { m_period = kep3::pi * 2. * std::sqrt(a * a * a / m_mu_central_body); } else { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5c7d9c0..8c3bfe7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -42,4 +42,5 @@ ADD_kep3_TESTCASE(lambert_problem_test) ADD_kep3_TESTCASE(stark_problem_test) ADD_kep3_TESTCASE(leg_sims_flanagan_test) ADD_kep3_TESTCASE(ta_stark_test) -ADD_kep3_TESTCASE(ta_cr3bp_test) \ No newline at end of file +ADD_kep3_TESTCASE(ta_cr3bp_test) +ADD_kep3_TESTCASE(flyby_test) \ No newline at end of file diff --git a/test/flyby_test.cpp b/test/flyby_test.cpp new file mode 100644 index 0000000..2e6f0d7 --- /dev/null +++ b/test/flyby_test.cpp @@ -0,0 +1,125 @@ +// 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 "catch.hpp" +#include "kep3/core_astro/constants.hpp" +#include "test_helpers.hpp" + +TEST_CASE("fb_con") +{ + const double mu = 398600441800000.0; + const double safe_radius = 7015800.000000001; + // Call from doubles + { + const std::array v_rel_in = {7000., 0., 0.}; + const std::array v_rel_out = {7100, 0., 0.}; + + const std::pair ground_truth{-1410000.0, -1.1335787944199696}; + auto res = kep3::fb_con(v_rel_in, v_rel_out, mu, safe_radius); + REQUIRE(kep3_tests::floating_point_error(res.first, ground_truth.first) < 1e-13); + REQUIRE(kep3_tests::floating_point_error(res.second, ground_truth.second) < 1e-13); + } + { + const std::array v_rel_in = {7200., -4567.7655, 1234.4233}; + const std::array v_rel_out = {7100., 220.123, -144.432}; + const std::pair ground_truth{23748967.80882015, -0.2801587391727032}; + auto res = kep3::fb_con(v_rel_in, v_rel_out, mu, safe_radius); + REQUIRE(kep3_tests::floating_point_error(res.first, ground_truth.first) < 1e-14); + REQUIRE(kep3_tests::floating_point_error(res.second, ground_truth.second) < 1e-14); + } + { + // Call from planet + // We use 2000-01-01 as a reference epoch for all these tests + const double ref_epoch = 0.; + // This is a circular orbit at 1 AU. + const std::array, 2> pos_vel_0{{{kep3::AU, 0., 0.}, {0., kep3::EARTH_VELOCITY, 0.}}}; + // A keplerian planet orbiting the Sun on such a perfectly circular orbit. + const kep3::udpla::keplerian udpla{ + kep3::epoch(ref_epoch), pos_vel_0, kep3::MU_SUN, "test", {mu, safe_radius, safe_radius}}; + const kep3::planet pl_test{udpla}; + const std::array v_rel_in = {7200., -4567.7655, 1234.4233}; + const std::array v_rel_out = {7100., 220.123, -144.432}; + const std::pair ground_truth{23748967.80882015, -0.2801587391727032}; + auto res = kep3::fb_con(v_rel_in, v_rel_out, pl_test); + REQUIRE(kep3_tests::floating_point_error(res.first, ground_truth.first) < 1e-14); + REQUIRE(kep3_tests::floating_point_error(res.second, ground_truth.second) < 1e-14); + } +} + +TEST_CASE("fb_dv") +{ + const double mu = 398600441800000.0; + const double safe_radius = 7015800.000000001; + // Call from doubles + { + const std::array v_rel_in = {7000., 0., 0.}; + const std::array v_rel_out = {7100, 0., 0.}; + const double ground_truth = 100.; + auto res = kep3::fb_dv(v_rel_in, v_rel_out, mu, safe_radius); + REQUIRE(kep3_tests::floating_point_error(res, ground_truth) < 1e-13); + } + { + const std::array v_rel_in = {7200., -4567.7655, 1234.4233}; + const std::array v_rel_out = {7100., 220.123, -144.432}; + const double ground_truth = 1510.704060449003; + auto res = kep3::fb_dv(v_rel_in, v_rel_out, mu, safe_radius); + REQUIRE(kep3_tests::floating_point_error(res, ground_truth) < 1e-14); + } + { + // Call from planet + // We use 2000-01-01 as a reference epoch for all these tests + const double ref_epoch = 0.; + // This is a circular orbit at 1 AU. + const std::array, 2> pos_vel_0{{{kep3::AU, 0., 0.}, {0., kep3::EARTH_VELOCITY, 0.}}}; + // A keplerian planet orbiting the Sun on such a perfectly circular orbit. + const kep3::udpla::keplerian udpla{ + kep3::epoch(ref_epoch), pos_vel_0, kep3::MU_SUN, "test", {mu, safe_radius, safe_radius}}; + const kep3::planet pl_test{udpla}; + const std::array v_rel_in = {7200., -4567.7655, 1234.4233}; + const std::array v_rel_out = {7100., 220.123, -144.432}; + const double ground_truth = 1510.704060449003; + auto res = kep3::fb_dv(v_rel_in, v_rel_out, mu, safe_radius); + REQUIRE(kep3_tests::floating_point_error(res, ground_truth) < 1e-14); + } +} + +TEST_CASE("fb_vout") +{ + { + const double rp = 2.; + const double beta = 3.1415 / 3.; + const double mu = 1.; + const std::array v_in = {1., 0., 0.}; + const std::array v_pla = {0., 1., 0.}; + const std::array ground_truth = {0.5805947972994727, -0.2594052027005272, 0.27714295365321495}; + auto res = kep3::fb_vout(v_in, v_pla, rp, beta, mu); + REQUIRE(kep3_tests::floating_point_error_vector(ground_truth, res) < 1e-14); + } + { + const double rp = 7600000.; + const double beta = 3.1415 / 3. * 2.; + const double mu = kep3::MU_EARTH; + const std::array v_in = {20000., 1234.5678, -12098.123565234}; + const std::array v_pla = {25000.243234322, 1., 0.}; + const std::array ground_truth = {15773.610960458058, 3862.490616725501, -8535.037364832286}; + auto res = kep3::fb_vout(v_in, v_pla, rp, beta, mu); + REQUIRE(kep3_tests::floating_point_error_vector(ground_truth, res) < 1e-14); + } +} \ No newline at end of file diff --git a/test/ic2eq2ic_test.cpp b/test/ic2eq2ic_test.cpp index 525cc4a..37b3cf6 100644 --- a/test/ic2eq2ic_test.cpp +++ b/test/ic2eq2ic_test.cpp @@ -13,8 +13,7 @@ #include #include -#include - +#include #include #include @@ -24,10 +23,9 @@ using Catch::Matchers::WithinRel; using kep3::eq2ic; using kep3::ic2eq; +using kep3::pi; -constexpr double pi{boost::math::constants::pi()}; - -TEST_CASE("ic2eq") + TEST_CASE("fb_con") { // Zero inclination and eccentricity { diff --git a/test/test_helpers.hpp b/test/test_helpers.hpp index 1c5f331..4243aaa 100644 --- a/test/test_helpers.hpp +++ b/test/test_helpers.hpp @@ -33,7 +33,7 @@ template double L_infinity_norm(T a, T b) { if (a.size() != b.size()) { - throw std::domain_error("Computing the L-infinity norm of two vectors having unequal size."); + throw std::domain_error("Subtracting two vectors having unequal size."); } double retval = 0.;