From a9359f5d8d1fb2a56a015bcb9cf2a9b551d495af Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Sun, 27 Aug 2023 22:29:14 +0200 Subject: [PATCH] propagate lagrangian tests on the way --- include/kep3/core_astro/ic2par.hpp | 6 ++-- src/core_astro/ic2par.cpp | 8 ++--- test/ic2par2ic_test.cpp | 20 ++++++------ test/propagate_lagrangian_test.cpp | 51 +++++++++++++++++++++++++++--- 4 files changed, 63 insertions(+), 22 deletions(-) diff --git a/include/kep3/core_astro/ic2par.hpp b/include/kep3/core_astro/ic2par.hpp index 21f431d0..9d4686b3 100644 --- a/include/kep3/core_astro/ic2par.hpp +++ b/include/kep3/core_astro/ic2par.hpp @@ -18,10 +18,10 @@ #include -#include +#include namespace kep3 { -kep3_DLL_PUBLIC std::array ic2par(const std::array &rin, - const std::array &vin, double mu); +kep3_DLL_PUBLIC std::array +ic2par(const std::array, 2> &pos_vel, double mu); } // namespace kep3 #endif // kep3_IC2PAR_H \ No newline at end of file diff --git a/src/core_astro/ic2par.cpp b/src/core_astro/ic2par.cpp index d3088093..aaa20fe1 100644 --- a/src/core_astro/ic2par.cpp +++ b/src/core_astro/ic2par.cpp @@ -28,14 +28,14 @@ namespace kep3 { // for hyperbolae. The anomalies W, w, f are in [0, 2pi]. Inclination is in [0, // pi]. -std::array ic2par(const std::array &rin, - const std::array &vin, double mu) { +std::array +ic2par(const std::array, 2> &pos_vel, double mu) { // Return value std::array retval{}; // 0 - We prepare a few xtensor constants. xt::xtensor_fixed> k{0.0, 0.0, 1.0}; - auto r0 = xt::adapt(rin); - auto v0 = xt::adapt(vin); + auto r0 = xt::adapt(pos_vel[0]); + auto v0 = xt::adapt(pos_vel[1]); // 1 - We compute the orbital angular momentum vector auto h = cross(r0, v0); // h = r0 x v0 diff --git a/test/ic2par2ic_test.cpp b/test/ic2par2ic_test.cpp index ec080b57..1f55075f 100644 --- a/test/ic2par2ic_test.cpp +++ b/test/ic2par2ic_test.cpp @@ -29,7 +29,7 @@ constexpr double pi{boost::math::constants::pi()}; TEST_CASE("ic2par") { // Singular orbit zero inclination and eccentricity { - auto par = ic2par({1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, 1.0); + auto par = ic2par({1.0, 0.0, 0.0, 0.0, 1.0, 0.0}, 1.0); REQUIRE(par[0] == 1.); // sma is 1 REQUIRE(par[1] == 0.); // ecc is zero REQUIRE(par[2] == 0.); // incl is zero @@ -39,7 +39,7 @@ TEST_CASE("ic2par") { } // Orbit at 90 degrees inclination { - auto par = ic2par({1.0, 0.0, 0.0}, {0.0, 0.0, 1.1}, 1.0); + auto par = ic2par({1.0, 0.0, 0.0, 0.0, 0.0, 1.1}, 1.0); REQUIRE_THAT(par[0], WithinRel(1.2658227848101269, 1e-14)); REQUIRE_THAT(par[1], WithinRel(0.21, 1e-14)); REQUIRE_THAT(par[2], WithinRel(pi / 2, 1e-14)); // incl at 90 degrees @@ -49,7 +49,7 @@ TEST_CASE("ic2par") { } // Orbit at 90 degrees inclination { - auto par = ic2par({1.0, 0.0, 0.0}, {0.0, 0.0, -1.1}, 1.0); + auto par = ic2par({1.0, 0.0, 0.0, 0.0, 0.0, -1.1}, 1.0); REQUIRE_THAT(par[0], WithinRel(1.2658227848101269, 1e-14)); REQUIRE_THAT(par[1], WithinRel(0.21, 1e-14)); REQUIRE_THAT(par[2], WithinRel(pi / 2, 1e-14)); // incl at 90 degrees @@ -59,7 +59,7 @@ TEST_CASE("ic2par") { } // Retrogade orbit { - auto par = ic2par({1.0, 0.0, 0.0}, {0.0, -1.0, 0.1}, 1.0); + auto par = ic2par({1.0, 0.0, 0.0, 0.0, -1.0, 0.1}, 1.0); REQUIRE_THAT(par[0], WithinRel(1.01010101010101, 1e-14)); REQUIRE_THAT(par[1], WithinRel(0.01, 1e-14)); REQUIRE_THAT(par[2], WithinRel(174.289406862500 / 180.0 * pi, @@ -70,8 +70,8 @@ TEST_CASE("ic2par") { } // A random orbit { - auto par = ic2par({-1.1823467354129, 0.0247369349235, -0.014848484784}, - {0.00232349642367, 1.1225625625625, -0.34678634567}, 1.0); + auto par = ic2par({-1.1823467354129, 0.0247369349235, -0.014848484784, + 0.00232349642367, 1.1225625625625, -0.34678634567}, 1.0); REQUIRE_THAT(par[0], WithinRel(3.21921322281178, 1e-14)); REQUIRE_THAT(par[1], WithinAbs(0.63283595179672, 1e-14)); REQUIRE_THAT(par[2], WithinRel(162.82902986040048 / 180.0 * pi, @@ -161,8 +161,8 @@ TEST_CASE("ic2par2ic") { double Omega = Omega_d(rng_engine); double omega = omega_d(rng_engine); double f = f_d(rng_engine); - auto [r, v] = par2ic({sma, ecc, incl, Omega, omega, f}, 1.); - auto par = ic2par(r, v, 1.0); + auto pos_vel = par2ic({sma, ecc, incl, Omega, omega, f}, 1.); + auto par = ic2par(pos_vel, 1.0); REQUIRE_THAT(par[0], WithinAbs(sma, 1e-10)); } } @@ -180,8 +180,8 @@ TEST_CASE("ic2par2ic") { if (std::cos(f) < -1 / ecc) { continue; } - auto [r, v] = par2ic({sma, ecc, incl, Omega, omega, f}, 1.); - auto par = ic2par(r, v, 1.0); + auto pos_vel = par2ic({sma, ecc, incl, Omega, omega, f}, 1.); + auto par = ic2par(pos_vel, 1.0); REQUIRE_THAT(par[0], WithinAbs(sma, 1e-10)); } } diff --git a/test/propagate_lagrangian_test.cpp b/test/propagate_lagrangian_test.cpp index 8825625e..90c9c2a6 100644 --- a/test/propagate_lagrangian_test.cpp +++ b/test/propagate_lagrangian_test.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 #include #include @@ -14,6 +15,7 @@ #include +#include #include #include @@ -26,10 +28,13 @@ using kep3::propagate_lagrangian_u; constexpr double pi{boost::math::constants::pi()}; -TEST_CASE("propagate_lagrangian") { +void test_propagate_lagrangian( + const std::function, 2> &, double, + double)> &propagate, + unsigned int N = 10000) { { // Elliptical circular motion xy std::array, 2> pos_vel = {1., 0, 0., 0., 1., 0.}; - propagate_lagrangian(pos_vel, pi / 2., 1.); + propagate(pos_vel, pi / 2., 1.); auto &[pos, vel] = pos_vel; REQUIRE_THAT(pos[0], WithinAbs(0., 1e-14)); @@ -41,7 +46,7 @@ TEST_CASE("propagate_lagrangian") { } { // Elliptical circular motion xy std::array, 2> pos_vel = {1., 0, 0., 0., 1., 0.}; - propagate_lagrangian(pos_vel, - pi / 2., 1.); + propagate(pos_vel, -pi / 2., 1.); auto &[pos, vel] = pos_vel; REQUIRE_THAT(pos[0], WithinAbs(0., 1e-14)); @@ -53,7 +58,7 @@ TEST_CASE("propagate_lagrangian") { } { // Elliptical circular motion xz std::array, 2> pos_vel = {1., 0, 0., 0., 0., 1.}; - propagate_lagrangian(pos_vel, pi / 2., 1.); + propagate(pos_vel, pi / 2., 1.); auto &[pos, vel] = pos_vel; REQUIRE_THAT(pos[0], WithinAbs(0., 1e-14)); @@ -65,7 +70,7 @@ TEST_CASE("propagate_lagrangian") { } { // Elliptical circular motion yz std::array, 2> pos_vel = {0., 1, 0., 0., 0., 1.}; - propagate_lagrangian(pos_vel, pi / 2., 1.); + propagate(pos_vel, pi / 2., 1.); auto &[pos, vel] = pos_vel; REQUIRE_THAT(pos[0], WithinAbs(0., 1e-14)); @@ -75,4 +80,40 @@ TEST_CASE("propagate_lagrangian") { REQUIRE_THAT(vel[1], WithinAbs(-1., 1e-14)); REQUIRE_THAT(vel[2], WithinAbs(0., 1e-14)); } + // We test orbital parameters are unchanged for random propagations + // Engines + // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp) + std::mt19937 rng_engine(12201203u); + // Distribution for the initial Cartesian state + std::uniform_real_distribution r_d(-2., 2.); + std::uniform_real_distribution v_d(-2., 2.); + std::uniform_real_distribution time(0.1, 20.); + // Testing on N random calls + for (auto i = 0u; i < N; ++i) { + std::array pos = {r_d(rng_engine), r_d(rng_engine), + r_d(rng_engine)}; + std::array vel = {v_d(rng_engine), v_d(rng_engine), + v_d(rng_engine)}; + std::array, 2> pos_vel = {pos, vel}; + auto par_before = kep3::ic2par(pos_vel, 1.0); + // We filter out cases of little significance (too close to singularity) + if (std::abs(par_before[0]) > 0.5 && std::abs(par_before[0]) < 10. && + std::abs(1 - par_before[1]) > 1e-1) { + propagate(pos_vel, time(rng_engine), 1.); + auto par_after = kep3::ic2par(pos_vel, 1.0); + fmt::print("\n\n{}", par_before); + fmt::print("\n{}", par_after); + + REQUIRE_THAT(par_before[0], WithinAbs(par_after[0], 1e-11)); + REQUIRE_THAT(par_before[1], WithinAbs(par_after[1], 1e-11)); + REQUIRE_THAT(par_before[2], WithinAbs(par_after[2], 1e-11)); + REQUIRE_THAT(par_before[3], WithinAbs(par_after[3], 1e-11)); + REQUIRE_THAT(par_before[4], WithinAbs(par_after[4], 1e-11)); + } + } +} + +TEST_CASE("propagate_lagrangian") { + //test_propagate_lagrangian(&propagate_lagrangian, 10000u); + // test_propagate_lagrangian(&propagate_lagrangian_u, 10000u); } \ No newline at end of file