Skip to content

Commit

Permalink
More testing, coverage.
Browse files Browse the repository at this point in the history
  • Loading branch information
@bluescarni
bluescarni committed Jul 4, 2024
1 parent 667b5d9 commit 426c1ce
Show file tree
Hide file tree
Showing 3 changed files with 216 additions and 10 deletions.
17 changes: 17 additions & 0 deletions include/heyoka/model/sgp4.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#ifndef HEYOKA_MODEL_SGP4_HPP
#define HEYOKA_MODEL_SGP4_HPP

#include <cmath>
#include <concepts>
#include <cstddef>
#include <cstdint>
Expand Down Expand Up @@ -83,6 +84,22 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS sgp4_propagator
for (auto i = 0u; i < 9u; ++i) {
for (std::size_t j = 0; j < in.extent(1); ++j) {
sat_buffer.push_back(in(i, j));

if (i == 7u) {
using std::abs;

const auto epoch_hi = in(7, j);
const auto epoch_lo = in(8, j);

// NOTE: the magnitude of the high half of the epoch cannot be less than the magnitude
// of the low half in order to use double-length arithmetic.
if (!(abs(epoch_hi) >= abs(epoch_lo))) [[unlikely]] {
throw std::invalid_argument(fmt::format(
"Invalid reference epoch detected for the satellite at index {}: the magnitude "
"of the Julian date ({}) is less than the magnitude of the fractional correction ({})",
j, epoch_hi, epoch_lo));
}
}
}
}

Expand Down
10 changes: 2 additions & 8 deletions src/model/sgp4.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -587,14 +587,8 @@ T sgp4_date_to_tdelta(SizeType i, Dates dates, const std::vector<T> &sat_buffer,
const auto epoch_hi = sat_buffer[static_cast<SizeType>(7) * n_sats + i];
const auto epoch_lo = sat_buffer[static_cast<SizeType>(8) * n_sats + i];

// NOTE: the magnitude of the high half cannot be less than the magnitude
// of the low half in order to use double-length arithmetic.
if (!(abs(epoch_hi) >= abs(epoch_lo))) [[unlikely]] {
throw std::invalid_argument(
fmt::format("Invalid reference epoch detected for the satellite at index {}: the magnitude of the Julian "
"date ({}) is less than the magnitude of the fractional correction ({})",
i, epoch_hi, epoch_lo));
}
// NOTE: this has been checked during construction.
assert(abs(epoch_hi) >= abs(epoch_lo));

Check warning on line 591 in src/model/sgp4.cpp

View check run for this annotation

Codecov / codecov/patch

src/model/sgp4.cpp#L591 

Added line #L591 was not covered by tests

// Normalise it into a double-length number.
const auto epoch = normalise(dfloat(epoch_hi, epoch_lo));
Expand Down
199 changes: 197 additions & 2 deletions test/model_sgp4.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include <cstddef>
#include <ranges>
#include <stdexcept>
#include <utility>
#include <vector>

#include <boost/math/constants/constants.hpp>
Expand All @@ -27,7 +28,7 @@ using namespace heyoka_test;
const auto revday2radmin = [](auto x) { return x * 2. * boost::math::constants::pi<double>() / 1440.; };
const auto deg2rad = [](auto x) { return x * 2. * boost::math::constants::pi<double>() / 360.; };

TEST_CASE("basic")
TEST_CASE("model expression")
{
detail::edb_disabler ed;

Expand Down Expand Up @@ -121,7 +122,58 @@ TEST_CASE("basic")
}
}

