Merge pull request #5 from darioizzo/lambert

Lambert

CMakeLists.txt

@@ -131,6 +131,7 @@ endif()
 set(kep3_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/epoch.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/planet.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/lambert_problem.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/planets/keplerian.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/ic2par2ic.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/ic2eq2ic.cpp"

benchmark/CMakeLists.txt

@@ -18,3 +18,5 @@ endfunction()
 
 ADD_kep3_BENCHMARK(convert_anomalies_benchmark)
 ADD_kep3_BENCHMARK(propagate_lagrangian_benchmark)
+ADD_kep3_BENCHMARK(lambert_problem_benchmark)
+

benchmark/lambert_problem_benchmark.cpp

@@ -0,0 +1,101 @@
+/*****************************************************************************
+ *   Copyright (C) 2004-2018 The pykep development team,                     *
+ *   Advanced Concepts Team (ACT), European Space Agency (ESA)               *
+ *                                                                           *
+ *   https://gitter.im/esa/pykep                                             *
+ *   https://github.com/esa/pykep                                            *
+ *                                                                           *
+ *   [email protected]                                                             *
+ *                                                                           *
+ *   This program is free software; you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation; either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program; if not, write to the                           *
+ *   Free Software Foundation, Inc.,                                         *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.               *
+ *****************************************************************************/
+
+#include <array>
+#include <chrono>
+#include <iomanip>
+#include <iostream>
+#include <random>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <kep3/core_astro/constants.hpp>
+#include <kep3/core_astro/propagate_lagrangian.hpp>
+#include <kep3/lambert_problem.hpp>
+#include <stdexcept>
+
+using std::chrono::duration_cast;
+using std::chrono::high_resolution_clock;
+using std::chrono::microseconds;
+
+int main() {
+  // Number of trials
+  const unsigned trials = 50000u;
+
+  std::array<std::array<double, 3>, trials> r1s{}, r2s{};
+  std::array<double, trials> tof{};
+  std::array<bool, trials> cw{};
+  std::array<double, trials> mu{};
+
+  // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
+  std::mt19937 rng_engine(122012203u);
+  std::uniform_int_distribution<unsigned> cw_d(0, 1);
+  std::uniform_real_distribution<double> r_d(-2, 2);
+  std::uniform_real_distribution<double> tof_d(2., 40.);
+  std::uniform_real_distribution<double> mu_d(0.9, 1.1);
+  unsigned revs_max = 20u;
+  unsigned long count = 0lu;
+
+  for (auto i = 0u; i < trials; ++i) {
+    // 1 - generate a random problem geometry
+    r1s[i][0] = r_d(rng_engine);
+    r1s[i][1] = r_d(rng_engine);
+    r1s[i][2] = r_d(rng_engine);
+    r2s[i][0] = r_d(rng_engine);
+    r2s[i][1] = r_d(rng_engine);
+    r2s[i][2] = r_d(rng_engine);
+    tof[i] = tof_d(rng_engine);
+    cw[i] = static_cast<bool>(cw_d(rng_engine));
+    mu[i] = mu_d(rng_engine);
+  }
+
+  auto start = high_resolution_clock::now();
+  for (auto i = 0u; i < trials; ++i) {
+    // 2 - Solve the lambert problem
+    kep3::lambert_problem lp(r1s[i], r2s[i], tof[i], mu[i], cw[i], revs_max);
+    count = count + lp.get_v1().size();
+  }
+  auto stop = high_resolution_clock::now();
+  auto duration = duration_cast<microseconds>(stop - start);
+  fmt::print("Lambert:\n{} solutions computed in {:.3f}s\n", count,
+             (static_cast<double>(duration.count()) / 1e6));
+  fmt::print("Projected number of solutions per second: {}\n",
+             static_cast<double>(count) / ((static_cast<double>(duration.count()) / 1e6)));
+
+  count = 0; //reset counter
+  start = high_resolution_clock::now();
+  for (auto i = 0u; i < trials; ++i) {
+    // 3 - Solve the lambert problem for the singe rev case
+    kep3::lambert_problem lp(r1s[i], r2s[i], tof[i], mu[i], cw[i], 0u);
+    count += 1u;
+  }
+  stop = high_resolution_clock::now();
+  duration = duration_cast<microseconds>(stop - start);
+  fmt::print("\nLambert (0 revs only):\n{} solutions computed in {:.3f}s\n", count,
+             (static_cast<double>(duration.count()) / 1e6));
+  fmt::print("Projected number of solutions per second: {}\n",
+             static_cast<double>(count) / ((static_cast<double>(duration.count()) / 1e6)));
+}

benchmark/propagate_lagrangian_benchmark.cpp

@@ -7,10 +7,10 @@
 // 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/ic2par2ic.hpp"
-#include "kep3/core_astro/propagate_lagrangian.hpp"
 #include <chrono>
 #include <iostream>
+#include <kep3/core_astro/ic2par2ic.hpp>
+#include <kep3/core_astro/propagate_lagrangian.hpp>
 #include <random>
 
 #include <fmt/core.h>
@@ -46,7 +46,7 @@ void perform_test_speed(
   std::uniform_real_distribution<double> Omega_d(0, 2 * pi);
   std::uniform_real_distribution<double> omega_d(0., 2 * pi);
   std::uniform_real_distribution<double> f_d(0, 2 * pi);
-  std::uniform_real_distribution<double> tof_d(0.1, 20.);
+  std::uniform_real_distribution<double> tof_d(10., 100.);
 
   // We generate the random dataset
   std::vector<std::array<std::array<double, 3>, 2>> pos_vels(N);
@@ -86,17 +86,17 @@ void perform_test_accuracy(
   // Engines
   //
   // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
-  std::mt19937 rng_engine(122012203u);
+  std::mt19937 rng_engine(53534535u);
   //
   // Distributions
   //
-  std::uniform_real_distribution<double> sma_d(0.5, 20.);
+  std::uniform_real_distribution<double> sma_d(0.5, 10.);
   std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
   std::uniform_real_distribution<double> incl_d(0., pi);
   std::uniform_real_distribution<double> Omega_d(0, 2 * pi);
   std::uniform_real_distribution<double> omega_d(0., 2 * pi);
   std::uniform_real_distribution<double> f_d(0, 2 * pi);
-  std::uniform_real_distribution<double> tof_d(0.1, 20.);
+  std::uniform_real_distribution<double> tof_d(10, 100.);
 
   // We generate the random dataset
   std::vector<std::array<std::array<double, 3>, 2>> pos_vels(N);
@@ -113,7 +113,8 @@ void perform_test_accuracy(
     }
 
     pos_vels[i] = kep3::par2ic({sma, ecc, incl_d(rng_engine),
-                                Omega_d(rng_engine), omega_d(rng_engine), f}, 1.);
+                                Omega_d(rng_engine), omega_d(rng_engine), f},
+                               1.);
     tofs[i] = tof_d(rng_engine);
   }
   // We log progress
