keplerian test

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/planets/keplerian.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/core_astro/constants.hpp

@@ -0,0 +1,48 @@
+// 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_CONSTANTS_H
+#define kep3_CONSTANTS_H
+
+#include <cmath>
+
+#include <boost/math/constants/constants.hpp>
+
+namespace kep3 {
+
+enum elements_type {
+  KEP_M,  // Keplerian Osculating (with Mean Anomaly)
+  KEP_F,  // Keplerian Osculating (with True Anomaly)
+  MEQ,    // Modified Equinoctial Elements
+  MEQ_R,  // Modified Equinoctial Elements (retrogade)
+  POSVEL, // position and Velocity
+};
+
+constexpr double pi = boost::math::constants::pi<double>();
+constexpr double half_pi = boost::math::constants::half_pi<double>();
+
+constexpr double AU = 149597870700.0;     // Astronomical Unit (m)
+constexpr double CAVENDISH = 73.6687e-11; // Cavendish constant
+constexpr double MU_SUN =
+    1.32712440018e20; // Sun's gravitational parameter (m^3/s^2 kg)
+constexpr double MU_EARTH =
+    398600441800000.0; // Earth's gravitational parameter (m^3/s^2 kg)
+constexpr double EARTH_VELOCITY = 2.97846918317e4; // Earth's velocity. (m/s)
+constexpr double EARTH_J2 = 1.08262668E-03;
+constexpr double EARTH_RADIUS = 6378137; // Earth's radius (m)
+constexpr double DEG2RAD = (pi / 180.0);
+constexpr double RAD2DEG = (180.0 / pi);
+constexpr double DAY2SEC = 86400.0;
+constexpr double SEC2DAY = (1. / DAY2SEC);
+constexpr double DAY2YEAR = (1. / 365.25);
+constexpr double G0 = 9.80665; // Acceleration at Earth's surface (m/s^2)
+
+} // namespace kep3
+
+#endif // kep3_CONSTANTS_H

include/kep3/detail/s11n.hpp

@@ -28,6 +28,7 @@
 #include <boost/serialization/unique_ptr.hpp>
 #include <boost/serialization/utility.hpp>
 #include <boost/serialization/vector.hpp>
+#include <boost/serialization/array.hpp>
 #include <boost/serialization/version.hpp>
 
 // Include the archives.

include/kep3/epoch.hpp

@@ -10,10 +10,13 @@
 #ifndef kep3_EPOCH_H
 #define kep3_EPOCH_H
 
+#include <iostream>
+
+#include <fmt/ostream.h>
+
 #include <boost/date_time/gregorian/gregorian.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
 #include <boost/lexical_cast.hpp>
-#include <iostream>
 
 #include <kep3/detail/s11n.hpp>
 #include <kep3/detail/visibility.hpp>
@@ -106,4 +109,6 @@ kep3_DLL_PUBLIC double operator-(const epoch &lhs, const epoch &rhs);
 
 } // end of namespace kep3
 
+template <> struct fmt::formatter<kep3::epoch> : ostream_formatter {};
+
 #endif // kep3_EPOCH_H

include/kep3/planets/keplerian.hpp

@@ -14,6 +14,9 @@
 #include <array>
 #include <utility>
 
+#include <fmt/ostream.h>
+
+#include <kep3/core_astro/constants.hpp>
 #include <kep3/detail/visibility.hpp>
 #include <kep3/epoch.hpp>
 #include <kep3/planet.hpp>
