Merge pull request #7 from darioizzo/jpl_lp

Jpl lp

CMakeLists.txt

@@ -133,6 +133,7 @@ set(kep3_SRC_FILES
     "${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/planets/jpl_lp.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"

include/kep3/planet.hpp

@@ -228,7 +228,7 @@ struct kep3_DLL_PUBLIC_INLINE_CLASS planet_inner final : planet_inner_base {
       // period from the energy at epoch
       auto [r, v] = eph(ep);
       double mu = get_mu_central_body();
-      double R = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+      double R = std::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
       double v2 = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
       double en = v2 / 2. - mu / R;
       if (en > 0) {

include/kep3/planets/jpl_lp.hpp

@@ -0,0 +1,87 @@
+// Copyright 2023, 2024 Dario Izzo ([email protected]), Francesco Biscani
+// ([email protected])
+//
+// 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 kep3_PLANET_JPL_LP_H
+#define kep3_PLANET_JPL_LP_H
+
+#include <array>
+
+#include <fmt/ostream.h>
+
+#include <kep3/core_astro/constants.hpp>
+#include <kep3/detail/s11n.hpp>
+#include <kep3/detail/visibility.hpp>
+#include <kep3/epoch.hpp>
+#include <kep3/planet.hpp>
+
+namespace kep3::udpla {
+
+/// Solar System Planet (jpl simplified ephemerides)
+/**
+ * This class allows to instantiate planets of
+ * the solar system by referring to their common names. The ephemeris used
+ * are low_precision ephemeris defined in the JPL pages:
+ * https://ssd.jpl.nasa.gov/planets/approx_pos.html and valid for the timeframe
+ * 1800AD - 2050 AD
+ */
+
+class kep3_DLL_PUBLIC jpl_lp {
+
+  std::array<double, 6> m_elements;
+  std::array<double, 6> m_elements_dot;
+  std::string m_name;
+  double m_mu_central_body;
+  double m_mu_self;
+  double m_radius;
+  double m_safe_radius;
+
+  friend class boost::serialization::access;
+  template <typename Archive> void serialize(Archive &ar, unsigned) {
+    ar &m_elements;
+    ar &m_elements_dot;
+    ar &m_name;
+    ar &m_mu_central_body;
+    ar &m_mu_self;
+    ar &m_radius;
+    ar &m_safe_radius;
+  }
+
+public:
+  // Constructor
+  explicit jpl_lp(const std::string & = "earth");
+  // Mandatory UDPLA methods
+  [[nodiscard]] std::array<std::array<double, 3>, 2> eph(const epoch &) const;
+
+  // Optional UDPLA methods
+  [[nodiscard]] std::string get_name() const;
+  [[nodiscard]] double get_mu_central_body() const;
+  [[nodiscard]] double get_mu_self() const;
+  [[nodiscard]] double get_radius() const;
+  [[nodiscard]] double get_safe_radius() const;
+  [[nodiscard]] std::string get_extra_info() const;
+
+  // Other methods
+  [[nodiscard]] std::array<double, 6>
+  elements(const kep3::epoch & = kep3::epoch(),
+           kep3::elements_type = kep3::elements_type::KEP_F) const;
+
+private:
+  [[nodiscard]] std::array<double, 6>
+  _f_elements(const kep3::epoch & = kep3::epoch()) const;
+};
+kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &,
+                                         const kep3::udpla::jpl_lp &);
+} // namespace kep3::udpla
+
+// fmt formatter redirecting to the stream operator
+template <> struct fmt::formatter<kep3::udpla::jpl_lp> : ostream_formatter {};
+// necessary for serialization
+kep3_S11N_PLANET_EXPORT_KEY(kep3::udpla::jpl_lp);
+
+#endif // KEP_TOOLBOX_PLANET_JPL_LP_H

include/kep3/planets/keplerian.hpp

@@ -1,4 +1,3 @@
-
 // Copyright 2023, 2024 Dario Izzo ([email protected]), Francesco Biscani
 // ([email protected])
 //
@@ -80,8 +79,10 @@ kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &,
                                          const kep3::udpla::keplerian &);
 } // namespace kep3::udpla
 
+// fmt formatter redirecting to the stream operator
 template <>
 struct fmt::formatter<kep3::udpla::keplerian> : ostream_formatter {};
+// necessary for serialization
 kep3_S11N_PLANET_EXPORT_KEY(kep3::udpla::keplerian);
 
-#endif // kep3_EPOCH_H
+#endif // kep3_UDPLA_KEPLERIAN_H

src/planets/jpl_lp.cpp