@@ -149,17 +150,29 @@ int main() {
   perform_test_accuracy(0, 0.5, 100000, &kep3::propagate_lagrangian);
   perform_test_accuracy(0.5, 0.9, 100000, &kep3::propagate_lagrangian);
   perform_test_accuracy(0.9, 0.99, 100000, &kep3::propagate_lagrangian);
-
-  fmt::print("\nComputes speed at different eccentricity ranges [Universal "
-             "Anomaly]:\n");
-  perform_test_speed(0, 0.5, 1000000, &kep3::propagate_lagrangian_u);
-  perform_test_speed(0.5, 0.9, 1000000, &kep3::propagate_lagrangian_u);
-  perform_test_speed(0.9, 0.99, 1000000, &kep3::propagate_lagrangian_u);
-  perform_test_speed(1.1, 10., 1000000, &kep3::propagate_lagrangian_u);
-
-  fmt::print("\nComputes error at different eccentricity ranges [Universal "
-             "Anomaly]:\n");
-  perform_test_accuracy(0, 0.5, 100000, &kep3::propagate_lagrangian_u);
-  perform_test_accuracy(0.5, 0.9, 100000, &kep3::propagate_lagrangian_u);
-  perform_test_accuracy(0.9, 0.99, 100000, &kep3::propagate_lagrangian_u);
+//
+  //fmt::print("\nComputes speed at different eccentricity ranges [Universal "
+  //           "Anomaly]:\n");
+  //perform_test_speed(0, 0.5, 1000000, &kep3::propagate_lagrangian_u);
+  //perform_test_speed(0.5, 0.9, 1000000, &kep3::propagate_lagrangian_u);
+  //perform_test_speed(0.9, 0.99, 1000000, &kep3::propagate_lagrangian_u);
+  //perform_test_speed(1.1, 10., 1000000, &kep3::propagate_lagrangian_u);
+//
+  //fmt::print("\nComputes error at different eccentricity ranges [Universal "
+  //           "Anomaly]:\n");
+  //perform_test_accuracy(0, 0.5, 100000, &kep3::propagate_lagrangian_u);
+  //perform_test_accuracy(0.5, 0.9, 100000, &kep3::propagate_lagrangian_u);
+  //perform_test_accuracy(0.9, 0.99, 100000, &kep3::propagate_lagrangian_u);
+//
+  //fmt::print("\nComputes speed at different eccentricity ranges [keplerian "
+  //           "propagation]:\n");
+  //perform_test_speed(0, 0.5, 1000000, &kep3::propagate_keplerian);
+  //perform_test_speed(0.5, 0.9, 1000000, &kep3::propagate_keplerian);
+  //perform_test_speed(0.9, 0.99, 1000000, &kep3::propagate_keplerian);
+//
+  //fmt::print("\nComputes error at different eccentricity ranges [keplerian "
+  //           "propagation]:\n");
+  //perform_test_accuracy(0, 0.5, 100000, &kep3::propagate_keplerian);
+  //perform_test_accuracy(0.5, 0.9, 100000, &kep3::propagate_keplerian);
+  //perform_test_accuracy(0.9, 0.99, 100000, &kep3::propagate_keplerian);
 }

