diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 2403ffc9..cd8b7449 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -1,4 +1,4 @@ -name: GitHub CI +name: C++ Library on: push: branches: diff --git a/.github/workflows/pykep.yml b/.github/workflows/pykep.yml new file mode 100644 index 00000000..05cd3627 --- /dev/null +++ b/.github/workflows/pykep.yml @@ -0,0 +1,17 @@ +name: Python Module +on: + push: + branches: + - main + tags: + - 'v[0-9]+.[0-9]+.[0-9]+' + pull_request: + branches: + - main +jobs: + linux_pykep: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - name: Build + run: bash tools/gha_linux_pykep.sh diff --git a/.gitignore b/.gitignore index f0f5eb89..6a63076d 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,6 @@ build* # config files include/kep3/config.hpp + +# python stuff +*__pycache__* diff --git a/CMakeLists.txt b/CMakeLists.txt index f836f9ff..5297f584 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -165,9 +165,12 @@ target_include_directories(kep3 PUBLIC $ $) # Boost. -# NOTE: need 1.73 for atomic_ref. -find_package(Boost 1.73 REQUIRED) +# NOTE: need 1.69 for safe numerics. +set(_kep3_MIN_BOOST_VERSION "1.69") +find_package(Boost ${_kep3_MIN_BOOST_VERSION} REQUIRED) target_link_libraries(kep3 PUBLIC Boost::boost) +# NOTE: quench warnings from Boost when building the library. +target_compile_definitions(kep3 PRIVATE BOOST_ALLOW_DEPRECATED_HEADERS) # fmt. find_package(fmt CONFIG REQUIRED) @@ -207,6 +210,17 @@ install(TARGETS kep3 ARCHIVE DESTINATION "${kep3_INSTALL_LIBDIR}" RUNTIME DESTINATION bin ) +# Setup of the CMake config file. +configure_file("${CMAKE_CURRENT_SOURCE_DIR}/kep3-config.cmake.in" "${CMAKE_CURRENT_BINARY_DIR}/kep3-config.cmake" @ONLY) +install(FILES "${CMAKE_CURRENT_BINARY_DIR}/kep3-config.cmake" + DESTINATION "${kep3_INSTALL_LIBDIR}/cmake/kep3") +install(EXPORT kep3_export NAMESPACE kep3:: DESTINATION "${kep3_INSTALL_LIBDIR}/cmake/kep3") +# Take care of versioning. +include(CMakePackageConfigHelpers) +write_basic_package_version_file("${CMAKE_CURRENT_BINARY_DIR}/kep3-config-version.cmake" COMPATIBILITY SameMinorVersion) +install(FILES "${CMAKE_CURRENT_BINARY_DIR}/kep3-config-version.cmake" DESTINATION "${kep3_INSTALL_LIBDIR}/cmake/kep3") + + if(kep3_BUILD_TESTS) enable_testing() add_subdirectory(test) @@ -214,4 +228,24 @@ endif() if(kep3_BUILD_BENCHMARKS) add_subdirectory(benchmark) -endif() \ No newline at end of file + +endif() + +if(kep3_BUILD_PYTHON_BINDINGS) + # Find python + find_package(Python3 REQUIRED COMPONENTS Interpreter Development) + message(STATUS "Python3 interpreter: ${Python3_EXECUTABLE}") + message(STATUS "Python3 installation directory: ${Python3_SITEARCH}") + set(PYKEP_INSTALL_PATH "" CACHE STRING "pykep module installation path") + mark_as_advanced(PYKEP_INSTALL_PATH) + + # pybind11. + find_package(pybind11 REQUIRED) + + # Finding kep3 (cpp) + find_package(kep3 ${kep3_VERSION} EXACT REQUIRED) + + # Build directory + add_subdirectory("${CMAKE_SOURCE_DIR}/pykep") +endif() + diff --git a/include/kep3/core_astro/constants.hpp b/include/kep3/core_astro/constants.hpp index d150cf18..aab04bd6 100644 --- a/include/kep3/core_astro/constants.hpp +++ b/include/kep3/core_astro/constants.hpp @@ -28,7 +28,7 @@ constexpr double pi = boost::math::constants::pi(); constexpr double half_pi = boost::math::constants::half_pi(); constexpr double AU = 149597870700.0; // Astronomical Unit (m) -constexpr double CAVENDISH = 73.6687e-11; // Cavendish constant +constexpr double CAVENDISH = 73.6687e-11; // Cavendish constant (N M^2 / kg^2) constexpr double MU_SUN = 1.32712440018e20; // Sun's gravitational parameter (m^3/s^2 kg) constexpr double MU_EARTH = diff --git a/include/kep3/core_astro/convert_anomalies.hpp b/include/kep3/core_astro/convert_anomalies.hpp index 2b75ede9..a1f8cd47 100644 --- a/include/kep3/core_astro/convert_anomalies.hpp +++ b/include/kep3/core_astro/convert_anomalies.hpp @@ -10,12 +10,13 @@ #ifndef kep3_CONVERT_ANOMALIES_H #define kep3_CONVERT_ANOMALIES_H -#include "kep3/core_astro/constants.hpp" #include +#include #include #include +#include #include namespace kep3 { @@ -24,6 +25,10 @@ namespace kep3 { // of revolutions. // NOLINTNEXTLINE(bugprone-easily-swappable-parameters) inline double m2e(double M, double ecc) { + if (ecc >= 1) { + throw std::domain_error( + "m2e: Eccentric anomaly is not defined for ecc >=1."); + } // We compute the sin and cos of the mean anomaly which are used also later // for the initial guess (3rd order expansion of Kepler's equation). double sinM = std::sin(M), cosM = std::cos(M); @@ -31,9 +36,7 @@ inline double m2e(double M, double ecc) { // This makes sure that for high value of M no catastrophic cancellation // occurs, as would be the case using std::fmod(M, 2pi) double M_cropped = std::atan2(sinM, cosM); - if (M_cropped < 0) { - M_cropped += 2 * kep3::pi; - } + // The Initial guess follows from a third order expansion of Kepler's // equation. double IG = M_cropped + ecc * sinM + ecc * ecc * sinM * cosM + @@ -57,38 +60,73 @@ inline double m2e(double M, double ecc) { return sol; } // eccentric to mean (only ellipses) e<1 -inline double e2m(double E, double ecc) { return (E - ecc * std::sin(E)); } +inline double e2m(double E, double ecc) { + if (ecc >= 1) { + throw std::domain_error( + "e2m: Eccentric anomaly is not defined for ecc >=1."); + } + return (E - ecc * std::sin(E)); +} // eccentric to true (only ellipses) e<1 (returns in range [-pi,pi]) inline double e2f(double E, double ecc) { + if (ecc >= 1) { + throw std::domain_error( + "e2f: Eccentric anomaly is not defined for ecc >=1."); + } return 2 * std::atan(std::sqrt((1 + ecc) / (1 - ecc)) * std::tan(E / 2)); } // true to eccentric (only ellipses) e<1 (returns in range [-pi,pi]) inline double f2e(double f, double ecc) { + if (ecc >= 1) { + throw std::domain_error( + "f2e: Eccentric anomaly is not defined for ecc >=1."); + } return 2 * std::atan(std::sqrt((1 - ecc) / (1 + ecc)) * std::tan(f / 2)); } // mean to true (only ellipses) e<1 (returns in range [-pi,pi]) -inline double m2f(double M, double ecc) { return e2f(m2e(M, ecc), ecc); } +inline double m2f(double M, double ecc) { + if (ecc >= 1) { + throw std::domain_error("m2f: Mean anomaly is not defined for ecc >=1."); + }; + return e2f(m2e(M, ecc), ecc); +} // true to mean (only ellipses) e<1 (returns in range [-pi,pi]) -inline double f2m(double f, double ecc) { return e2m(f2e(f, ecc), ecc); } +inline double f2m(double f, double ecc) { + if (ecc >= 1) { + throw std::domain_error("f2m: Mean anomaly is not defined for ecc >=1."); + }; + return e2m(f2e(f, ecc), ecc); +} // gudermannian to true (only hyperbolas) e>1 (returns in range [-pi,pi]) inline double zeta2f(double f, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "zeta2f: zeta (Gudermannian) is not defined for ecc <=1."); + }; return 2 * std::atan(std::sqrt((1 + ecc) / (ecc - 1)) * std::tan(f / 2)); } // true to gudermannian (only hyperbolas) e>1 (returns in range [-pi,pi]) inline double f2zeta(double zeta, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "f2zeta: zeta (Gudermannian) is not defined for ecc <=1."); + }; return 2 * std::atan(std::sqrt((ecc - 1) / (1 + ecc)) * std::tan(zeta / 2)); } // mean to hyperbolic (only hyperbolas) e>1. // NOLINTNEXTLINE(bugprone-easily-swappable-parameters) inline double n2h(double N, double ecc) { - + if (ecc <= 1) { + throw std::domain_error( + "n2h: hyperbolic anomaly is not defined for ecc <=1."); + }; // The Initial guess (TODO(darioizo) improve) double IG = 1.; @@ -110,23 +148,49 @@ inline double n2h(double N, double ecc) { } // hyperbolic H to hyperbolic mean N (only hyperbolas) e>1 -inline double h2n(double H, double ecc) { return (ecc * std::sinh(H) - H); } +inline double h2n(double H, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "h2n: hyperbolic anomaly is not defined for ecc <=1."); + }; + return (ecc * std::sinh(H) - H); +} // hyperbolic H to true (only hyperbolas) e>1 (returns in range [-pi,pi]) inline double h2f(double H, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "h2f: hyperbolic anomaly is not defined for ecc <=1."); + }; return 2 * std::atan(std::sqrt((1 + ecc) / (ecc - 1)) * std::tanh(H / 2)); } // true to hyperbolic H (only hyperbolas) e>1 (returns in range [-pi,pi]) inline double f2h(double f, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "f2h: hyperbolic anomaly is not defined for ecc <=1."); + }; return 2 * std::atanh(std::sqrt((ecc - 1) / (1 + ecc)) * std::tan(f / 2)); } // mean hyperbolic to true (only hyperbolas) e>1 (returns in range [-pi,pi]) -inline double n2f(double N, double ecc) { return h2f(n2h(N, ecc), ecc); } +inline double n2f(double N, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "n2f: mean hyperbolic anomaly is not defined for ecc <=1."); + }; + return h2f(n2h(N, ecc), ecc); +} // true to mean hyperbolic (only hyperbolas) e>1 (returns in range [-pi,pi]) -inline double f2n(double f, double ecc) { return h2n(f2h(f, ecc), ecc); } +inline double f2n(double f, double ecc) { + if (ecc <= 1) { + throw std::domain_error( + "f2n: mean hyperbolic anomaly is not defined for ecc <=1."); + }; + return h2n(f2h(f, ecc), ecc); +} } // namespace kep3 #endif // kep3_TOOLBOX_M2E_H diff --git a/include/kep3/planets/jpl_lp.hpp b/include/kep3/planets/jpl_lp.hpp index 78e92e38..dcb95e73 100644 --- a/include/kep3/planets/jpl_lp.hpp +++ b/include/kep3/planets/jpl_lp.hpp @@ -54,7 +54,7 @@ class kep3_DLL_PUBLIC jpl_lp { public: // Constructor - explicit jpl_lp(const std::string & = "earth"); + explicit jpl_lp(std::string = "earth"); // Mandatory UDPLA methods [[nodiscard]] std::array, 2> eph(const epoch &) const; diff --git a/kep3-config.cmake.in b/kep3-config.cmake.in new file mode 100644 index 00000000..3a21bef0 --- /dev/null +++ b/kep3-config.cmake.in @@ -0,0 +1,14 @@ +# Mandatory public dependencies on Boost and fmt. +find_package(Boost @_kep3_MIN_BOOST_VERSION@ REQUIRED serialization) +find_package(fmt REQUIRED CONFIG) +find_package(heyoka REQUIRED CONFIG) +find_package(xtensor REQUIRED CONFIG) +find_package(xtensor-blas REQUIRED CONFIG) + +# Get current dir. +get_filename_component(_kep3_CONFIG_SELF_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH) + +include(${_kep3_CONFIG_SELF_DIR}/kep3_export.cmake) + +# Clean up. +unset(_kep3_CONFIG_SELF_DIR) \ No newline at end of file diff --git a/kep3_devel.yml b/kep3_devel.yml index c7871bf5..d337fea0 100644 --- a/kep3_devel.yml +++ b/kep3_devel.yml @@ -12,4 +12,5 @@ dependencies: - spdlog - xtensor - xtensor-blas + - pybind11 diff --git a/pykep/CMakeLists.txt b/pykep/CMakeLists.txt new file mode 100644 index 00000000..369f644b --- /dev/null +++ b/pykep/CMakeLists.txt @@ -0,0 +1,61 @@ +# Configure the version file. +configure_file("${CMAKE_CURRENT_SOURCE_DIR}/_version.py.in" "${CMAKE_CURRENT_SOURCE_DIR}/_version.py" @ONLY) + +# The list of pydsyre's Python files. +set(kep3_PYTHON_FILES __init__.py _version.py test.py) + +# Copy the python files in the current binary dir, +# so that we can import pydesyre from the build dir. +# NOTE: importing from the build dir will work +# only on single-configuration generators. +foreach(kep3_PYTHON_FILE ${kep3_PYTHON_FILES}) + configure_file("${CMAKE_CURRENT_SOURCE_DIR}/${kep3_PYTHON_FILE}" + "${CMAKE_CURRENT_BINARY_DIR}/${kep3_PYTHON_FILE}" COPYONLY) +endforeach() + +# Core module. +Python3_add_library(core MODULE WITH_SOABI + core.cpp + docstrings.cpp +) + +target_link_libraries(core PRIVATE kep3::kep3) +target_link_libraries(core PRIVATE "${pybind11_LIBRARIES}") + +target_include_directories(core SYSTEM PRIVATE "${pybind11_INCLUDE_DIR}" "${Python3_INCLUDE_DIRS}") +target_compile_definitions(core PRIVATE "${pybind11_DEFINITIONS}") +target_compile_options(core PRIVATE + "$<$:${KEP3_CXX_FLAGS_DEBUG}>" + "$<$:${KEP3_CXX_FLAGS_RELEASE}>" + "$<$:${KEP3_CXX_FLAGS_RELEASE}>" + "$<$:${KEP3_CXX_FLAGS_RELEASE}>" + +) +target_include_directories(core PUBLIC + $ + $ + $) +set_target_properties(core PROPERTIES CXX_VISIBILITY_PRESET hidden) +set_target_properties(core PROPERTIES VISIBILITY_INLINES_HIDDEN TRUE) +target_compile_features(core PRIVATE cxx_std_17) +set_property(TARGET core PROPERTY CXX_EXTENSIONS NO) + +# Installation setup. +if(PYKEP_INSTALL_PATH STREQUAL "") + message(STATUS "pykep will be installed in the default location: ${Python3_SITEARCH}") + set(_PYKEP_INSTALL_DIR "${Python3_SITEARCH}/pykep") +else() + message(STATUS "pykep will be installed in the custom location: ${PYKEP_INSTALL_PATH}") + set(_PYKEP_INSTALL_DIR "${PYKEP_INSTALL_PATH}/pykep") +endif() + +# Install the core module. +install(TARGETS core + RUNTIME DESTINATION ${_PYKEP_INSTALL_DIR} + LIBRARY DESTINATION ${_PYKEP_INSTALL_DIR} +) + +# Install the Python files. +install(FILES ${PYKEP_PYTHON_FILES} DESTINATION ${_PYKEP_INSTALL_DIR}) + +unset(_PYKEP_INSTALL_DIR) \ No newline at end of file diff --git a/pykep/__init__.py b/pykep/__init__.py new file mode 100644 index 00000000..a7967db6 --- /dev/null +++ b/pykep/__init__.py @@ -0,0 +1,18 @@ +"""Pykep is a coolbox for interplanetary trajectory design developed by ESA's Advanced Concepts Team. Its main + purpose is fast prototyping of reseacrh ideas, and is not intended for operational usage. + + Some important conventions followed: + 1 - All units expected are S.I. (m,sec,kg,N) unless explicitly stated. + 2 - The default set of osculating orbital parameters is, in this order: [sma, ecc, incl, W, w, f], where f is the true anomaly + 3 - The default option to represent epochs as floats is the modified julian date 2000 (MJD2000).""" + +# Version setup. +from ._version import __version__ +del _version + + +# Importing cpp functionalities +from .core import * + +# We import the unit test submodule +from . import test \ No newline at end of file diff --git a/pykep/_version.py b/pykep/_version.py new file mode 100644 index 00000000..f102a9ca --- /dev/null +++ b/pykep/_version.py @@ -0,0 +1 @@ +__version__ = "0.0.1" diff --git a/pykep/_version.py.in b/pykep/_version.py.in new file mode 100644 index 00000000..ecc1e1a0 --- /dev/null +++ b/pykep/_version.py.in @@ -0,0 +1 @@ +__version__ = "@kep3_VERSION@" \ No newline at end of file diff --git a/pykep/core.cpp b/pykep/core.cpp new file mode 100644 index 00000000..5b7e7132 --- /dev/null +++ b/pykep/core.cpp @@ -0,0 +1,75 @@ +// 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 +#include +#include + +#include "docstrings.hpp" + +namespace py = pybind11; +namespace pk = pykep; + +PYBIND11_MODULE(core, m) { + py::options options; + options.disable_function_signatures(); + m.doc() = pk::core_module_doc(); + + // We expose various global constants: + m.attr("AU") = py::float_(kep3::AU); + m.attr("CAVENDISH") = py::float_(kep3::CAVENDISH); + m.attr("MU_SUN") = py::float_(kep3::MU_SUN); + m.attr("MU_EARTH") = py::float_(kep3::MU_EARTH); + m.attr("EARTH_VELOCITY") = py::float_(kep3::EARTH_VELOCITY); + m.attr("EARTH_J2") = py::float_(kep3::EARTH_J2); + m.attr("EARTH_RADIUS") = py::float_(kep3::EARTH_RADIUS); + m.attr("RAD2DEG") = py::float_(kep3::RAD2DEG); + m.attr("DAY2SEC") = py::float_(kep3::DAY2SEC); + m.attr("SEC2DAY") = py::float_(kep3::SEC2DAY); + m.attr("DAY2YEAR") = py::float_(kep3::DAY2YEAR); + m.attr("G0") = py::float_(kep3::G0); + + // We expose here global enums: + py::enum_(m, "elements_type", "") + .value("KEP_M", kep3::KEP_M, + "Keplerian Elements a,e,i,W,w,M (Mean anomaly)") + .value("KEP_F", kep3::KEP_F, + "Keplerian Elements a,e,i,W,w,f (True anomaly)") + .value("MEQ", kep3::MEQ, + "Modified Equinoctial Elements p,f,g,h,k,L (Mean Longitude)") + .value("MEQ_R", kep3::MEQ_R, + "Modified Equinoctial Elements (retrograde) p,f,g,h,k,L (Mean " + "Longitude)") + .value("POSVEL", kep3::POSVEL, "Position and Velocity") + .export_values(); + + // We expose the various anomaly conversions + m.def("m2e", &kep3::m2e, pk::m2e_doc().c_str()); + m.def("e2m", &kep3::e2m, pk::e2m_doc().c_str()); + m.def("m2f", &kep3::m2f, pk::m2f_doc().c_str()); + m.def("f2m", &kep3::f2m, pk::f2m_doc().c_str()); + m.def("e2f", &kep3::e2f, pk::e2f_doc().c_str()); + m.def("f2e", &kep3::f2e, pk::f2e_doc().c_str()); + m.def("n2h", &kep3::n2h, pk::n2h_doc().c_str()); + m.def("h2n", &kep3::h2n, pk::h2n_doc().c_str()); + m.def("n2f", &kep3::n2f, pk::n2f_doc().c_str()); + m.def("f2n", &kep3::f2n, pk::f2n_doc().c_str()); + m.def("h2f", &kep3::h2f, pk::h2f_doc().c_str()); + m.def("f2h", &kep3::f2h, pk::f2h_doc().c_str()); + m.def("zeta2f", &kep3::zeta2f, pk::zeta2f_doc().c_str()); + m.def("f2zeta", &kep3::f2zeta, pk::f2zeta_doc().c_str()); +} \ No newline at end of file diff --git a/pykep/docstrings.cpp b/pykep/docstrings.cpp new file mode 100644 index 00000000..edf76599 --- /dev/null +++ b/pykep/docstrings.cpp @@ -0,0 +1,317 @@ +// 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 "docstrings.hpp" + +namespace pykep { + +std::string core_module_doc() { + return R"(core is the Pykep module that contains most of its core routines efficiently coded in c++ +)"; +} + +std::string m2e_doc() { + return R"(m2e(M, ecc) + + Converts from Mean to Eccentric anomaly. Requires ecc < 1. + + Args: + M (float): the Mean anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Eccentric anomaly in [-pi, pi] (rad.) + + Examples: + >>> import pykep as pk + >>> M = 1.2 + >>> ecc = 0.1 + >>> pk.m2e(M, ecc) + 1.1929459841174401 +)"; +} + +std::string e2m_doc() { + return R"(e2m(E, ecc) + + Converts from Eccentric to Mean anomaly. Requires ecc < 1. + + Args: + E (float): the Eccentric anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Mean anomaly (rad.) + + Examples: + >>> import pykep as pk + >>> E = 0.5 + >>> ecc = 0.1 + >>> pk.e2m(E, ecc) + 0.4520574461395797 +)"; +} + +std::string e2f_doc() { + return R"(e2f(E, ecc) + + Converts from eccentric to true anomaly. Requires ecc < 1. + + Args: + E (float): the Eccentric anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the True anomaly in [-pi, pi] (rad.) + + Examples: + >>> import pykep as pk + >>> E = 0.5 + >>> ecc = 0.1 + >>> pk.e2f(E, ecc) + 0.5502639747136633 +)"; +} + +std::string f2e_doc() { + return R"(f2e(f, ecc) + + Converts from True to Eccentric anomaly. Requires ecc < 1. + + Args: + f (float): the True anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Eccentric anomaly in [-pi, pi] (rad.) + + Examples: + >>> import pykep as pk + >>> f = 1.2 + >>> ecc = 0.1 + >>> pk.f2e(f, ecc) + 1.1082931139529482 +)"; +} + +std::string f2m_doc() { + return R"(f2m(f, ecc) + + Converts from True to Mean anomaly. Requires ecc < 1. + + Args: + f (float): the True anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Mean anomaly in [-pi, pi] (rad.) + + Examples: + >>> import pykep as pk + >>> f = -0.34 + >>> ecc = 0.67 + >>> pk.f2m(f, ecc) + -0.8380280766377411 +)"; +} + +std::string m2f_doc() { + return R"(m2f(M, ecc) + + Converts from Mean to True anomaly. Requires ecc < 1. + + Args: + M (float): the Mean anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the True anomaly in [-pi, pi] (rad.) + + Examples: + >>> import pykep as pk + >>> M = 0.32 + >>> ecc = 0.65 + >>> pk.m2f(M, ecc) + 1.4497431281728277 +)"; +} + +std::string h2n_doc() { + return R"(h2n(H, ecc) + + Converts from Hyperbolic to Hyperbolic Mean anomaly. Requires ecc > 1. + + Args: + H (float): the Hyperbolic anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Hyperbolic Mean anomaly (rad.) + + Examples: + >>> import pykep as pk + >>> H = 1.2 + >>> ecc = 10.32 + >>> pk.h2n(H, ecc) + 14.377641187853621 +)"; +} + +std::string n2h_doc() { + return R"(n2h(N, ecc) + + Converts from Hyperbolic Mean to Hyperbolic anomaly. Requires ecc > 1. + + Args: + N (float): the Hyperbolic Mean anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Hyperbolic anomaly (rad.) + + Examples: + >>> import pykep as pk + >>> N = 1.2 + >>> ecc = 10.32 + >>> pk.n2h(N, ecc) + 0.12836469743916526 +)"; +} + +std::string h2f_doc() { + return R"(h2f(H, ecc) + + Converts from Hyperbolic to True anomaly. Requires ecc > 1. + + Args: + H (float): the Hyperbolic anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the True anomaly in [-pi, pi] (rad.) + + Examples: + >>> import pykep as pk + >>> H = 10.32 + >>> ecc = 4.5 + >>> pk.h2f(H, ecc) + 1.7948251330114304 +)"; +} + +std::string f2h_doc() { + return R"(f2h(f, ecc) + + Converts from True to Hyperbolic anomaly. Requires ecc > 1. + + Args: + f (float): the True anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Hyperbolic anomaly + + Examples: + >>> import pykep as pk + >>> f = 1.2 + >>> ecc = 0.1 + >>> pk.f2h(f, ecc) + 1.1929459841174401 +)"; +} + +std::string f2n_doc() { + return R"(f2n(f, ecc) + + Converts from True to Hyperbolic Mean anomaly. Requires ecc > 1. + + Args: + f (float): the True anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Hyperbolic Mean anomaly + + Examples: + >>> import pykep as pk + >>> f = 1.2 + >>> ecc = 5.7 + >>> pk.f2n(f, ecc) + 8.421335633880908 +)"; +} + +std::string n2f_doc() { + return R"(n2f(N, ecc) + + Converts from Hyperbolic Mean to True anomaly. Requires ecc > 1. + + Args: + N (float): the Hyperbolic Mean anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the True anomaly + + Examples: + >>> import pykep as pk + >>> N = 10.32 + >>> ecc = 13.45 + >>> pk.m2e(M, ecc) + 0.7373697968359353 +)"; +} + +std::string zeta2f_doc() { + return R"(zeta2f(zeta, ecc) + + Converts from Gudermannian to True anomaly. Requires ecc > 1. + + See Battin: "An Introduction to the Mathematics and Methods of Astrodynamics" for a + definition of zeta and the treatment of the resulting equations. + + Args: + zeta (float): the Gudermannian (rad.) + ecc (float): the eccentricity + + Returns: + float: the True anomaly + + Examples: + >>> import pykep as pk + >>> zeta = 8.2 + >>> ecc = 2.2 + >>> pk.zeta2f(zeta, ecc) + 2.3290929552114266 +)"; +} + +std::string f2zeta_doc() { + return R"(f2zeta(f, ecc) + + Converts from True anomaly to Gudermannian. Requires ecc > 1. + + Args: + f (float): the True anomaly (rad.) + ecc (float): the eccentricity + + Returns: + float: the Gudermannian + + Examples: + >>> import pykep as pk + >>> f = 0.5 + >>> ecc = 3.3 + >>> pk.f2zeta(f, ecc) + 0.36923933496389816 +)"; +} +} // namespace pykep \ No newline at end of file diff --git a/pykep/docstrings.hpp b/pykep/docstrings.hpp new file mode 100644 index 00000000..8a56c11d --- /dev/null +++ b/pykep/docstrings.hpp @@ -0,0 +1,38 @@ +// 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/. + +#ifndef PYKEP_DOCSTRINGS_HPP +#define PYKEP_DOCSTRINGS_HPP + +#include + +namespace pykep { +// Modules +std::string core_module_doc(); +// Anomaly conversions +std::string m2e_doc(); +std::string e2m_doc(); +std::string m2f_doc(); +std::string f2m_doc(); +std::string e2f_doc(); +std::string f2e_doc(); + +std::string n2h_doc(); +std::string h2n_doc(); +std::string n2f_doc(); +std::string f2n_doc(); +std::string h2f_doc(); +std::string f2h_doc(); + +std::string zeta2f_doc(); +std::string f2zeta_doc(); + +} // namespace pykep + +#endif \ No newline at end of file diff --git a/pykep/test.py b/pykep/test.py new file mode 100644 index 00000000..4b2f572f --- /dev/null +++ b/pykep/test.py @@ -0,0 +1,63 @@ +## Copyright 2023, 2024 Dario Izzo (dario.izzo@gmail.com), Francesco Biscani +## (bluescarni@gmail.com) +## +## This file is part of the pykep 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/. + +import unittest as _ut + + +def float_abs_error(a: float, b: float): + return abs(a - b) + +class anomaly_conversions_tests(_ut.TestCase): + def test_m2e(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.m2e(pk.e2m(0.1, 0.1), 0.1), 0.1) < 1e-14) + + def test_m2f(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.m2f(pk.f2m(0.1, 0.1), 0.1), 0.1) < 1e-14) + + def test_f2e(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.f2e(pk.e2f(0.1, 0.1), 0.1), 0.1) < 1e-14) + + def test_n2h(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.n2h(pk.h2n(0.1, 10.1), 10.1), 0.1) < 1e-14) + + def test_n2f(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.n2f(pk.f2n(0.1, 10.1), 10.1), 0.1) < 1e-14) + + def test_f2h(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.f2h(pk.h2f(0.1, 10.1), 10.1), 0.1) < 1e-14) + + def test_f2zeta(self): + import pykep as pk + + self.assertTrue(float_abs_error(pk.f2zeta(pk.zeta2f(0.1, 10.1), 10.1), 0.1) < 1e-14) + + +def run_test_suite(): + suite = _ut.TestSuite() + suite.addTest(anomaly_conversions_tests("test_m2e")) + suite.addTest(anomaly_conversions_tests("test_m2f")) + suite.addTest(anomaly_conversions_tests("test_f2e")) + suite.addTest(anomaly_conversions_tests("test_n2h")) + suite.addTest(anomaly_conversions_tests("test_n2f")) + suite.addTest(anomaly_conversions_tests("test_f2h")) + suite.addTest(anomaly_conversions_tests("test_f2zeta")) + + test_result = _ut.TextTestRunner(verbosity=2).run(suite) diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index f4f0bd2c..1ba83b31 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -7,8 +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/. -#include "kep3/core_astro/eq2par2eq.hpp" -#include "kep3/epoch.hpp" #include #include #include @@ -19,6 +17,8 @@ #include #include +#include +#include #include #include @@ -45,18 +45,17 @@ static const std::array neptune_el = {30.06992276, 0.00859048, 1.7700 static const std::array neptune_el_dot = {0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664}; // clang-format on -jpl_lp::jpl_lp(const std::string &name) - : m_elements(), m_elements_dot(), m_name(name), - m_mu_central_body(kep3::MU_SUN) { +// NOLINTNEXTLINE(cert-err58-cpp) +const std::unordered_map mapped_planets = { + {"mercury", 1}, {"venus", 2}, {"earth", 3}, {"mars", 4}, {"jupiter", 5}, + {"saturn", 6}, {"uranus", 7}, {"neptune", 8}, {"pluto", 9}}; - // Cannot be marked const as operator[] is not and we use it. - static std::unordered_map mapped_planets = { - {"mercury", 1}, {"venus", 2}, {"earth", 3}, - {"mars", 4}, {"jupiter", 5}, {"saturn", 6}, - {"uranus", 7}, {"neptune", 8}, {"pluto", 9}}; +jpl_lp::jpl_lp(std::string name) + : m_elements(), m_elements_dot(), m_name(std::move(name)), + m_mu_central_body(kep3::MU_SUN) { boost::algorithm::to_lower(m_name); - switch (mapped_planets[m_name]) { + switch (mapped_planets.at(m_name)) { case (1): { m_elements = mercury_el; m_elements_dot = mercury_el_dot; @@ -114,7 +113,7 @@ jpl_lp::jpl_lp(const std::string &name) m_mu_self = 6836529e9; } break; default: { - throw std::logic_error("unknown planet name: "); + throw std::logic_error("unknown planet name: "); // LCOV_EXCL_LINE } } } diff --git a/test/convert_anomalies_test.cpp b/test/convert_anomalies_test.cpp index 9f4e4d88..9ebcc759 100644 --- a/test/convert_anomalies_test.cpp +++ b/test/convert_anomalies_test.cpp @@ -11,6 +11,7 @@ #include #include +#include #include "catch.hpp" @@ -133,3 +134,22 @@ TEST_CASE("zeta2e") { Approx(std::cos(true_anom)).epsilon(0.).margin(1e-14)); } } + +TEST_CASE("throws") { + REQUIRE_THROWS_AS(kep3::m2e(0.3,1.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::e2m(0.3,1.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::m2f(0.3,1.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::f2m(0.3,1.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::e2f(0.3,1.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::f2e(0.3,1.1), std::domain_error); + + REQUIRE_THROWS_AS(kep3::n2h(0.3,0.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::h2n(0.3,0.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::n2f(0.3,0.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::f2n(0.3,0.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::h2f(0.3,0.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::f2h(0.3,0.1), std::domain_error); + + REQUIRE_THROWS_AS(kep3::zeta2f(0.3,0.1), std::domain_error); + REQUIRE_THROWS_AS(kep3::f2zeta(0.3,0.1), std::domain_error); +} \ No newline at end of file diff --git a/test/lambert_problem_test.cpp b/test/lambert_problem_test.cpp index 53d7b6db..11286ac6 100644 --- a/test/lambert_problem_test.cpp +++ b/test/lambert_problem_test.cpp @@ -66,7 +66,6 @@ TEST_CASE("delta_guidance") { tof = tof_d(rng_engine); bool cw = static_cast(cw_d(rng_engine)); double mu = mu_d(rng_engine); - // 2 - Solve the lambert problem kep3::lambert_problem lp(r1, r2, tof, mu, cw, revs_max); diff --git a/tools/gha_linux_pykep.sh b/tools/gha_linux_pykep.sh new file mode 100644 index 00000000..d0a922e7 --- /dev/null +++ b/tools/gha_linux_pykep.sh @@ -0,0 +1,35 @@ +#!/usr/bin/env bash + +# Echo each command +set -x + +# Exit on error. +set -e + +# Core deps. +sudo apt-get install wget + +# Install conda+deps. +wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-Linux-x86_64.sh -O mambaforge.sh +export deps_dir=$HOME/local +export PATH="$HOME/mambaforge/bin:$PATH" +bash mambaforge.sh -b -p $HOME/mambaforge +mamba env create -f kep3_devel.yml -q -p $deps_dir +source activate $deps_dir + +# First we build and install the kep3 library +mkdir build +cd build +cmake -G "Ninja" ../ -DCMAKE_INSTALL_PREFIX=$deps_dir -DCMAKE_PREFIX_PATH=$deps_dir -DCMAKE_BUILD_TYPE=Release -Dkep3_BUILD_TESTS=no -Dkep3_BUILD_BENCHMARKS=no -Dkep3_BUILD_PYTHON_BINDINGS=no -DBoost_NO_BOOST_CMAKE=ON +cmake --build . --target=install --config=Release -- -j 2 +# Then we build and install pykep +cd .. +mkdir build_pykep +cd build_pykep +cmake -G "Ninja" ../ -DCMAKE_INSTALL_PREFIX=$deps_dir -DCMAKE_PREFIX_PATH=$deps_dir -DCMAKE_BUILD_TYPE=Release -Dkep3_BUILD_TESTS=no -Dkep3_BUILD_BENCHMARKS=no -Dkep3_BUILD_PYTHON_BINDINGS=yes -DBoost_NO_BOOST_CMAKE=ON +cmake --build . --target=install --config=Release -- -j 2 + +python -c "import pykep.test; pykep.test.run_test_suite()" + +set +e +set +x