@@ -22,8 +25,6 @@ namespace kep3::udpla {
 
 class kep3_DLL_PUBLIC keplerian {
 
-  static const std::array<double, 6> default_elements;
-
   kep3::epoch m_ref_epoch;
   std::array<std::array<double, 3>, 2> m_pos_vel_0;
   std::string m_name;
@@ -48,8 +49,7 @@ class kep3_DLL_PUBLIC keplerian {
 public:
   // NOTE: in here elem is a,e,i,W,w,M (Mean anomaly, not true anomaly)
   // NOTE: added_param contains mu_self, radius and safe_radius
-  explicit keplerian(const epoch &ref_epoch = kep3::epoch(0),
-                     const std::array<double, 6> &elem = default_elements,
+  explicit keplerian(const epoch &ref_epoch, const std::array<double, 6> &par,
                      double mu_central_body = 1., std::string name = "Unknown",
                      std::array<double, 3> added_params = {-1., -1., -1.});
   explicit keplerian(const epoch &ref_epoch = kep3::epoch(0),
@@ -69,10 +69,16 @@ class kep3_DLL_PUBLIC keplerian {
   [[nodiscard]] std::string get_extra_info() const;
 
   // Other methods
-  [[nodiscard]] std::string get_ref_epoch() const;
-  [[nodiscard]] std::string get_elem() const;
+  [[nodiscard]] kep3::epoch get_ref_epoch() const;
+  [[nodiscard]] std::array<double, 6>
+      elements(kep3::elements_type = kep3::elements_type::KEP_F) const;
 };
+kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &,
+                                         const kep3::udpla::keplerian &);
 } // namespace kep3::udpla
+
+template <>
+struct fmt::formatter<kep3::udpla::keplerian> : ostream_formatter {};
 kep3_S11N_PLANET_EXPORT_KEY(kep3::udpla::keplerian);
 
 #endif // kep3_EPOCH_H

src/core_astro/eq2par2eq.cpp

@@ -10,15 +10,11 @@
 #include <array>
 #include <cmath>
 
-#include <boost/math/constants/constants.hpp>
-
+#include <kep3/core_astro/constants.hpp>
 #include <kep3/core_astro/eq2par2eq.hpp>
 
 namespace kep3 {
 
-constexpr double half_pi{boost::math::constants::half_pi<double>()};
-constexpr double pi{boost::math::constants::pi<double>()};
-
 std::array<double, 6> eq2par(const std::array<double, 6> &eq, bool retrogade) {
   std::array<double, 6> retval{};
   int I = 1;

src/core_astro/ic2eq2ic.cpp

@@ -10,24 +10,20 @@
 #include <array>
 #include <cmath>
 
-#include <boost/math/constants/constants.hpp>
-
 #include <fmt/core.h>
 #include <fmt/ranges.h>
 
-
 #include <xtensor-blas/xlinalg.hpp>
 #include <xtensor/xadapt.hpp>
 
+#include <kep3/core_astro/constants.hpp>
 #include <kep3/core_astro/ic2eq2ic.hpp>
 
 using xt::linalg::cross;
 using xt::linalg::dot;
 
 namespace kep3 {
 
-constexpr double half_pi{boost::math::constants::half_pi<double>()};
-
 // Implementation following:
 // Cefola: Equinoctial orbit elements - Application to artificial satellite
 // orbitsCefola, P., 1972, September. Equinoctial orbit elements-Application to
@@ -108,54 +104,52 @@ std::array<double, 6> ic2eq(const std::array<std::array<double, 3>, 2> &pos_vel,
   }
 }
 
-std::array<std::array<double, 3>, 2>
-eq2ic(const std::array<double, 6> &eq, double mu, bool retrogade)
-{
-    std::array<std::array<double, 3>, 2> retval{};
-    int I = 0;
-    if (retrogade) {
-        I = -1;
-    } else {
-        I = 1;
-    }
+std::array<std::array<double, 3>, 2> eq2ic(const std::array<double, 6> &eq,
+                                           double mu, bool retrogade) {
+  std::array<std::array<double, 3>, 2> retval{};
+  int I = 0;
+  if (retrogade) {
+    I = -1;
+  } else {
+    I = 1;
+  }
 
-    // p = a (1-e^2) will be negative for eccentricities > 1, we here need a positive number for the following
-    // computations
-    // to make sense
-    double par = std::abs(eq[0]);
-    double f = eq[1];
-    double g = eq[2];
-    double h = eq[3];
-    double k = eq[4];
-    double L = eq[5];
-
-    // We compute the equinoctial reference frame
-    double den = k * k + h * h + 1;
-    double fx = (1 - k * k + h * h) / den;
-    double fy = (2 * k * h) / den;
-    double fz = (-2 * I * k) / den;
-
-    double gx = (2 * I * k * h) / den;
-    double gy = (1 + k * k - h * h) * I / den;
-    double gz = (2 * h) / den;
-
-    // Auxiliary
-    double radius = par / (1 + g * std::sin(L) + f * std::cos(L));
-    // In the equinoctial reference frame
-    double X = radius * std::cos(L);
-    double Y = radius * std::sin(L);
-    double VX = -std::sqrt(mu / par) * (g + std::sin(L));
-    double VY = std::sqrt(mu / par) * (f + std::cos(L));
-
-    // Results
-    retval[0][0] = X * fx + Y * gx;
-    retval[0][1] = X * fy + Y * gy;
-    retval[0][2] = X * fz + Y * gz;
-
-    retval[1][0] = VX * fx + VY * gx;
-    retval[1][1] = VX * fy + VY * gy;
-    retval[1][2] = VX * fz + VY * gz;
-
-    return retval;
+  // p = a (1-e^2) will be negative for eccentricities > 1, we here need a
+  // positive number for the following computations to make sense
+  double par = std::abs(eq[0]);
+  double f = eq[1];
+  double g = eq[2];
+  double h = eq[3];
+  double k = eq[4];
+  double L = eq[5];
+
+  // We compute the equinoctial reference frame
+  double den = k * k + h * h + 1;
+  double fx = (1 - k * k + h * h) / den;
+  double fy = (2 * k * h) / den;
+  double fz = (-2 * I * k) / den;
+
+  double gx = (2 * I * k * h) / den;
+  double gy = (1 + k * k - h * h) * I / den;
+  double gz = (2 * h) / den;
+
+  // Auxiliary
+  double radius = par / (1 + g * std::sin(L) + f * std::cos(L));
+  // In the equinoctial reference frame
+  double X = radius * std::cos(L);
+  double Y = radius * std::sin(L);
+  double VX = -std::sqrt(mu / par) * (g + std::sin(L));
+  double VY = std::sqrt(mu / par) * (f + std::cos(L));
+
+  // Results
+  retval[0][0] = X * fx + Y * gx;
+  retval[0][1] = X * fy + Y * gy;
+  retval[0][2] = X * fz + Y * gz;
+
+  retval[1][0] = VX * fx + VY * gx;
+  retval[1][1] = VX * fy + VY * gy;
+  retval[1][2] = VX * fz + VY * gz;
+
+  return retval;
 }
 } // namespace kep3

src/core_astro/ic2par2ic.cpp

@@ -10,20 +10,17 @@
 #include <array>
 #include <cmath>
 
-#include <boost/math/constants/constants.hpp>
-
 #include <xtensor-blas/xlinalg.hpp>
 #include <xtensor/xadapt.hpp>
 
+#include <kep3/core_astro/constants.hpp>
 #include <kep3/core_astro/ic2par2ic.hpp>
 
 using xt::linalg::cross;
 using xt::linalg::dot;
 
 namespace kep3 {
 
-constexpr double pi{boost::math::constants::pi<double>()};
-
 // r,v,mu -> keplerian osculating elements [a,e,i,W,w,f]. The last
 // is the true anomaly. The semi-major axis a is positive for ellipses, negative
 // for hyperbolae. The anomalies W, w, f are in [0, 2pi]. Inclination is in [0,

src/core_astro/propagate_lagrangian.cpp

@@ -9,19 +9,17 @@
 
 #include <array>
 #include <cmath>
+#include <stdexcept>
 
-#include <boost/math/constants/constants.hpp>
 #include <boost/math/tools/roots.hpp>
 
+#include <kep3/core_astro/constants.hpp>
 #include <kep3/core_astro/kepler_equations.hpp>
 #include <kep3/core_astro/propagate_lagrangian.hpp>
 #include <kep3/core_astro/special_functions.hpp>
-#include <stdexcept>
 
 namespace kep3 {
 
-constexpr double pi{boost::math::constants::pi<double>()};
-
 /// Lagrangian propagation
 /**
  * This function propagates an initial Cartesian state for a time t assuming a

src/planet.cpp

@@ -7,6 +7,7 @@
 // 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/detail/exceptions.hpp>
 #include <kep3/planet.hpp>