ic2par and opposite started

CMakeLists.txt

@@ -131,6 +131,8 @@ endif()
 set(kep3_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/epoch.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/planet.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/ic2par.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/par2ic.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/detail/type_name.cpp"
 )
 
@@ -178,6 +180,14 @@ endif()
 find_package(spdlog CONFIG REQUIRED)
 target_link_libraries(kep3 PRIVATE spdlog::spdlog)
 
+# xtensor.
+find_package(xtensor CONFIG REQUIRED)
+target_link_libraries(kep3 PRIVATE xtensor)
+
+# xtensor.
+find_package(xtensor-blas CONFIG REQUIRED)
+target_link_libraries(kep3 PRIVATE xtensor-blas)
+
 # Configure config.hpp.
 configure_file("${CMAKE_CURRENT_SOURCE_DIR}/config.hpp.in" "${CMAKE_CURRENT_SOURCE_DIR}/include/kep3/config.hpp" @ONLY)
 

include/kep3/core_astro/convert_anomalies.hpp

@@ -48,8 +48,8 @@ inline double f2e(double f, double e) {
 }
 
 // gudermannian to true (only hyperbolas) e>1 (returns in range [-pi,pi])
-inline double zeta2f(double E, double e) {
-  return 2 * std::atan(std::sqrt((1 + e) / (e - 1)) * std::tan(E / 2));
+inline double zeta2f(double f, double e) {
+  return 2 * std::atan(std::sqrt((1 + e) / (e - 1)) * std::tan(f / 2));
 }
 // true to gudermannian (only hyperbolas) e>1 (returns in range [-pi,pi])
 inline double f2zeta(double zeta, double e) {

include/kep3/core_astro/ic2par.hpp

@@ -0,0 +1,27 @@
+// 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/.
+
+/// From cartesian to osculating Keplerian
+/**
+ * Transforms cartesian coordinates (r,v) to Keplerian elements (a,e,i,W,w,E).
+ * Note that we use the eccentric anomaly (or Gudermannian if e > 1)
+ */
+
+#ifndef kep3_IC2PAR_H
+#define kep3_IC2PAR_H
+
+#include <array>
+
+#include<kep3/detail/visibility.hpp>
+
+namespace kep3 {
+kep3_DLL_PUBLIC std::array<double, 6> ic2par(const std::array<double, 3> &rin,
+                             const std::array<double, 3> &vin, double mu);
+} // namespace kep3
+#endif // kep3_IC2PAR_H

include/kep3/core_astro/par2ic.hpp

@@ -0,0 +1,27 @@
+// 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/.
+
+/// From cartesian to osculating Keplerian
+/**
+ * Transforms cartesian coordinates (r,v) to Keplerian elements (a,e,i,W,w,E).
+ * Note that we use the eccentric anomaly (or Gudermannian if e > 1)
+ */
+
+#ifndef kep3_PAR2IC_H
+#define kep3_PAR2IC_H
+
+#include <array>
+
+#include<kep3/detail/visibility.hpp>
+
+namespace kep3 {
+kep3_DLL_PUBLIC std::array<std::array<double, 3>, 2> par2ic(const std::array<double, 6> &par,
+                                            double mu);
+} // namespace kep3
+#endif // kep3_PAR2IC_H

kep3_devel.yml

@@ -10,4 +10,5 @@ dependencies:
   - fmt
   - heyoka >=0.21.0
   - spdlog
+  - xtensor
 

src/core_astro/ic2par.cpp

@@ -0,0 +1,103 @@
+// 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 <cmath>
+
+#include <boost/math/constants/constants.hpp>
+
+#include <xtensor-blas/xlinalg.hpp>
+#include <xtensor/xadapt.hpp>
+#include <xtensor/xtensor.hpp>
+
+#include <kep3/core_astro/ic2par.hpp>
+
+using boost::math::constants::pi;
+using xt::linalg::cross;
+using xt::linalg::dot;
+
+namespace kep3 {
+
+// r,v,mu -> keplerian osculating elements [a,e,i,W,w,E]. The last
+// is the eccentric anomaly or the Gudermannian according to e. The
+// semi-major axis a is positive for ellipses, negative for hyperbolae.
+// The anomalies W, w and E are in [0, 2pi]. Inclination is in [0, pi].
+
+std::array<double, 6> ic2par(const std::array<double, 3> &rin,
+                             const std::array<double, 3> &vin, double mu) {
+  // Return value
+  std::array<double, 6> retval{};
+  // 0 - We prepare a few xtensor constants.
+  xt::xtensor_fixed<double, xt::xshape<3>> k{0.0, 0.0, 1.0};
+  xt::xtensor_fixed<double, xt::xshape<3>> r0 = xt::adapt(rin);
+  xt::xtensor_fixed<double, xt::xshape<3>> v0 = xt::adapt(vin);
+
+  // 1 - We compute the orbital angular momentum vector
+  auto h = cross(r0, v0); // h = r0 x v0
+
+  // 2 - We compute the orbital parameter
+  auto p = dot(h, h) / mu; // p = h^2 / mu
+
+  // 3 - We compute the vector of the node line
+  // This operation is singular when inclination is zero, in which case the
+  // Keplerian orbital parameters are not well defined
+  auto n = cross(k, h);
+  n = n / xt::linalg::norm(n); // n = (k x h) / |k x h|
+
+  // 4 - We compute the eccentricity vector
+  auto R0 = xt::linalg::norm(r0);
+  auto evett = cross(v0, h) / mu - r0 / R0; // e = (v x h)/mu - r0/R0;
+
+  // The eccentricity is calculated and stored as the second orbital element
+  retval[1] = xt::linalg::norm(evett);
+
+  // The semi-major axis (positive for ellipses, negative for hyperbolas) is
+  // calculated and stored as the first orbital element a = p / (1 - e^2)
+  retval[0] = p(0) / (1 - retval[1] * retval[1]);
+
+  // Inclination is calculated and stored as the third orbital element
+  // i = acos(hy/h)
+  retval[2] = std::acos(h(2) / xt::linalg::norm(h));
+
+  // Argument of pericentrum is calculated and stored as the fifth orbital
+  // elemen.t w = acos(n.e)\|n||e|
+  auto temp = dot(n, evett);
+  retval[4] = std::acos(temp(0) / retval[1]);
+  if (evett(2) < 0) {
+    retval[4] = 2 * pi<double>() - retval[4];
+  }
+  // Argument of longitude is calculated and stored as the fourth orbital
+  // element
+  retval[3] = std::acos(n(0));
+  if (n(1) < 0) {
+    retval[3] = 2 * boost::math::constants::pi<double>() - retval[3];
+  }
+
+  // 4 - We compute ni: the true anomaly (in 0, 2*PI)
+  temp = dot(evett, r0);
+  auto ni = std::acos(temp(0) / retval[1] / R0);
+
+  temp = dot(r0, v0);
+  if (temp(0) < 0.0) {
+    ni = 2 * boost::math::constants::pi<double>() - ni;
+  }
+
+  // Eccentric anomaly or the gudermannian is calculated and stored as the
+  // sixth orbital element
+  if (retval[1] < 1.0) {
+    retval[5] = 2.0 * atan(sqrt((1 - retval[1]) / (1 + retval[1])) *
+                           tan(ni / 2.0)); // algebraic Kepler's equation
+  } else {
+    retval[5] =
+        2.0 * atan(sqrt((retval[1] - 1) / (retval[1] + 1)) *
+                   tan(ni / 2.0)); // algebraic equivalent of Kepler's
+                                   // equation in terms of the Gudermannian
+  }
+  return retval;
+}
+} // namespace kep3

src/core_astro/par2ic.cpp

@@ -0,0 +1,108 @@
+// 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 <cmath>
+
+#include <xtensor-blas/xlinalg.hpp>
+#include <xtensor/xadapt.hpp>
+#include <xtensor/xarray.hpp>
+
+#include <boost/math/constants/constants.hpp>
+
+#include <kep3/core_astro/par2ic.hpp>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+namespace kep3 {
+
+constexpr double pi4{boost::math::constants::quarter_pi<double>()};
+
+// keplerian osculating elements [a,e,i,W,w,E] -> r,v.
+// The last osculating elements needs to be the eccentric anomaly or
+// the Gudermannian according to e. The semi-major axis a needs to be positive
+// for ellipses, negative for hyperbolae.
+// The anomalies W, w and E must be in [0, 2pi] and inclination in [0, pi].
+
+std::array<std::array<double, 3>, 2> par2ic(const std::array<double, 6> &par,
+                                            double mu) {
+  // Return values
+  std::array<double, 3> pos{};
+  std::array<double, 3> vel{};
+  auto pos_xt = xt::adapt(pos, {3u, 1u});
+  auto vel_xt = xt::adapt(vel, {3u, 1u});
+
+  // Rename some variables for readibility
+  double acc = par[0];
+  double ecc = par[1];
+  double inc = par[2];
+  double omg = par[3];
+  double omp = par[4];
+  double EA = par[5];
+
+  // TODO(darioizzo): Check a<0 if e>1
+
+  // 1 - We start by evaluating position and velocity in the perifocal reference
+  // system
+  double xper = 0., yper = 0., xdotper = 0., ydotper = 0.;
+  if (ecc < 1.0) // EA is the eccentric anomaly
+  {
+    double b = acc * std::sqrt(1 - ecc * ecc);
+    double n = std::sqrt(mu / (acc * acc * acc));
+    xper = acc * (std::cos(EA) - ecc);
+    yper = b * std::sin(EA);
+    xdotper = -(acc * n * std::sin(EA)) / (1 - ecc * std::cos(EA));
+    ydotper = (b * n * std::cos(EA)) / (1 - ecc * std::cos(EA));
+  } else // EA is the Gudermannian
+  {
+    double b = -acc * std::sqrt(ecc * ecc - 1);
+    double n = std::sqrt(-mu / (acc * acc * acc));
+
+    double dNdZeta = ecc * (1 + std::tan(EA) * std::tan(EA)) -
+                     (0.5 + 0.5 * std::pow(std::tan(0.5 * EA + pi4), 2)) /
+                         std::tan(0.5 * EA + pi4);
+
+    xper = acc / std::cos(EA) - acc * ecc;
+    yper = b * std::tan(EA);
+
+    xdotper = acc * tan(EA) / cos(EA) * n / dNdZeta;
+    ydotper = b / pow(cos(EA), 2) * n / dNdZeta;
+  }
+
+  // 2 - We then built the rotation matrix from perifocal reference frame to
+  // inertial
+
+  double cosomg = std::cos(omg);
+  double cosomp = std::cos(omp);
+  double sinomg = std::sin(omg);
+  double sinomp = std::sin(omp);
+  double cosi = std::cos(inc);
+  double sini = std::sin(inc);
+
+  xt::xtensor_fixed<double, xt::xshape<3, 3>> R = {
+      {cosomg * cosomp - sinomg * sinomp * cosi,
+       -cosomg * sinomp - sinomg * cosomp * cosi, sinomg * sini},
+      {sinomg * cosomp + cosomg * sinomp * cosi,
+       -sinomg * sinomp + cosomg * cosomp * cosi, -cosomg * sini},
+      {sinomp * sini, cosomp * sini, cosi}};
+
+  // 3 - We end by transforming according to this rotation matrix
+  xt::xtensor_fixed<double, xt::xshape<3, 1>> temp1{{xper}, {yper}, {0.0}};
+  xt::xtensor_fixed<double, xt::xshape<3, 1>> temp2{
+      {xdotper}, {ydotper}, {0.0}};
+
+  pos_xt = R * temp1;
+  vel_xt = R * temp2;
+  fmt::print("\npar2ic: {}, {}", pos_xt, vel_xt);
+  fmt::print("\npar2ic: {}, {}", pos, vel);
+
+  return {pos, vel};
+}
+
+} // namespace kep3

test/CMakeLists.txt

@@ -28,4 +28,5 @@ endfunction()
 
 ADD_kep3_TESTCASE(convert_anomalies_test)
 ADD_kep3_TESTCASE(epoch_test)
-ADD_kep3_TESTCASE(planet_test)
+ADD_kep3_TESTCASE(planet_test)
+ADD_kep3_TESTCASE(ic2par2ic_test)

test/convert_anomalies_test.cpp

@@ -14,7 +14,13 @@
 
 #include "catch.hpp"
 
-using namespace kep3;
+using kep3::e2m;
+using kep3::m2e;
+using kep3::e2f;
+using kep3::f2e;
+using kep3::zeta2f;
+using kep3::f2zeta;
+
 
 TEST_CASE("m2e") {
   using Catch::Detail::Approx;

test/ic2par2ic_test.cpp

@@ -0,0 +1,24 @@
+// 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 <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <kep3/core_astro/ic2par.hpp>
+#include <kep3/core_astro/par2ic.hpp>
+
+#include "catch.hpp"
+
+TEST_CASE("ic2par") {
+    fmt::print("\nResult: {}", kep3::ic2par({1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, 1.0));
+    fmt::print("\nResult: {}", kep3::par2ic({1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 1.0));
+    fmt::print("\nResult: {}", kep3::ic2par({1.0, 0.0, 0.0}, {0.0, 0.0, 1.1}, 1.0));
+    fmt::print("\nResult: {}", kep3::par2ic({1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 1.0));
+
+}