Skip to content

Commit

Permalink
eq2ic etc.. tests finished
Browse files Browse the repository at this point in the history
  • Loading branch information
darioizzo committed Aug 31, 2023
1 parent e262488 commit 7422a43
Show file tree
Hide file tree
Showing 8 changed files with 321 additions and 28 deletions.
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ set(kep3_SRC_FILES
"${CMAKE_CURRENT_SOURCE_DIR}/src/planet.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"
"${CMAKE_CURRENT_SOURCE_DIR}/src/core_astro/propagate_lagrangian.cpp"
"${CMAKE_CURRENT_SOURCE_DIR}/src/detail/type_name.cpp"
)
Expand Down
6 changes: 3 additions & 3 deletions include/kep3/core_astro/eq2par2eq.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@
* Note that we use the eccentric anomaly (or Gudermannian if e > 1)
*/

#ifndef kep3_EQ2PAR_H
#define kep3_EQ2PAR_H
#ifndef kep3_EQ2PAR2EQ_H
#define kep3_EQ2PAR2EQ_H

#include <array>

Expand All @@ -29,4 +29,4 @@ kep3_DLL_PUBLIC std::array<double, 6> par2eq(const std::array<double, 6> &par,
bool retrogade = false);

} // namespace kep3
#endif // kep3_EQ2PAR_H
#endif // kep3_EQ2PAR2EQ_H
36 changes: 23 additions & 13 deletions src/core_astro/eq2par2eq.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,33 +12,44 @@

#include <boost/math/constants/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 = false) {
std::array<double, 6> eq2par(const std::array<double, 6> &eq, bool retrogade) {
std::array<double, 6> retval{};
int I = 1;
if (retrogade) {
I = -1;
}
auto ecc = std::sqrt(eq[1] * eq[1] + eq[2] * eq[2]);
auto tmp = std::sqrt(eq[3] * eq[3] + eq[4] * eq[4]);
auto zita = std::atan2(eq[2] / ecc, eq[1] / ecc);
double ecc = std::sqrt(eq[1] * eq[1] + eq[2] * eq[2]);
double tmp = std::sqrt(eq[3] * eq[3] + eq[4] * eq[4]);
double zita = std::atan2(eq[2] / ecc, eq[1] / ecc); // [-pi, pi]
if (zita < 0) {
zita += 2 * pi; // [0, 2*pi]
}

retval[1] = ecc;
retval[0] = eq[0] / (1. - ecc * ecc);
retval[2] = half_pi * (1. - I) +
2. * I * std::atan(tmp);
retval[3] = std::atan2(eq[4] / tmp, eq[3] / tmp);
retval[4] = zita - I * retval[3];
retval[5] = eq[5] - zita;
retval[2] = half_pi * (1. - I) + 2. * I * std::atan(tmp);
retval[3] = std::atan2(eq[4] / tmp, eq[3] / tmp); // [-pi, pi]
if (retval[3] < 0) {
retval[3] += 2 * pi; // [0, 2*pi]
}
retval[4] = zita - I * retval[3]; //
if (retval[4] < 0) {
retval[4] += 2 * pi;
} else if (retval[4] > 2 * pi) {
retval[4] -= 2 * pi;
}
retval[5] = eq[5] - I * retval[3] - retval[4];
return retval;
}

