diff --git a/include/heyoka/model/sgp4.hpp b/include/heyoka/model/sgp4.hpp index 955476074..6c64d7793 100644 --- a/include/heyoka/model/sgp4.hpp +++ b/include/heyoka/model/sgp4.hpp @@ -9,6 +9,7 @@ #ifndef HEYOKA_MODEL_SGP4_HPP #define HEYOKA_MODEL_SGP4_HPP +#include #include #include #include @@ -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)); + } + } } } diff --git a/src/model/sgp4.cpp b/src/model/sgp4.cpp index a198c15dd..8022a588a 100644 --- a/src/model/sgp4.cpp +++ b/src/model/sgp4.cpp @@ -587,14 +587,8 @@ T sgp4_date_to_tdelta(SizeType i, Dates dates, const std::vector &sat_buffer, const auto epoch_hi = sat_buffer[static_cast(7) * n_sats + i]; const auto epoch_lo = sat_buffer[static_cast(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)); // Normalise it into a double-length number. const auto epoch = normalise(dfloat(epoch_hi, epoch_lo)); diff --git a/test/model_sgp4.cpp b/test/model_sgp4.cpp index 316205a6e..cc25206e6 100644 --- a/test/model_sgp4.cpp +++ b/test/model_sgp4.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -27,7 +28,7 @@ using namespace heyoka_test; const auto revday2radmin = [](auto x) { return x * 2. * boost::math::constants::pi() / 1440.; }; const auto deg2rad = [](auto x) { return x * 2. * boost::math::constants::pi() / 360.; }; -TEST_CASE("basic") +TEST_CASE("model expression") { detail::edb_disabler ed; @@ -121,7 +122,58 @@ TEST_CASE("basic") } } -TEST_CASE("propagator") +TEST_CASE("propagator basics") +{ + detail::edb_disabler ed; + + using prop_t = model::sgp4_propagator; + + REQUIRE_NOTHROW(prop_t{}); + + // Copy construction. + using md_input_t = mdspan>; + + const std::vector 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 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; @@ -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 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{{{2460486.5 + 1, 0.6478633000000116}, {0., 1.}}}; + prop_t::in_1d date_in{dates.data(), 2}; + + std::vector 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 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{{{2460486.5 + 1, 0.6478633000000116}, + {2458826.5, 0.6933954099999937}, + {2460486.5, 0.6478633000000116}, + {2458826.5 + 1, 0.6933954099999937}}}; + prop_t::in_2d date_b{dates_batch.data(), 1, 2}; + + std::vector 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") @@ -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>; + using prop_t = model::sgp4_propagator; + + std::vector 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>{ins.data(), tot_N}}; + + std::vector times; + times.resize(tot_N, prop_t::date{2460486.5, 0.6478633000000116}); + + std::vector outs; + outs.resize(tot_N * 6); + prop_t::out_2d out_span{outs.data(), 6, tot_N}; + + prop(out_span, prop_t::in_1d{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{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.)); + } + } +}