After egg+phl's review.

astronomy/orbit_analysis_test.cpp

@@ -173,11 +173,12 @@ class OrbitAnalysisTest : public ::testing::Test {
   std::tuple<OrbitalElements, OrbitRecurrence, OrbitGroundTrack>
   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();

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> OrbitalElements::ForTrajectory(
   auto const osculating_elements =
       [&primary, &secondary, &trajectory](
           Instant const& time) -> KeplerianElements<PrimaryCentred> {
-    DegreesOfFreedom<PrimaryCentred> const primary_dof{
+    DegreesOfFreedom<PrimaryCentred> const primary_degrees_of_freedom{
         PrimaryCentred::origin, PrimaryCentred::unmoving};
     DegreesOfFreedom<PrimaryCentred> const degrees_of_freedom =
         trajectory.EvaluateDegreesOfFreedom(time);
     return KeplerOrbit<PrimaryCentred>(
-               primary, secondary, degrees_of_freedom - primary_dof, time)
+               primary,
+               secondary,
+               degrees_of_freedom - primary_degrees_of_freedom,
+               time)
         .elements_at_epoch();
   };
 
@@ -137,14 +141,14 @@ absl::StatusOr<OrbitalElements> 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
     // sidereal periods (from aberrant trajectories, see #2811).
     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> 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<Instant,
                                                      Product<Length, Time>,
@@ -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));
         }

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(

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<double, stages> c{{{0, 1}}};
-  static constexpr FixedStrictlyLowerTriangularMatrix<double, stages> a{
-      {{0.5}}};
-  static constexpr FixedVector<double, stages> b̂{{{0.5, 0.5}}};
-  static constexpr FixedVector<double, stages> b{{{1, 0}}};
-};
-
 struct DormandPrince1986RK547FC :  EmbeddedExplicitRungeKutta {
   static constexpr int higher_order = 5;
   static constexpr int lower_order = 4;

ksp_plugin/integrators.hpp

@@ -21,8 +21,8 @@ DiscreteTrajectorySegment<Barycentric>::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<Barycentric>::DownsamplingParameters
 OrbitAnalyserDownsamplingParameters();