std::array<double, 6> par2eq(const std::array<double, 6> &par,
bool retrogade = false) {
std::array<double, 6> par2eq(const std::array<double, 6> &par, bool retrogade) {
std::array<double, 6> eq{};
int I = 0;
if (retrogade) {
Expand All @@ -56,5 +67,4 @@ std::array<double, 6> par2eq(const std::array<double, 6> &par,
eq[5] = par[5] + par[4] + I * par[3];
return eq;
}

} // namespace kep3
1 change: 1 addition & 0 deletions test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -31,4 +31,5 @@ ADD_kep3_TESTCASE(epoch_test)
ADD_kep3_TESTCASE(planet_test)
ADD_kep3_TESTCASE(ic2par2ic_test)
ADD_kep3_TESTCASE(ic2eq2ic_test)
ADD_kep3_TESTCASE(eq2par2eq_test)
ADD_kep3_TESTCASE(propagate_lagrangian_test)
174 changes: 174 additions & 0 deletions test/eq2par2eq_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
// 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 <stdexcept>

#include <fmt/core.h>
#include <fmt/ranges.h>

#include <boost/math/constants/constants.hpp>

#include <kep3/core_astro/eq2par2eq.hpp>

#include "catch.hpp"
#include "test_helpers.hpp"

constexpr double pi{boost::math::constants::pi<double>()};

TEST_CASE("eq2par2eq") {
// Engines
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
// Distributions for the elements
std::uniform_real_distribution<double> sma_d(1.1, 100.);
std::uniform_real_distribution<double> ecc_d(0, 0.99);
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> ni_d(0, 2 * pi);

{
// Testing on N random calls on ellipses
unsigned N = 10000;
for (auto i = 0u; i < N; ++i) {
// We sample randomly on the Keplerian space
double sma = sma_d(rng_engine);
double ecc = ecc_d(rng_engine);
double incl = incl_d(rng_engine);
double Omega = Omega_d(rng_engine);
double omega = omega_d(rng_engine);
double ni = ni_d(rng_engine);
// Compute the initial eq
auto eq = kep3::par2eq({sma, ecc, incl, Omega, omega, ni});
// Test eq2par2eq
auto par = kep3::eq2par(eq);

// Here we do not use catch matchers to test floating point as for small
// numbers (<1) we care about absolute while for large (>1) we care for
// relative error.
REQUIRE(kep3_tests::floating_point_error(sma, par[0]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(ecc, par[1]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(incl, par[2]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(Omega, par[3]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(omega, par[4]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(ni, par[5]) < 1e-13);
}
}
{
// Testing on N random calls on hyperbolas
unsigned N = 10000;
for (auto i = 0u; i < N; ++i) {
// We sample randomly on the Keplerian space
double sma = -sma_d(rng_engine);
double ecc = ecc_d(rng_engine) + 1.;
double incl = incl_d(rng_engine);
double Omega = Omega_d(rng_engine);
double omega = omega_d(rng_engine);
double ni = ni_d(rng_engine);
// Skipping if true anomaly is way out of asymptotes
if (std::cos(ni) < -1 / ecc + 0.1) {
continue;
}
// Compute the initial eq
auto eq = kep3::par2eq({sma, ecc, incl, Omega, omega, ni});
// Test eq2par2eq
auto par = kep3::eq2par(eq);

// Here we do not use catch matchers to test floating point as for small
// numbers (<1) we care about absolute while for large (>1) we care for
// relative error.
REQUIRE(kep3_tests::floating_point_error(sma, par[0]) <
1e-10); // errors arise since p = a * (1-e^2) and a = p / (1-e^2)
// [when e is close to 1 and a is high]
REQUIRE(kep3_tests::floating_point_error(ecc, par[1]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(incl, par[2]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(Omega, par[3]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(omega, par[4]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(ni, par[5]) < 1e-13);
}
}
}

TEST_CASE("eq2par2eq_retrogade") {
// Engines
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
// Distributions for the elements
std::uniform_real_distribution<double> sma_d(1.1, 100.);
std::uniform_real_distribution<double> ecc_d(0, 0.99);
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> ni_d(0, 2 * pi);

{
// Testing on N random calls on ellipses
unsigned N = 10000;
for (auto i = 0u; i < N; ++i) {
// We sample randomly on the Keplerian space
double sma = sma_d(rng_engine);
double ecc = ecc_d(rng_engine);
double incl = incl_d(rng_engine);
double Omega = Omega_d(rng_engine);
double omega = omega_d(rng_engine);
double ni = ni_d(rng_engine);
// Compute the initial eq
auto eq = kep3::par2eq({sma, ecc, incl, Omega, omega, ni}, true);
// Test eq2par2eq
auto par = kep3::eq2par(eq, true);

// Here we do not use catch matchers to test floating point as for small
// numbers (<1) we care about absolute while for large (>1) we care for
// relative error.
REQUIRE(kep3_tests::floating_point_error(sma, par[0]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(ecc, par[1]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(incl, par[2]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(Omega, par[3]) < 1e-13);
if (kep3_tests::floating_point_error(omega, par[4]) > 1e-13) {
fmt::print("\n{}, {}", omega / pi * 180, par[4] / pi * 180);
}
REQUIRE(kep3_tests::floating_point_error(omega, par[4]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(ni, par[5]) < 1e-13);
}
}
{
// Testing on N random calls on hyperbolas
unsigned N = 10000;
for (auto i = 0u; i < N; ++i) {
// We sample randomly on the Keplerian space
double sma = -sma_d(rng_engine);
double ecc = ecc_d(rng_engine) + 1.;
double incl = incl_d(rng_engine);
double Omega = Omega_d(rng_engine);
double omega = omega_d(rng_engine);
double ni = ni_d(rng_engine);
// Skipping if true anomaly is way out of asymptotes
if (std::cos(ni) < -1 / ecc + 0.1) {
continue;
}
// Compute the initial eq
auto eq = kep3::par2eq({sma, ecc, incl, Omega, omega, ni}, true);
// Test eq2par2eq
auto par = kep3::eq2par(eq, true);

// Here we do not use catch matchers to test floating point as for small
// numbers (<1) we care about absolute while for large (>1) we care for
// relative error.
REQUIRE(kep3_tests::floating_point_error(sma, par[0]) <
1e-10); // errors arise since p = a * (1-e^2) and a = p / (1-e^2)
// [when e is close to 1 and a is high]
REQUIRE(kep3_tests::floating_point_error(ecc, par[1]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(incl, par[2]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(Omega, par[3]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(omega, par[4]) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(ni, par[5]) < 1e-13);
}
}
}
90 changes: 89 additions & 1 deletion test/ic2eq2ic_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include <kep3/core_astro/ic2par2ic.hpp>

#include "catch.hpp"
#include "test_helpers.hpp"

using Catch::Matchers::WithinRel;
using kep3::eq2ic;
Expand Down Expand Up @@ -68,7 +69,7 @@ TEST_CASE("ic2eq2ic") {
// Distributions for the elements
std::uniform_real_distribution<double> sma_d(1.1, 100.);
std::uniform_real_distribution<double> ecc_d(0, 0.99);
std::uniform_real_distribution<double> incl_d(0., 3.); // excluding 180 degrees
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., pi);
std::uniform_real_distribution<double> ni_d(0, 2 * pi);
Expand Down Expand Up @@ -151,3 +152,90 @@ TEST_CASE("ic2eq2ic") {
}
}
}

TEST_CASE("ic2eq2ic_retrogade") {
// Engines
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
// Distributions for the elements
std::uniform_real_distribution<double> sma_d(1.1, 100.);
std::uniform_real_distribution<double> ecc_d(0, 0.99);
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., pi);
std::uniform_real_distribution<double> ni_d(0, 2 * pi);

{
// Testing on N random calls on ellipses
unsigned N = 10000;
for (auto i = 0u; i < N; ++i) {
// We sample randomly on the Keplerian space
double sma = sma_d(rng_engine);
double ecc = ecc_d(rng_engine);
double incl = incl_d(rng_engine);
double Omega = Omega_d(rng_engine);
double omega = omega_d(rng_engine);
double ni = ni_d(rng_engine);
// Compute the initial r,v
auto pos_vel = kep3::par2ic({sma, ecc, incl, Omega, omega, ni}, 1.);
// Test ic2eq2ic
auto eq = ic2eq(pos_vel, 1., true);
auto pos_vel_new = eq2ic(eq, 1.0, true);
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 V = std::sqrt(pos_vel[1][0] * pos_vel[1][0] +
pos_vel[1][1] * pos_vel[1][1] +
pos_vel[1][2] * pos_vel[1][2]);
double R_new = std::sqrt(pos_vel_new[0][0] * pos_vel_new[0][0] +
pos_vel_new[0][1] * pos_vel_new[0][1] +
pos_vel_new[0][2] * pos_vel_new[0][2]);
double V_new = std::sqrt(pos_vel_new[1][0] * pos_vel_new[1][0] +
pos_vel_new[1][1] * pos_vel_new[1][1] +
pos_vel_new[1][2] * pos_vel_new[1][2]);
// Here we do not use catch matchers to test floating point as for small numbers (<1) we care about absolute
// while for large (>1) we care for relative error.
REQUIRE(kep3_tests::floating_point_error(R, R_new) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(V, V_new) < 1e-13);
}
}
{
// Testing on N random calls on ellipses
unsigned N = 10000;
for (auto i = 0u; i < N; ++i) {
// We sample randomly on the Keplerian space
double sma = -sma_d(rng_engine);
double ecc = ecc_d(rng_engine) + 1.1;
double incl = incl_d(rng_engine);
double Omega = Omega_d(rng_engine);
double omega = omega_d(rng_engine);
double ni = ni_d(rng_engine);
// Skipping if true anomaly is way out of asymptotes
if (std::cos(ni) < -1 / ecc + 0.1) {
continue;
}
// Compute the initial r,v
auto pos_vel = kep3::par2ic({sma, ecc, incl, Omega, omega, ni}, 1.);
// Test ic2eq2ic
auto eq = ic2eq(pos_vel, 1., true);
auto pos_vel_new = eq2ic(eq, 1.0, true);

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 V = std::sqrt(pos_vel[1][0] * pos_vel[1][0] +
pos_vel[1][1] * pos_vel[1][1] +
pos_vel[1][2] * pos_vel[1][2]);
double R_new = std::sqrt(pos_vel_new[0][0] * pos_vel_new[0][0] +
pos_vel_new[0][1] * pos_vel_new[0][1] +
pos_vel_new[0][2] * pos_vel_new[0][2]);
double V_new = std::sqrt(pos_vel_new[1][0] * pos_vel_new[1][0] +
pos_vel_new[1][1] * pos_vel_new[1][1] +
pos_vel_new[1][2] * pos_vel_new[1][2]);
// Here we do not use catch matchers to test floating point as for small numbers (<1) we care about absolute
// while for large (>1) we care for relative error.
REQUIRE(kep3_tests::floating_point_error(R, R_new) < 1e-13);
REQUIRE(kep3_tests::floating_point_error(V, V_new) < 1e-13);
}
}
}
Loading

0 comments on commit 7422a43

Please sign in to comment.