From 26b212489961b266cd34e05f51c048f1fb42105d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 1 Jan 2025 12:00:51 +0100 Subject: [PATCH] Implement the UTC -> TAI conversion in the sgp4 propagator. --- include/heyoka/model/sgp4.hpp | 5 ++- src/model/sgp4.cpp | 60 ++++++++++++++++++++++------------- 2 files changed, 40 insertions(+), 25 deletions(-) diff --git a/include/heyoka/model/sgp4.hpp b/include/heyoka/model/sgp4.hpp index e031727d5..928c2b341 100644 --- a/include/heyoka/model/sgp4.hpp +++ b/include/heyoka/model/sgp4.hpp @@ -163,9 +163,8 @@ class HEYOKA_DLL_PUBLIC_INLINE_CLASS sgp4_propagator // - the reference epoch (as a Julian date), // - a fractional correction to the epoch (in Julian days). // - // Note that UTC Julian dates will result in slightly incorrect results when - // propagating across leap seconds. See the Python tutorial for an explanation. - // TAI Julian dates can be used instead for accurate propagation across leap seconds. + // Julian dates are to be provided in the UTC scale of time. Internal conversion + // to TAI will ensure correct propagation across leap seconds. template requires(!igor::has_unnamed_arguments()) explicit sgp4_propagator( diff --git a/src/model/sgp4.cpp b/src/model/sgp4.cpp index 8d753895a..e38130cfb 100644 --- a/src/model/sgp4.cpp +++ b/src/model/sgp4.cpp @@ -851,33 +851,51 @@ namespace detail namespace { +// The erfa function for UTC to TAI conversion. +extern "C" int eraUtctai(double, double, double *, double *); + // Helper to convert a Julian date into a time delta suitable for use in the SGP4 algorithm. // The reference epochs for all satellites are stored in sat_buffer, the dates are passed in // via the 'dates' argument, n_sats is the total number of satellites and i is the index -// of the satellite on which we are operating. The conversion is done in double-length arithmetic, -// using the fractional parts of the date and epochs as lower halves of the double-length -// numbers. +// of the satellite on which we are operating. The double-length UTC Julian dates are first converted +// to TAI Julian dates using an erfa routine. The resulting TAI dates are then normalised and +// subtracted to yield the final time delta. template T sgp4_date_to_tdelta(SizeType i, Dates dates, const std::vector &sat_buffer, std::uint32_t n_sats) { using dfloat = heyoka::detail::dfloat; - // Normalise the propagation date into a double-length number. + // Convert the target propagation date to TAI using erfa. + // NOTE: the docs of the erfa routine indicate that we do not have + // to normalise the input UTC double-length Julian dates. + double tai_prop_date_hi{}, tai_prop_date_lo{}; + auto res = eraUtctai(dates(i).jd, dates(i).frac, &tai_prop_date_hi, &tai_prop_date_lo); + if (res == -1) [[unlikely]] { + throw std::invalid_argument( + fmt::format("An invalid UTC propagation date ({}, {}) was detected for the satellite at index {}", + dates(i).jd, dates(i).frac, i)); + } + + // Convert the satellite epoch to TAI. + double tai_epoch_hi{}, tai_epoch_lo{}; + const auto utc_epoch_hi = sat_buffer[static_cast(7) * n_sats + i]; + const auto utc_epoch_lo = sat_buffer[static_cast(8) * n_sats + i]; + res = eraUtctai(utc_epoch_hi, utc_epoch_lo, &tai_epoch_hi, &tai_epoch_lo); + if (res == -1) [[unlikely]] { + throw std::invalid_argument(fmt::format( + "An invalid UTC epoch ({}, {}) was detected for the satellite at index {}", utc_epoch_hi, utc_epoch_lo, i)); + } + + // Normalise the TAI propagation date into a double-length number. // NOTE: we use Knuth's EFT here in order to ensure correctness // even if the fractional part of the date is greater in magnitude // than the main part. - const auto norm_prop_date = heyoka::detail::eft_add_knuth(dates(i).jd, dates(i).frac); + const auto norm_prop_date + = heyoka::detail::eft_add_knuth(static_cast(tai_prop_date_hi), static_cast(tai_prop_date_lo)); const auto date = dfloat(norm_prop_date.first, norm_prop_date.second); - // Load the reference epoch for the i-th satellite. - const auto epoch_hi = sat_buffer[static_cast(7) * n_sats + i]; - const auto epoch_lo = sat_buffer[static_cast(8) * n_sats + i]; - - // Normalise it into a double-length number. - // NOTE: we use Knuth's EFT here in order to ensure correctness - // even if the fractional part of the epoch is greater in magnitude - // than the main part. - const auto norm_epoch = heyoka::detail::eft_add_knuth(epoch_hi, epoch_lo); + // Do the same for the epoch. + const auto norm_epoch = heyoka::detail::eft_add_knuth(static_cast(tai_epoch_hi), static_cast(tai_epoch_lo)); const auto epoch = dfloat(norm_epoch.first, norm_epoch.second); // Compute the time delta in double-length arithmetic, truncate to a single-length @@ -909,12 +927,10 @@ void sgp4_propagator::operator()(out_2d out, in_1d dates) using tms_vec_size_t = decltype(tms_vec.size()); tms_vec.resize(boost::numeric_cast(dates.extent(0))); - // NOTE: we have a rough estimate of <~50 flops for a single execution of sgp4_date_to_tdelta(). + // NOTE: we have a rough estimate of ~200 flops for a single execution of sgp4_date_to_tdelta(). // Taking the usual figure of ~10'000 clock cycles as minimum threshold to parallelise, we enable - // parallelisation if we have 200 satellites or more. - // NOTE: as an alternative, if we had double-length add/sub in the expression system, - // we could do the transformation on-line as we evaluate the sgp4 propagation. - if (n_sats >= 200u) { + // parallelisation if we have 50 satellites or more. + if (n_sats >= 50u) { oneapi::tbb::parallel_for(oneapi::tbb::blocked_range(0, tms_vec.size()), [&tms_vec, dates, this, n_sats](const auto &range) { for (auto i = range.begin(); i != range.end(); ++i) { @@ -1029,10 +1045,10 @@ void sgp4_propagator::operator()(out_3d out, in_2d dates) using tms_vec_size_t = decltype(tms_vec.size()); tms_vec.resize(boost::safe_numerics::safe(n_evals) * dates.extent(1)); - // NOTE: we have a rough estimate of <~50 flops for a single execution of sgp4_date_to_tdelta(). + // NOTE: we have a rough estimate of ~200 flops for a single execution of sgp4_date_to_tdelta(). // Taking the usual figure of ~10'000 clock cycles as minimum threshold to parallelise, we enable - // parallelisation if tms_vec.size() (i.e., n_evals * n_sats) is 200 or more. - if (tms_vec.size() >= 200u) { + // parallelisation if tms_vec.size() (i.e., n_evals * n_sats) is 50 or more. + if (tms_vec.size() >= 50u) { oneapi::tbb::parallel_for( oneapi::tbb::blocked_range2d(0, dates.extent(0), 0, dates.extent(1)), [&tms_vec, dates, this, n_sats](const auto &range) {