Merge pull request #352 from bluescarni/pr/kepE_improvements

kepF

CMakeLists.txt

@@ -11,7 +11,7 @@ if(NOT CMAKE_BUILD_TYPE)
 	FORCE)
 endif()
 
-project(heyoka VERSION 3.0.0 LANGUAGES CXX C)
+project(heyoka VERSION 3.1.0 LANGUAGES CXX C)
 
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/yacma")
 
@@ -263,6 +263,8 @@ set(HEYOKA_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/cr3bp.cpp"
     # Math functions.
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepE.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepF.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepDE.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/cos.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/exp.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/log.cpp"
@@ -302,7 +304,7 @@ if(HEYOKA_WITH_SLEEF)
 endif()
 
 # Setup the heyoka ABI version number.
-set(HEYOKA_ABI_VERSION 24)
+set(HEYOKA_ABI_VERSION 25)
 
 if(HEYOKA_BUILD_STATIC_LIBRARY)
     # Setup of the heyoka static library.

doc/ad_notes.rst

@@ -588,10 +588,12 @@ and :eq:`eq_leibniz_00` to obtain, for :math:`n > 0`:
 Celestial mechanics
 -------------------
 
-Inverse of Kepler's elliptic equation
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. _kepE_ad:
 
-The inverse of Kepler's elliptic equation is the bivariate function :math:`E = E\left( e, M \right)` implicitly defined by the
+Kepler's eccentric anomaly
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The `eccentric anomaly <https://en.wikipedia.org/wiki/Eccentric_anomaly>`__ is the bivariate function :math:`E = E\left( e, M \right)` implicitly defined by the
 trascendental equation
 
 .. math::
@@ -661,6 +663,59 @@ and :eq:`eq_leibniz_00` and re-arrange to obtain, for :math:`n > 0`:
    \right)
    \right].
 
+Eccentric longitude
+^^^^^^^^^^^^^^^^^^^
+
+The `eccentric longitude <https://articles.adsabs.harvard.edu//full/1972CeMec...5..303B/0000309.000.html>`__
+is the trivariate function :math:`F = F\left( h, k, \lambda \right)` implicitly defined by the
+trascendental equation
+
+.. math::
+   :label:
+
+   \lambda = F + h\cos F - k\sin F,
+
+with :math:`h^2+k^2 < 1`. Given :math:`a\left( t \right) = F\left( h\left( t \right), k \left( t \right), \lambda \left( t \right) \right)`,
+we have
+
+.. math::
+   :label:
+
+   a^\prime\left( t \right) = \frac{k^\prime \left( t \right)\sin a\left(t \right) - h^\prime \left( t \right)\cos a\left(t \right)
+   + \lambda^\prime\left( t \right)}{1-h\left( t \right)\sin a\left(t \right) -k\left( t \right)\cos a\left(t \right)}.
+
+After the introduction of the auxiliary functions
+
+.. math::
+   :label:
+
+   \begin{cases}
+   c\left( t \right) = h\left( t \right) \sin a\left(t \right), \\
+   d\left( t \right) = k\left( t \right) \cos a\left(t \right), \\
+   e\left( t \right) = \sin a\left(t \right), \\
+   f\left( t \right) = \cos a\left(t \right), \\
+   \end{cases}
+
+we can then proceed in the same way as explained for the :ref:`eccentric anomaly <kepE_ad>`.
+The final result, for :math:`n > 0`, is:
+
+.. math::
+   :label:
+
+   a^{\left[ n \right]}\left( t \right) = \frac{1}{n \left(1 - c^{\left[ 0 \right]}\left( t \right) -d^{\left[ 0 \right]}\left( t \right)\right)}
+   \left\{
+   n\left( k^{\left[ n \right]}\left( t \right) e^{\left[ 0 \right]}\left( t \right)
+   - h^{\left[ n \right]}\left( t \right) f^{\left[ 0 \right]}\left( t \right)
+   + \lambda^{\left[ n \right]}\left( t \right)\right) +
+   \sum_{j=1}^{n-1}j\left[
+      a^{\left[ j \right]}\left( t \right) \left(
+         c^{\left[ n - j \right]}\left( t \right) + d^{\left[ n - j \right]}\left( t \right)
+      \right)
+      + k^{\left[ j \right]}\left( t \right) e^{\left[ n - j \right]}\left( t \right)
+      - h^{\left[ j \right]}\left( t \right) f^{\left[ n - j \right]}\left( t \right)
+   \right]
+   \right\}.
+
 Time functions
 --------------
 

