Skip to content

Commit

Permalink
propagate lagrangian tests on the way
Browse files Browse the repository at this point in the history
  • Loading branch information
@darioizzo
darioizzo committed Aug 27, 2023
1 parent 8a7f6bc commit a9359f5
Show file tree
Hide file tree
Showing 4 changed files with 63 additions and 22 deletions.
6 changes: 3 additions & 3 deletions include/kep3/core_astro/ic2par.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,10 +18,10 @@

#include <array>

#include<kep3/detail/visibility.hpp>
#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);
kep3_DLL_PUBLIC std::array<double, 6>
ic2par(const std::array<std::array<double, 3>, 2> &pos_vel, double mu);
} // namespace kep3
#endif // kep3_IC2PAR_H
8 changes: 4 additions & 4 deletions src/core_astro/ic2par.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,14 +28,14 @@ namespace kep3 {
// for hyperbolae. The anomalies W, w, f 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) {
std::array<double, 6>
ic2par(const std::array<std::array<double, 3>, 2> &pos_vel, 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};
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
Expand Down
20 changes: 10 additions & 10 deletions test/ic2par2ic_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ constexpr double pi{boost::math::constants::pi<double>()};
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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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));
}
}
Expand All @@ -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));
}
}
Expand Down
51 changes: 46 additions & 5 deletions test/propagate_lagrangian_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,15 @@
// 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 <functional>
#include <stdexcept>

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

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

#include <kep3/core_astro/ic2par.hpp>
#include <kep3/core_astro/kepler_equations.hpp>
#include <kep3/core_astro/propagate_lagrangian.hpp>

Expand All @@ -26,10 +28,13 @@ using kep3::propagate_lagrangian_u;

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

TEST_CASE("propagate_lagrangian") {
void test_propagate_lagrangian(
const std::function<void(std::array<std::array<double, 3>, 2> &, double,
double)> &propagate,
unsigned int N = 10000) {
{ // Elliptical circular motion xy
std::array<std::array<double, 3>, 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));
Expand All @@ -41,7 +46,7 @@ TEST_CASE("propagate_lagrangian") {
}
{ // Elliptical circular motion xy
std::array<std::array<double, 3>, 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));
Expand All @@ -53,7 +58,7 @@ TEST_CASE("propagate_lagrangian") {
}
{ // Elliptical circular motion xz
std::array<std::array<double, 3>, 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));
Expand All @@ -65,7 +70,7 @@ TEST_CASE("propagate_lagrangian") {
}
{ // Elliptical circular motion yz
std::array<std::array<double, 3>, 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));
Expand All @@ -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<double> r_d(-2., 2.);
std::uniform_real_distribution<double> v_d(-2., 2.);
std::uniform_real_distribution<double> time(0.1, 20.);
// Testing on N random calls
for (auto i = 0u; i < N; ++i) {
std::array<double, 3> pos = {r_d(rng_engine), r_d(rng_engine),
r_d(rng_engine)};
std::array<double, 3> vel = {v_d(rng_engine), v_d(rng_engine),
v_d(rng_engine)};
std::array<std::array<double, 3>, 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);
}

0 comments on commit a9359f5

Please sign in to comment.