TEST_CASE("propagator")
TEST_CASE("propagator basics")
{
detail::edb_disabler ed;

using prop_t = model::sgp4_propagator<double>;

REQUIRE_NOTHROW(prop_t{});

// Copy construction.
using md_input_t = mdspan<const double, extents<std::size_t, 9, std::dynamic_extent>>;

const std::vector<double> ins = {revday2radmin(13.75091047972192),
revday2radmin(15.50103472202482),
0.0024963,
0.0007417,
deg2rad(90.2039),
deg2rad(51.6439),
deg2rad(55.5633),
deg2rad(211.2001),
deg2rad(320.5956),
deg2rad(17.6667),
deg2rad(91.4738),
deg2rad(85.6398),
0.75863e-3,
.38792e-4,
2460486.5,
2458826.5,
0.6478633000000116,
0.6933954099999937};

const auto tm = std::array{1440., 0.};
const prop_t::in_1d<double> tm_in{tm.data(), 2};

prop_t prop{md_input_t{ins.data(), 2}};
auto prop2 = prop;
REQUIRE(prop2.get_n_sats() == 2u);

// Move construction.
auto prop3 = std::move(prop2);
REQUIRE(prop3.get_n_sats() == 2u);

// Revive prop2 via copy assignment.
prop2 = prop3;
REQUIRE(prop2.get_n_sats() == 2u);

// Revive via move assignment.
prop_t prop4;
prop4 = std::move(prop2);
REQUIRE(prop4.get_n_sats() == 2u);
}

TEST_CASE("propagator single")
{
detail::edb_disabler ed;

Expand Down Expand Up @@ -307,6 +359,85 @@ TEST_CASE("error handling")

REQUIRE_THROWS_MATCHES((prop_t{md_input_t{input.data(), 0}}), std::invalid_argument,
Message("Cannot initialise an sgp4_propagator with an empty list of satellites"));

// Invalid jd+frac for the epochs.
std::vector<double> ins = {revday2radmin(13.75091047972192),
revday2radmin(15.50103472202482),
0.0024963,
0.0007417,
deg2rad(90.2039),
deg2rad(51.6439),
deg2rad(55.5633),
deg2rad(211.2001),
deg2rad(320.5956),
deg2rad(17.6667),
deg2rad(91.4738),
deg2rad(85.6398),
0.75863e-3,
.38792e-4,
2460486.5,
0.,
0.6478633000000116,
1.};

REQUIRE_THROWS_MATCHES(
(prop_t{md_input_t{ins.data(), 2}}), std::invalid_argument,
Message("Invalid reference epoch detected for the satellite at index 1: the magnitude "
"of the Julian date (0) is less than the magnitude of the fractional correction (1)"));

ins = std::vector{revday2radmin(13.75091047972192),
revday2radmin(15.50103472202482),
0.0024963,
0.0007417,
deg2rad(90.2039),
deg2rad(51.6439),
deg2rad(55.5633),
deg2rad(211.2001),
deg2rad(320.5956),
deg2rad(17.6667),
deg2rad(91.4738),
deg2rad(85.6398),
0.75863e-3,
.38792e-4,
2460486.5,
2458826.5,
0.6478633000000116,
0.6933954099999937};

prop_t prop{md_input_t{ins.data(), 2}};

auto dates = std::array<prop_t::date, 2>{{{2460486.5 + 1, 0.6478633000000116}, {0., 1.}}};
prop_t::in_1d<prop_t::date> date_in{dates.data(), 2};

std::vector<double> outs(12u);
prop_t::out_2d out{outs.data(), 6, 2};

REQUIRE_THROWS_MATCHES(
prop(out, date_in), std::invalid_argument,
Message("Invalid propagation date detected for the satellite at index 1: the magnitude of the Julian "
"date (0) is less than the magnitude of the fractional correction (1)"));

prop_t::in_1d<prop_t::date> date_in2{dates.data(), 1};

REQUIRE_THROWS_MATCHES(
prop(out, date_in2), std::invalid_argument,
Message("Invalid array of dates passed to the call operator of an sgp4_propagator: the number of "
"satellites is 2, while the number of dates is 1"));

auto dates_batch = std::array<prop_t::date, 4>{{{2460486.5 + 1, 0.6478633000000116},
{2458826.5, 0.6933954099999937},
{2460486.5, 0.6478633000000116},
{2458826.5 + 1, 0.6933954099999937}}};
prop_t::in_2d<prop_t::date> date_b{dates_batch.data(), 1, 2};

std::vector<double> outs_batch(24u);
prop_t::out_3d out_batch{outs.data(), 2, 6, 2};

REQUIRE_THROWS_MATCHES(
prop(out_batch, date_b), std::invalid_argument,
Message("Invalid dimensions detected in batch-mode sgp4 propagation: the number of evaluations "
"inferred from the output array is 2, which is not consistent with the number of evaluations "
"inferred from the times array (1)"));
}

