From b30a10b4beb015f3098cc29be80391b8327f72a4 Mon Sep 17 00:00:00 2001 From: Pascal Leroy Date: Sun, 12 Feb 2023 18:14:32 +0800 Subject: [PATCH] After egg+phl's review. --- astronomy/orbit_analysis_test.cpp | 11 +++++---- astronomy/orbital_elements_body.hpp | 38 +++++++++++++++++------------ astronomy/orbital_elements_test.cpp | 6 ++--- integrators/methods.hpp | 14 ----------- ksp_plugin/integrators.hpp | 4 +-- 5 files changed, 33 insertions(+), 40 deletions(-) diff --git a/astronomy/orbit_analysis_test.cpp b/astronomy/orbit_analysis_test.cpp index cf80759b00..8fafeb479b 100644 --- a/astronomy/orbit_analysis_test.cpp +++ b/astronomy/orbit_analysis_test.cpp @@ -173,11 +173,12 @@ class OrbitAnalysisTest : public ::testing::Test { std::tuple ElementsAndRecurrence(SP3Orbit const& orbit) { auto const earth_centred_trajectory = EarthCentredTrajectory(orbit); - auto elements = - OrbitalElements::ForTrajectory(*earth_centred_trajectory, - earth_, - MasslessBody{}, - /*sample_dense_elements=*/true).value(); + auto const elements = OrbitalElements::ForTrajectory( + *earth_centred_trajectory, + earth_, + MasslessBody{}, + /*fill_osculating_equinoctial_elements=*/true) + .value(); { auto const identifier = (std::stringstream() << orbit.satellite).str(); diff --git a/astronomy/orbital_elements_body.hpp b/astronomy/orbital_elements_body.hpp index b88e23f3a1..2a027e3105 100644 --- a/astronomy/orbital_elements_body.hpp +++ b/astronomy/orbital_elements_body.hpp @@ -9,7 +9,7 @@ #include "base/jthread.hpp" #include "base/status_utilities.hpp" #include "numerics/quadrature.hpp" -#include "integrators/explicit_runge_kutta_integrator.hpp" +#include "integrators/embedded_explicit_runge_kutta_integrator.hpp" #include "integrators/methods.hpp" #include "physics/discrete_trajectory.hpp" #include "physics/kepler_orbit.hpp" @@ -46,6 +46,7 @@ using quantities::si::Metre; using quantities::si::Milli; using quantities::si::Radian; +constexpr int osculating_equinoctial_elements_per_sidereal_period = 256; constexpr double max_clenshaw_curtis_relative_error = 1.0e-6; constexpr int max_clenshaw_curtis_points = 2000; // Carefully tuned based on MercuryOrbiter test. @@ -71,12 +72,15 @@ absl::StatusOr OrbitalElements::ForTrajectory( auto const osculating_elements = [&primary, &secondary, &trajectory]( Instant const& time) -> KeplerianElements { - DegreesOfFreedom const primary_dof{ + DegreesOfFreedom const primary_degrees_of_freedom{ PrimaryCentred::origin, PrimaryCentred::unmoving}; DegreesOfFreedom const degrees_of_freedom = trajectory.EvaluateDegreesOfFreedom(time); return KeplerOrbit( - primary, secondary, degrees_of_freedom - primary_dof, time) + primary, + secondary, + degrees_of_freedom - primary_degrees_of_freedom, + time) .elements_at_epoch(); }; @@ -137,7 +141,6 @@ absl::StatusOr OrbitalElements::ForTrajectory( osculating_equinoctial_elements, trajectory.t_min(), trajectory.t_max()); RETURN_IF_ERROR(sidereal_period); orbital_elements.sidereal_period_ = sidereal_period.value(); - if (!IsFinite(orbital_elements.sidereal_period_) || orbital_elements.sidereal_period_ <= Time{}) { // Guard against NaN sidereal periods (from hyperbolic orbits) or negative @@ -145,6 +148,7 @@ absl::StatusOr OrbitalElements::ForTrajectory( return absl::OutOfRangeError( "sidereal period is " + DebugString(orbital_elements.sidereal_period_)); } + auto mean_equinoctial_elements = MeanEquinoctialElements(osculating_equinoctial_elements, trajectory.t_min(), @@ -155,8 +159,10 @@ absl::StatusOr OrbitalElements::ForTrajectory( std::move(mean_equinoctial_elements).value(); if (fill_osculating_equinoctial_elements) { - for (Instant t = trajectory.t_min(); t <= trajectory.t_max(); - t += orbital_elements.sidereal_period_ / 256) { + for (Instant t = trajectory.t_min(); + t <= trajectory.t_max(); + t += orbital_elements.sidereal_period_ / + osculating_equinoctial_elements_per_sidereal_period) { orbital_elements.osculating_equinoctial_elements_.push_back( osculating_equinoctial_elements(t)); } @@ -301,7 +307,6 @@ OrbitalElements::MeanEquinoctialElements( // directly, we precompute the integrals from |t_min|. The integral from // |t - period / 2| to |t + period / 2| is then computed by subtraction. - static constexpr int points_per_period = 113; using ODE = ExplicitFirstOrderOrdinaryDifferentialEquation, @@ -386,11 +391,11 @@ OrbitalElements::MeanEquinoctialElements( /*safety_factor=*/0.9)); RETURN_IF_ERROR(instance->Solve(t_max)); - // TODO(egg): Find a nice way to do linear interpolation. + // TODO(egg): Integrate the moving average directly. auto const evaluate_integrals = [&integrals](Instant const& t) -> IntegratedEquinoctialElements { CHECK_LE(t, integrals.back().t); - auto it = std::partition_point( + auto const it = std::partition_point( integrals.begin(), integrals.end(), [&t](IntegratedEquinoctialElements const& elements) { @@ -401,7 +406,7 @@ OrbitalElements::MeanEquinoctialElements( if (it == integrals.begin()) { return high; } else { - IntegratedEquinoctialElements const& low = *--it; + IntegratedEquinoctialElements const& low = *std::prev(it); double const α = (t - low.t) / (high.t - low.t); auto const interpolate = [α, &low, &high](auto const element) { return low.*element + α * (high.*element - low.*element); @@ -499,16 +504,17 @@ inline absl::Status OrbitalElements::ComputePeriodsAndPrecession() { auto const interpolate_function_of_mean_classical_element = [this](auto const f, Instant const& t) { CHECK_LE(t, mean_classical_elements_.back().time); - auto it = std::partition_point(mean_classical_elements_.begin(), - mean_classical_elements_.end(), - [&t](ClassicalElements const& elements) { - return elements.time < t; - }); + auto const it = + std::partition_point(mean_classical_elements_.begin(), + mean_classical_elements_.end(), + [&t](ClassicalElements const& elements) { + return elements.time < t; + }); ClassicalElements const& high = *it; if (it == mean_classical_elements_.begin()) { return f(high); } else { - ClassicalElements const& low = *--it; + ClassicalElements const& low = *std::prev(it); double const α = (t - low.time) / (high.time - low.time); return f(low) + α * (f(high) - f(low)); } diff --git a/astronomy/orbital_elements_test.cpp b/astronomy/orbital_elements_test.cpp index f1542b41b2..7e56ab10a5 100644 --- a/astronomy/orbital_elements_test.cpp +++ b/astronomy/orbital_elements_test.cpp @@ -170,7 +170,7 @@ TEST_F(OrbitalElementsTest, KeplerOrbit) { J2000 + 10 * Day, *ephemeris), spherical_earth, MasslessBody{}, - /*sample_dense_elements=*/true); + /*fill_osculating_equinoctial_elements=*/true); ASSERT_THAT(status_or_elements, IsOk()); OrbitalElements const& elements = status_or_elements.value(); EXPECT_THAT( @@ -268,7 +268,7 @@ TEST_F(OrbitalElementsTest, J2Perturbation) { *ephemeris), oblate_earth, MasslessBody{}, - /*sample_dense_elements=*/true); + /*fill_osculating_equinoctial_elements=*/true); ASSERT_THAT(status_or_elements, IsOk()); OrbitalElements const& elements = status_or_elements.value(); EXPECT_THAT( @@ -372,7 +372,7 @@ TEST_F(OrbitalElementsTest, RealPerturbation) { initial_osculating, J2000, J2000 + mission_duration, *ephemeris), earth, MasslessBody{}, - /*sample_dense_elements=*/true); + /*fill_osculating_equinoctial_elements=*/true); ASSERT_THAT(status_or_elements, IsOk()); OrbitalElements const& elements = status_or_elements.value(); EXPECT_THAT( diff --git a/integrators/methods.hpp b/integrators/methods.hpp index 16ff803ebd..3da4dc15d2 100644 --- a/integrators/methods.hpp +++ b/integrators/methods.hpp @@ -471,20 +471,6 @@ struct DormandالمكاوىPrince1986RKN434FM : {13.0 / 21.0, -20.0 / 27.0, 275.0 / 189.0, -1.0 / 3.0}}}; }; -struct HeunEuler : EmbeddedExplicitRungeKutta { - static constexpr int higher_order = 2; - static constexpr int lower_order = 1; - static constexpr int stages = 2; - static constexpr bool first_same_as_last = false; - static constexpr serialization::AdaptiveStepSizeIntegrator::Kind kind = - serialization::AdaptiveStepSizeIntegrator::DORMAND_PRINCE_1986_RK_547FC; - static constexpr FixedVector c{{{0, 1}}}; - static constexpr FixedStrictlyLowerTriangularMatrix a{ - {{0.5}}}; - static constexpr FixedVector b̂{{{0.5, 0.5}}}; - static constexpr FixedVector b{{{1, 0}}}; -}; - struct DormandPrince1986RK547FC : EmbeddedExplicitRungeKutta { static constexpr int higher_order = 5; static constexpr int lower_order = 4; diff --git a/ksp_plugin/integrators.hpp b/ksp_plugin/integrators.hpp index c0925e6e8b..89781d9fb1 100644 --- a/ksp_plugin/integrators.hpp +++ b/ksp_plugin/integrators.hpp @@ -21,8 +21,8 @@ DiscreteTrajectorySegment::DownsamplingParameters DefaultDownsamplingParameters(); // Parameters for the orbit analyser. Finer-grained than the default to obtain -// reasonable element. The 10 times smaller tolerance results in trajectories -// that are about twice as big. +// reasonable element. The 10'000 times smaller tolerance results in +// trajectories that are about 10 times as big. DiscreteTrajectorySegment::DownsamplingParameters OrbitAnalyserDownsamplingParameters();