doc/changelog.rst

@@ -1,6 +1,24 @@
 Changelog
 =========
 
+3.1.0 (unreleased)
+------------------
+
+New
+~~~
+
+- Implement the eccentric longitude :math:`F` in the expression
+  system (`#352 <https://github.com/bluescarni/heyoka/pull/352>`__).
+- Implement the delta eccentric anomaly :math:`\Delta E` in the expression
+  system (`#352 <https://github.com/bluescarni/heyoka/pull/352>`__).
+  Taylor derivatives are not implemented yet.
+
+Fix
+~~~
+
+- Improve the numerical stability of the Kepler solver
+  (`#352 <https://github.com/bluescarni/heyoka/pull/352>`__).
+
 3.0.0 (2023-10-07)
 ------------------
 

doc/index.rst

@@ -120,5 +120,6 @@ Dario Izzo (European Space Agency).
    ad_notes.rst
    changelog.rst
    breaking_changes.rst
+   known_issues.rst
    acknowledgement.rst
    bibliography.rst

doc/known_issues.rst

@@ -0,0 +1,8 @@
+Known issues
+============
+
+* Due to an upstream bug, if you compile heyoka linking statically against LLVM 17
+  while enabling the ``HEYOKA_HIDE_LLVM_SYMBOLS`` option (see the
+  :ref:`installation instructions <installation_from_source>`), you may experience
+  runtime errors. This problem should be fixed in LLVM 18. A patch fixing the issue in LLVM 17
+  is available `here <https://github.com/llvm/llvm-project/commit/122ebe3b500190b1f408e2e6db753853e297ba28>`__.

doc/variable.rst

@@ -20,7 +20,7 @@ The :cpp:class:`~heyoka::variable` class
       Variable names beginning with a double underscore (``__``) are reserved
       for internal use.
 
-   .. cpp:function:: variable()
+   .. cpp:function:: variable() noexcept
 
       Default constructor.
 

include/heyoka/detail/llvm_helpers.hpp

@@ -123,6 +123,7 @@ HEYOKA_DLL_PUBLIC llvm::Value *llvm_square(llvm_state &, llvm::Value *);
 HEYOKA_DLL_PUBLIC llvm::Constant *llvm_constantfp(llvm_state &, llvm::Type *, double);
 
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_fcmp_ult(llvm_state &, llvm::Value *, llvm::Value *);
+HEYOKA_DLL_PUBLIC llvm::Value *llvm_fcmp_uge(llvm_state &, llvm::Value *, llvm::Value *);
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_fcmp_oge(llvm_state &, llvm::Value *, llvm::Value *);
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_fcmp_ole(llvm_state &, llvm::Value *, llvm::Value *);
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_fcmp_olt(llvm_state &, llvm::Value *, llvm::Value *);
@@ -170,9 +171,11 @@ HEYOKA_DLL_PUBLIC std::pair<llvm::Value *, llvm::Value *>
 llvm_penc_cargo_shisha(llvm_state &, llvm::Type *, llvm::Value *, std::uint32_t, llvm::Value *, std::uint32_t);
 
 llvm::Function *llvm_add_inv_kep_E(llvm_state &, llvm::Type *, std::uint32_t);