TEST_CASE("derivatives")
Expand Down Expand Up @@ -371,3 +502,67 @@ TEST_CASE("derivatives")
}
}
}

// A test with several satellites to test parallelisation.
TEST_CASE("large")
{
const auto tot_N = 1000;

detail::edb_disabler ed;

using md_input_t = mdspan<double, extents<std::size_t, 9, std::dynamic_extent>>;
using prop_t = model::sgp4_propagator<double>;

std::vector<double> ins;
ins.resize(tot_N * 9);
md_input_t in(ins.data(), tot_N);
for (auto i = 0; i < tot_N; ++i) {
in(0, i) = revday2radmin(13.75091047972192);
in(1, i) = 0.0024963;
in(2, i) = deg2rad(90.2039);
in(3, i) = deg2rad(55.5633);
in(4, i) = deg2rad(320.5956);
in(5, i) = deg2rad(91.4738);
in(6, i) = 0.75863e-3;
in(7, i) = 2460486.5;
in(8, i) = 0.6478633000000116;
}

prop_t prop{mdspan<const double, extents<std::size_t, 9, std::dynamic_extent>>{ins.data(), tot_N}};

std::vector<prop_t::date> times;
times.resize(tot_N, prop_t::date{2460486.5, 0.6478633000000116});

std::vector<double> outs;
outs.resize(tot_N * 6);
prop_t::out_2d out_span{outs.data(), 6, tot_N};

prop(out_span, prop_t::in_1d<prop_t::date>{times.data(), tot_N});

for (auto i = 0; i < tot_N; ++i) {
REQUIRE(out_span(0, i) == approximately(2561.223660636298, 10000.));
REQUIRE(out_span(1, i) == approximately(3698.797144057697, 10000.));
REQUIRE(out_span(2, i) == approximately(5818.772215708888, 10000.));
REQUIRE(out_span(3, i) == approximately(-3.276142513618007, 10000.));
REQUIRE(out_span(4, i) == approximately(-4.806489082829041, 10000.));
REQUIRE(out_span(5, i) == approximately(4.511134501638151, 10000.));
}

const auto n_evals = 100;
times.resize(tot_N * n_evals, prop_t::date{2460486.5, 0.6478633000000116});
outs.resize(tot_N * 6 * n_evals);
prop_t::out_3d out_span_batch{outs.data(), n_evals, 6, tot_N};

prop(out_span_batch, prop_t::in_2d<prop_t::date>{times.data(), n_evals, tot_N});

for (auto i = 0; i < tot_N; ++i) {
for (auto k = 0; k < n_evals; ++k) {
REQUIRE(out_span_batch(k, 0, i) == approximately(2561.223660636298, 10000.));
REQUIRE(out_span_batch(k, 1, i) == approximately(3698.797144057697, 10000.));
REQUIRE(out_span_batch(k, 2, i) == approximately(5818.772215708888, 10000.));
REQUIRE(out_span_batch(k, 3, i) == approximately(-3.276142513618007, 10000.));
REQUIRE(out_span_batch(k, 4, i) == approximately(-4.806489082829041, 10000.));
REQUIRE(out_span_batch(k, 5, i) == approximately(4.511134501638151, 10000.));
}
}
}

0 comments on commit 426c1ce

Please sign in to comment.