tmp

include/kep3/planets/keplerian.hpp

@@ -18,30 +18,31 @@
 #include <kep3/epoch.hpp>
 #include <kep3/planet.hpp>
 
-
 namespace kep3::udpla {
 
 class kep3_DLL_PUBLIC keplerian {
 
   static const std::array<double, 6> default_elements;
 
   kep3::epoch m_ref_epoch;
-  std::array<double, 6> m_elem;
+  std::array<std::array<double, 3>, 2> m_pos_vel_0;
   std::string m_name;
   double m_mu_central_body;
   double m_mu_self;
   double m_radius;
   double m_safe_radius;
+  bool m_ellipse;
 
   friend class boost::serialization::access;
   template <typename Archive> void serialize(Archive &ar, unsigned) {
     ar &m_ref_epoch;
-    ar &m_elem;
+    ar &m_pos_vel_0;
     ar &m_name;
     ar &m_mu_central_body;
     ar &m_mu_self;
     ar &m_radius;
     ar &m_safe_radius;
+    ar &m_ellipse;
   }
 
 public:
@@ -51,9 +52,13 @@ class kep3_DLL_PUBLIC keplerian {
                      const std::array<double, 6> &elem = default_elements,
                      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),
+                     const std::array<std::array<double, 3>, 2> &pos_vel =
+                         {{{1.0, 0.0, 0.0}, {0., 1.0, 0.0}}},
+                     double mu_central_body = 1., std::string name = "Unknown",
+                     std::array<double, 3> added_params = {-1., -1., -1.});
   // Mandatory UDPLA methods
-  static std::array<std::array<double, 3>, 2> eph(const epoch &);
+  [[nodiscard]] std::array<std::array<double, 3>, 2> eph(const epoch &) const;
 
   // Optional UDPLA methods
   [[nodiscard]] std::string get_name() const;

src/planets/keplerian.cpp

@@ -0,0 +1,57 @@
+// 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 <string>
+
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/roots.hpp>
+
+#include <kep3/core_astro/kepler_equations.hpp>
+#include <kep3/core_astro/propagate_lagrangian.hpp>
+#include <kep3/epoch.hpp>
+#include <kep3/planets/keplerian.hpp>
+#include <utility>
+
+namespace kep3::udpla {
+
+constexpr double pi{boost::math::constants::pi<double>()};
+const double DAY2SEC = 24.*60.*60.;
+
+keplerian::keplerian(const epoch &ref_epoch,
+                     const std::array<std::array<double, 3>, 2> &pos_vel,
+                     double mu_central_body, std::string name,
+                     std::array<double, 3> added_params)
+    : m_ref_epoch(ref_epoch), m_name(std::move(name)),
+      m_mu_central_body(mu_central_body), m_mu_self(added_params[0]),
+      m_radius(added_params[1]), m_ellipse(), m_safe_radius(added_params[2]),
+      m_pos_vel_0(pos_vel) {
+  double R =
+      std::sqrt(pos_vel[0][0] * pos_vel[0][0] + pos_vel[0][1] * pos_vel[0][1] +
+                pos_vel[0][2] * pos_vel[0][2]);
+  double en = (pos_vel[1][0] * pos_vel[1][0] + pos_vel[1][1] * pos_vel[1][1] +
+              pos_vel[1][2] * pos_vel[1][2]) /
+                 2. -
+             mu_central_body / R;
+  (en>0) ? m_ellipse=false : m_ellipse=true;
+}
+
+std::array<std::array<double, 3>, 2>
+keplerian::eph(const kep3::epoch &ep) const {
+  // 1 - We compute the dt
+  double dt = (ep.mjd2000() - m_ref_epoch.mjd2000()) * DAY2SEC;
+  // 2 - We propagate (make a copy as we do not want to change m_pos_vel_0)
+  auto retval(m_pos_vel_0);
+  kep3::propagate_lagrangian(retval, dt, m_mu_central_body);
+  return retval;
+}
+
+const std::array<double, 6> keplerian::default_elements = {1.0, 0.0, 0.0,
+                                                           0.0, 0.0, 0.0};
+
+} // namespace kep3::udpla