-
 HEYOKA_DLL_PUBLIC void llvm_add_inv_kep_E_wrapper(llvm_state &, llvm::Type *, std::uint32_t, const std::string &);
 
+llvm::Function *llvm_add_inv_kep_F(llvm_state &, llvm::Type *, std::uint32_t);
+llvm::Function *llvm_add_inv_kep_DE(llvm_state &, llvm::Type *, std::uint32_t);
+
 HEYOKA_DLL_PUBLIC llvm::Value *llvm_add_bc_array(llvm_state &, llvm::Type *, std::uint32_t);
 
 // Primitives for double-length arithmetic.

include/heyoka/detail/real_helpers.hpp

@@ -34,6 +34,7 @@ llvm::Value *llvm_real_fneg(llvm_state &, llvm ::Value *);
 llvm::Function *real_nary_op(llvm_state &, llvm::Type *, const std::string &, unsigned);
 std::pair<llvm::Value *, llvm::Value *> llvm_real_sincos(llvm_state &, llvm::Value *);
 llvm::Value *llvm_real_fcmp_ult(llvm_state &, llvm::Value *, llvm::Value *);
+llvm::Value *llvm_real_fcmp_uge(llvm_state &, llvm::Value *, llvm::Value *);
 llvm::Value *llvm_real_fcmp_oge(llvm_state &, llvm::Value *, llvm::Value *);
 llvm::Value *llvm_real_fcmp_ole(llvm_state &, llvm::Value *, llvm::Value *);
 llvm::Value *llvm_real_fcmp_olt(llvm_state &, llvm::Value *, llvm::Value *);

include/heyoka/math.hpp

@@ -21,7 +21,9 @@
 #include <heyoka/math/cosh.hpp>
 #include <heyoka/math/erf.hpp>
 #include <heyoka/math/exp.hpp>
+#include <heyoka/math/kepDE.hpp>
 #include <heyoka/math/kepE.hpp>
+#include <heyoka/math/kepF.hpp>
 #include <heyoka/math/log.hpp>
 #include <heyoka/math/pow.hpp>
 #include <heyoka/math/prod.hpp>

include/heyoka/math/kepDE.hpp