codecov.yml

@@ -4,4 +4,9 @@ coverage:
 codecov: 
  token: 60883dbc-866d-4da4-9306-8df7ac9b1360
 
-comment: off
+comment:                  # this is a top-level key
+  layout: " diff, flags, files"
+  behavior: default
+  require_changes: false  # if true: only post the comment if coverage changes
+  require_base: false        # [true :: must have a base report to post]
+  require_head: true       # [true :: must have a head report to post]

include/kep3/core_astro/convert_anomalies.hpp

@@ -34,15 +34,16 @@ inline double m2e(double M, double ecc) {
   if (M_cropped < 0) {
     M_cropped += 2 * kep3::pi;
   }
-  // The Initial guess follows from a third order expansion of Kepler's equation.
+  // The Initial guess follows from a third order expansion of Kepler's
+  // equation.
   double IG = M_cropped + ecc * sinM + ecc * ecc * sinM * cosM +
               ecc * ecc * ecc * sinM * (1.5 * cosM * cosM - 0.5);
 
   const int digits = std::numeric_limits<double>::digits;
   std::uintmax_t max_iter = 100u;
 
-  // Newton-raphson or Halley iterates can be used here. Similar performances, thus
-  // we choose the simplest algorithm.
+  // Newton-raphson or Halley iterates can be used here. Similar performances,
+  // thus we choose the simplest algorithm.
   double sol = boost::math::tools::newton_raphson_iterate(
       [M_cropped, ecc](double E) {
         return std::make_tuple(kepE(E, M_cropped, ecc), d_kepE(E, ecc));
@@ -62,23 +63,70 @@ inline double e2m(double E, double ecc) { return (E - ecc * std::sin(E)); }
 inline double e2f(double E, double ecc) {
   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) {
   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); }
+
 // 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); }
 
 // gudermannian to true (only hyperbolas) e>1 (returns in range [-pi,pi])
 inline double zeta2f(double f, double ecc) {
   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) {
   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) {
+
+  // The Initial guess (TODO(darioizo) improve)
+  double IG = 1.;
+
+  const int digits = std::numeric_limits<double>::digits;
+  std::uintmax_t max_iter = 100u;
+
+  // Newton-raphson iterates.
+  double sol = boost::math::tools::newton_raphson_iterate(
+      [N, ecc](double H) {
+        return std::make_tuple(kepH(H, N, ecc), d_kepH(H, ecc));
+      },
+      IG, IG - 20 * kep3::pi, IG + 20 * kep3::pi, digits, max_iter);
+  if (max_iter == 100u) {
+    throw std::domain_error(
+        "Maximum number of iterations exceeded when solving Kepler's "
+        "equation for the hyperbolic anomaly in m2h.");
+  }
+  return sol;
+}
+
+// hyperbolic H to hyperbolic mean N (only hyperbolas) e>1
+inline double h2n(double H, double ecc) { 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) {
+  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) {
+  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); }
+
+// 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); }
+
 } // namespace kep3
 #endif // kep3_TOOLBOX_M2E_H

include/kep3/core_astro/propagate_lagrangian.hpp

@@ -13,6 +13,7 @@
 #include <array>
 #include <cmath>
 
+
 #include<kep3/detail/visibility.hpp>
 
 namespace kep3 {
@@ -30,6 +31,9 @@ kep3_DLL_PUBLIC void propagate_lagrangian(std::array<std::array<double, 3>, 2> &
 
 kep3_DLL_PUBLIC void propagate_lagrangian_u(std::array<std::array<double, 3>, 2> &pos_vel,
                           double dt, double mu);
+
+kep3_DLL_PUBLIC void propagate_keplerian(std::array<std::array<double, 3>, 2> &pos_vel,
+                          double dt, double mu);
 } // namespace kep3
 
 #endif // KEP_TOOLBOX_PROPAGATE_LAGRANGIAN_H