@@ -0,0 +1,212 @@
+// Copyright 2023, 2024 Dario Izzo ([email protected]), Francesco Biscani
+// ([email protected])
+//
+// 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 "kep3/core_astro/eq2par2eq.hpp"
+#include "kep3/epoch.hpp"
+#include <cmath>
+#include <stdexcept>
+#include <unordered_map>
+
+#include <boost/algorithm/string.hpp>
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <kep3/core_astro/constants.hpp>
+#include <kep3/core_astro/convert_anomalies.hpp>
+#include <kep3/planet.hpp>
+#include <kep3/planets/jpl_lp.hpp>
+
+namespace kep3::udpla {
+
+// Data from: https://ssd.jpl.nasa.gov/planets/approx_pos.html
+// a,e,i,L,W,w
+// clang-format off
+static const std::array<double, 6> mercury_el = {0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593};
+static const std::array<double, 6> mercury_el_dot = {0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081};
+static const std::array<double, 6> venus_el = {0.72333566, 0.00677672, 3.39467605, 181.97909950, 131.60246718, 76.67984255};
+static const std::array<double, 6> venus_el_dot = {0.00000390, -0.00004107, -0.00078890, 58517.81538729, 0.00268329, -0.27769418};
+static const std::array<double, 6> earth_moon_el = {1.00000261, 0.01671123, -0.00001531, 100.46457166, 102.93768193, 0.0};
+static const std::array<double, 6> earth_moon_el_dot = {0.00000562, -0.00004392, -0.01294668, 35999.37244981, 0.32327364, 0.0};
+static const std::array<double, 6> mars_el = {1.52371034, 0.09339410, 1.84969142, -4.55343205, -23.94362959, 49.55953891};
+static const std::array<double, 6> mars_el_dot = {0.00001847, 0.00007882, -0.00813131, 19140.30268499, 0.44441088, -0.29257343};
+static const std::array<double, 6> jupiter_el = {5.20288700, 0.04838624, 1.30439695, 34.39644051, 14.72847983, 100.47390909};
+static const std::array<double, 6> jupiter_el_dot = {-0.00011607, -0.00013253, -0.00183714, 3034.74612775, 0.21252668, 0.20469106};
+static const std::array<double, 6> saturn_el = {9.53667594, 0.05386179, 2.48599187, 49.95424423, 92.59887831, 113.66242448};
+static const std::array<double, 6> saturn_el_dot = {-0.00125060, -0.00050991, 0.00193609, 1222.49362201, -0.41897216, -0.28867794};
+static const std::array<double, 6> uranus_el = {19.18916464, 0.04725744, 0.77263783, 313.23810451, 170.95427630, 74.01692503};
+static const std::array<double, 6> uranus_el_dot = {-0.00196176, -0.00004397, -0.00242939, 428.48202785, 0.40805281, 0.04240589};
+static const std::array<double, 6> neptune_el = {30.06992276, 0.00859048, 1.77004347, -55.12002969, 44.96476227, 131.78422574};
+static const std::array<double, 6> 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) {
+
+  // Cannot be marked const as operator[] is not and we use it.
+  static std::unordered_map<std::string, int> mapped_planets = {
+      {"mercury", 1}, {"venus", 2},   {"earth", 3},
+      {"mars", 4},    {"jupiter", 5}, {"saturn", 6},
+      {"uranus", 7},  {"neptune", 8}, {"pluto", 9}};
+
+  boost::algorithm::to_lower(m_name);
+  switch (mapped_planets[m_name]) {
+  case (1): {
+    m_elements = mercury_el;
+    m_elements_dot = mercury_el_dot;
+    m_radius = 2440000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 22032e9;
+  } break;
+  case (2): {
+    m_elements = venus_el;
+    m_elements_dot = venus_el_dot;
+    m_radius = 6052000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 324859e9;
+  } break;
+  case (3): {
+    m_elements = earth_moon_el;
+    m_elements_dot = earth_moon_el_dot;
+    m_radius = 6378000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 398600.4418e9;
+  } break;
+  case (4): {
+    m_elements = mars_el;
+    m_elements_dot = mars_el_dot;
+    m_radius = 3397000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 42828e9;
+  } break;
+  case (5): {
+    m_elements = jupiter_el;
+    m_elements_dot = jupiter_el_dot;
+    m_radius = 71492000.;
+    m_safe_radius = 9. * m_radius;
+    m_mu_self = 126686534e9;
+  } break;
+  case (6): {
+    m_elements = saturn_el;
+    m_elements_dot = saturn_el_dot;
+    m_radius = 60330000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 37931187e9;
+  } break;
+  case (7): {
+    m_elements = uranus_el;
+    m_elements_dot = uranus_el_dot;
+    m_radius = 25362000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 5793939e9;
+  } break;
+  case (8): {
+    m_elements = neptune_el;
+    m_elements_dot = neptune_el_dot;
+    m_radius = 24622000.;
+    m_safe_radius = 1.1 * m_radius;
+    m_mu_self = 6836529e9;
+  } break;
+  default: {
+    throw std::logic_error("unknown planet name: ");
+  }
+  }
+}
+
+// Computes the kep3::KEP_F elements (osculating with true anomaly) at epoch.
+std::array<double, 6> jpl_lp::_f_elements(const kep3::epoch &ep) const {
+  double mjd2000 = ep.mjd2000();
+  if (mjd2000 <= -73048.0 || mjd2000 >= 18263.0) {
+    throw std::domain_error("Low precision Ephemeris are only valid in the "
+                            "range range [1800-2050]");
+  }
+  // algorithm from https://ssd.jpl.nasa.gov/planets/approx_pos.html accessed
+  // 2023.
+  std::array<double, 6> elements_updated{}, elements_f{};
+  double dt = mjd2000 / 36525.; // Number of centuries passed since J2000.0
+  for (unsigned int i = 0; i < 6; ++i) {
+    elements_updated[i] = (m_elements[i] + m_elements_dot[i] * dt);
+  }
+  elements_f[0] = elements_updated[0] * kep3::AU;
+  elements_f[1] = elements_updated[1];
+  elements_f[2] = elements_updated[2] * kep3::DEG2RAD;
+  elements_f[3] = elements_updated[5] * kep3::DEG2RAD;
+  elements_f[4] = (elements_updated[4] - elements_updated[5]) * kep3::DEG2RAD;
+  elements_f[5] = (elements_updated[3] - elements_updated[4]) * kep3::DEG2RAD;
+  elements_f[5] = kep3::m2f(elements_f[5], elements_f[1]);
+  return elements_f;
+}
+
+std::array<std::array<double, 3>, 2> jpl_lp::eph(const epoch &ep) const {
+  auto elements_f = _f_elements(ep);
+  return par2ic(elements_f, get_mu_central_body());
+}
+
+std::array<double, 6> jpl_lp::elements(const kep3::epoch &ep,
+                                       kep3::elements_type el_type) const {
+  auto elements = _f_elements(ep);
+  switch (el_type) {
+  case kep3::elements_type::KEP_F:
+    break;
+  case kep3::elements_type::KEP_M:
+    elements[5] = kep3::f2m(elements[5], elements[1]);
+    break;
+  case kep3::elements_type::MEQ:
+    elements = kep3::par2eq(elements, false);
+    break;
+  case kep3::elements_type::MEQ_R:
+    elements = kep3::par2eq(elements, true);
+    break;
+  default:
+    throw std::logic_error("You should not go here!");
+  }
+  return elements;
+}
+
+std::string jpl_lp::get_name() const { return m_name; }
+
+double jpl_lp::get_mu_central_body() const { return m_mu_central_body; }
+
+double jpl_lp::get_mu_self() const { return m_mu_self; }
+
+double jpl_lp::get_radius() const { return m_radius; }
+
+double jpl_lp::get_safe_radius() const { return m_safe_radius; }
+
+std::string jpl_lp::get_extra_info() const {
+  kep3::epoch ep{0., kep3::epoch::MJD2000};
+  auto par = elements(ep);
+  auto pos_vel = eph(ep);
+
+  std::string retval =
+      fmt::format("\nLow-precision planet elements: \n") +
+      fmt::format("Semi major axis (AU): {}\n", par[0] / kep3::AU) +
+      fmt::format("Eccentricity: {}\n", par[1]) +
+      fmt::format("Inclination (deg.): {}\n", par[2] * kep3::RAD2DEG) +
+      fmt::format("Big Omega (deg.): {}\n", par[3] * kep3::RAD2DEG) +
+      fmt::format("Small omega (deg.): {}\n", par[4] * kep3::RAD2DEG) +
+      fmt::format("True anomly (deg.): {}\n", par[5] * kep3::RAD2DEG);
+  retval += fmt::format("Mean anomly (deg.): {}\n",
+                        kep3::f2m(par[5], par[1]) * kep3::RAD2DEG);
+  retval +=
+      fmt::format("Elements reference epoch (MJD2000): {}\n", ep.mjd2000()) +
+      fmt::format("Elements reference epoch (date): {}\n", ep) +
+      fmt::format("r at ref. = {}\n", pos_vel[0]) +
+      fmt::format("v at ref. = {}\n", pos_vel[1]);
+  return retval;
+}
+
+std::ostream &operator<<(std::ostream &os, const kep3::udpla::jpl_lp &udpla) {
+  os << udpla.get_extra_info() << std::endl;
+  return os;
+}
+
+} // namespace kep3::udpla
+
+kep3_S11N_PLANET_IMPLEMENT(kep3::udpla::jpl_lp)

src/planets/keplerian.cpp

@@ -7,7 +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 <cmath>
 #include <limits>
 #include <stdexcept>
@@ -18,6 +17,7 @@
 
 #include <kep3/core_astro/constants.hpp>
 #include <kep3/core_astro/convert_anomalies.hpp>
+#include <kep3/core_astro/eq2par2eq.hpp>
 #include <kep3/core_astro/ic2eq2ic.hpp>
 #include <kep3/core_astro/ic2par2ic.hpp>
 #include <kep3/core_astro/propagate_lagrangian.hpp>

test/CMakeLists.txt

@@ -30,6 +30,7 @@ ADD_kep3_TESTCASE(convert_anomalies_test)
 ADD_kep3_TESTCASE(epoch_test)
 ADD_kep3_TESTCASE(planet_test)
 ADD_kep3_TESTCASE(planet_keplerian_test)
+ADD_kep3_TESTCASE(planet_jpl_lp_test)
 ADD_kep3_TESTCASE(ic2par2ic_test)
 ADD_kep3_TESTCASE(ic2eq2ic_test)
 ADD_kep3_TESTCASE(eq2par2eq_test)