@@ -0,0 +1,96 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani ([email protected]), Dario Izzo ([email protected])
+//
+// This file is part of the heyoka 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 HEYOKA_MATH_KEPDE_HPP
+#define HEYOKA_MATH_KEPDE_HPP
+
+#include <heyoka/config.hpp>
+
+#include <cstdint>
+#include <vector>
+
+#if defined(HEYOKA_HAVE_REAL128)
+
+#include <mp++/real128.hpp>
+
+#endif
+
+#if defined(HEYOKA_HAVE_REAL)
+
+#include <mp++/real.hpp>
+
+#endif
+
+#include <heyoka/detail/func_cache.hpp>
+#include <heyoka/detail/fwd_decl.hpp>
+#include <heyoka/detail/llvm_fwd.hpp>
+#include <heyoka/detail/visibility.hpp>
+#include <heyoka/func.hpp>
+#include <heyoka/s11n.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+namespace detail
+{
+
+class HEYOKA_DLL_PUBLIC kepDE_impl : public func_base
+{
+    friend class boost::serialization::access;
+    template <typename Archive>
+    HEYOKA_DLL_LOCAL void serialize(Archive &, unsigned);
+
+    template <typename T>
+    HEYOKA_DLL_LOCAL expression diff_impl(funcptr_map<expression> &, const T &) const;
+
+public:
+    kepDE_impl();
+    explicit kepDE_impl(expression, expression, expression);
+
+    expression diff(funcptr_map<expression> &, const std::string &) const;
+    expression diff(funcptr_map<expression> &, const param &) const;
+
+    [[nodiscard]] llvm::Value *llvm_eval(llvm_state &, llvm::Type *, const std::vector<llvm::Value *> &, llvm::Value *,
+                                         llvm::Value *, llvm::Value *, std::uint32_t, bool) const;
+
+    [[nodiscard]] llvm::Function *llvm_c_eval_func(llvm_state &, llvm::Type *, std::uint32_t, bool) const;
+};
+
+} // namespace detail
+
+HEYOKA_DLL_PUBLIC expression kepDE(expression, expression, expression);
+
+#define HEYOKA_DECLARE_KEPDE_OVERLOADS(type)                                                                           \
+    HEYOKA_DLL_PUBLIC expression kepDE(expression, type, type);                                                        \
+    HEYOKA_DLL_PUBLIC expression kepDE(type, expression, type);                                                        \
+    HEYOKA_DLL_PUBLIC expression kepDE(type, type, expression);                                                        \
+    HEYOKA_DLL_PUBLIC expression kepDE(expression, expression, type);                                                  \
+    HEYOKA_DLL_PUBLIC expression kepDE(expression, type, expression);                                                  \
+    HEYOKA_DLL_PUBLIC expression kepDE(type, expression, expression)
+
+HEYOKA_DECLARE_KEPDE_OVERLOADS(double);
+HEYOKA_DECLARE_KEPDE_OVERLOADS(long double);
+
+#if defined(HEYOKA_HAVE_REAL128)
+
+HEYOKA_DECLARE_KEPDE_OVERLOADS(mppp::real128);
+
+#endif
+
+#if defined(HEYOKA_HAVE_REAL)
+
+HEYOKA_DECLARE_KEPDE_OVERLOADS(mppp::real);
+
+#endif
+
+#undef HEYOKA_DECLARE_KEPDE_OVERLOADS
+
+HEYOKA_END_NAMESPACE
+
+HEYOKA_S11N_FUNC_EXPORT_KEY(heyoka::detail::kepDE_impl)
+
+#endif

include/heyoka/math/kepE.hpp

@@ -73,36 +73,27 @@ class HEYOKA_DLL_PUBLIC kepE_impl : public func_base
 
 HEYOKA_DLL_PUBLIC expression kepE(expression, expression);
 
-HEYOKA_DLL_PUBLIC expression kepE(expression, double);
-HEYOKA_DLL_PUBLIC expression kepE(expression, long double);
+#define HEYOKA_DECLARE_KEPE_OVERLOADS(type)                                                                            \
+    HEYOKA_DLL_PUBLIC expression kepE(expression, type);                                                               \
+    HEYOKA_DLL_PUBLIC expression kepE(type, expression);
 
-#if defined(HEYOKA_HAVE_REAL128)
-
-HEYOKA_DLL_PUBLIC expression kepE(expression, mppp::real128);
-
-#endif
-
-#if defined(HEYOKA_HAVE_REAL)
-
-HEYOKA_DLL_PUBLIC expression kepE(expression, mppp::real);
-
-#endif
-
-HEYOKA_DLL_PUBLIC expression kepE(double, expression);
-HEYOKA_DLL_PUBLIC expression kepE(long double, expression);
+HEYOKA_DECLARE_KEPE_OVERLOADS(double);
+HEYOKA_DECLARE_KEPE_OVERLOADS(long double);
 
 #if defined(HEYOKA_HAVE_REAL128)
 
-HEYOKA_DLL_PUBLIC expression kepE(mppp::real128, expression);
+HEYOKA_DECLARE_KEPE_OVERLOADS(mppp::real128);
 
 #endif
 
 #if defined(HEYOKA_HAVE_REAL)
 
-HEYOKA_DLL_PUBLIC expression kepE(mppp::real, expression);
+HEYOKA_DECLARE_KEPE_OVERLOADS(mppp::real);
 
 #endif
 
+#undef HEYOKA_DECLARE_KEPE_OVERLOADS
+
 HEYOKA_END_NAMESPACE
 
 HEYOKA_S11N_FUNC_EXPORT_KEY(heyoka::detail::kepE_impl)