Merge pull request #3785 from pleroy/3569

Don't allow the adaptive step integrator to use steps that are too long

astronomy/orbit_analysis_test.cpp

@@ -370,9 +370,9 @@ TEST_F(OrbitAnalysisTest, 北斗MEO) {
   EXPECT_THAT(elements.mean_inclination_interval().midpoint(),
               IsNear(55.10_(1) * Degree));
   EXPECT_THAT(elements.mean_eccentricity_interval().midpoint(),
-              IsNear(0.000557_(1)));
+              IsNear(0.000554_(1)));
   EXPECT_THAT(elements.mean_argument_of_periapsis_interval().midpoint(),
-              IsNear(0.4737_(1) * Degree));
+              IsNear(0.7849_(1) * Degree));
 }
 
 // COSPAR ID 2016-030A.
@@ -404,20 +404,20 @@ TEST_F(OrbitAnalysisTest, GalileoNominalSlot) {
             RelativeErrorFrom(nominal_nodal_precession, IsNear(0.011_(1)))));
   EXPECT_THAT(2 * π * Radian / elements.anomalistic_period(),
               AllOf(AbsoluteErrorFrom(nominal_anomalistic_mean_motion,
-                                      IsNear(0.63_(1) * Degree / Day)),
+                                      IsNear(0.66_(1) * Degree / Day)),
                     RelativeErrorFrom(nominal_anomalistic_mean_motion,
-                                      IsNear(0.00102_(1)))));
+                                      IsNear(0.00108_(1)))));
 
   EXPECT_THAT(elements.mean_semimajor_axis_interval().midpoint(),
               AbsoluteErrorFrom(29'599.8 * Kilo(Metre),
                                 IsNear(0.33_(1) * Kilo(Metre))));
   EXPECT_THAT(elements.mean_semimajor_axis_interval().measure(),
-              IsNear(00'000.084_(1) * Kilo(Metre)));
+              IsNear(00'000.089_(1) * Kilo(Metre)));
 
   EXPECT_THAT(elements.mean_eccentricity_interval().midpoint(),
               IsNear(0.000'17_(1)));  // Nominal: 0.0.
   EXPECT_THAT(elements.mean_eccentricity_interval().measure(),
-              IsNear(0.000'020_(1)));
+              IsNear(0.000'018_(1)));
 
   EXPECT_THAT(elements.mean_inclination_interval().midpoint(),
               AbsoluteErrorFrom(56.0 * Degree, IsNear(0.61_(1) * Degree)));
@@ -436,7 +436,7 @@ TEST_F(OrbitAnalysisTest, GalileoNominalSlot) {
   EXPECT_THAT(elements.mean_argument_of_periapsis_interval().midpoint(),
               IsNear(89_(1) * Degree));
   EXPECT_THAT(elements.mean_argument_of_periapsis_interval().measure(),
-              IsNear(7.7_(1) * Degree));
+              IsNear(7.8_(1) * Degree));
 
   // Since the reference parameters conventionally set ω = 0, the given mean
   // anomaly is actually the mean argument of latitude; in order to get numbers
@@ -485,9 +485,9 @@ TEST_F(OrbitAnalysisTest, GalileoExtendedSlot) {
 
   EXPECT_THAT(elements.mean_semimajor_axis_interval().midpoint(),
               AbsoluteErrorFrom(27'977.6 * Kilo(Metre),
-                                IsNear(0.0519_(1) * Kilo(Metre))));
+                                IsNear(0.0485_(1) * Kilo(Metre))));
   EXPECT_THAT(elements.mean_semimajor_axis_interval().measure(),
-              IsNear(00'000.099_(1) * Kilo(Metre)));
+              IsNear(00'000.101_(1) * Kilo(Metre)));
 
   EXPECT_THAT(elements.mean_eccentricity_interval().midpoint(),
               AbsoluteErrorFrom(0.162, IsNear(0.0041_(1))));
@@ -581,7 +581,7 @@ TEST_F(OrbitAnalysisTest, TOPEXPoséidon) {
   // bit less than 3 m above the nominal value around that time.
   EXPECT_THAT(
       elements.mean_semimajor_axis_interval().midpoint(),
-      DifferenceFrom(7714.42938 * Kilo(Metre), IsNear(2.41_(1) * Metre)));
+      DifferenceFrom(7714.42938 * Kilo(Metre), IsNear(2.63_(1) * Metre)));
   // Reference inclination from the legend of figure 9 of [BSFL98]; that
   // value is given as 66.040° in table 1 of [BSFL98], 66.039° in [BS96], and
   // 66.04° in [Ben97].
@@ -607,10 +607,10 @@ TEST_F(OrbitAnalysisTest, TOPEXPoséidon) {
   // theoretical and observed mean e and ω vary between 40 ppm and 140 ppm, and
   // between 60° and 120°, respectively.
   EXPECT_THAT(elements.mean_eccentricity_interval(),
-              AllOf(Field(&Interval<double>::min, IsNear(85e-6_(1))),
+              AllOf(Field(&Interval<double>::min, IsNear(83e-6_(1))),
                     Field(&Interval<double>::max, IsNear(109e-6_(1)))));
   EXPECT_THAT(elements.mean_argument_of_periapsis_interval(),
-              AllOf(Field(&Interval<Angle>::min, IsNear(74.4_(1) * Degree)),
+              AllOf(Field(&Interval<Angle>::min, IsNear(73.8_(1) * Degree)),
                     Field(&Interval<Angle>::max, IsNear(99.2_(1) * Degree))));
 
   // Nominal longitude of the equatorial crossing of the first ascending pass
@@ -693,7 +693,7 @@ TEST_F(OrbitAnalysisTest, SPOT5) {
   EXPECT_THAT(elements.mean_eccentricity_interval().midpoint(),
               IsNear(0.0012_(1)));
   EXPECT_THAT(elements.mean_argument_of_periapsis_interval().midpoint(),
-              IsNear(89.69_(1) * Degree));
+              IsNear(89.63_(1) * Degree));
 
   // The nominal mean solar times of the nodes are 22:30 ascending, 10:30
   // descending.
@@ -728,7 +728,7 @@ TEST_F(OrbitAnalysisTest, Sentinel3A) {
   EXPECT_THAT(elements.mean_eccentricity_interval().midpoint(),
               IsNear(0.0011_(1)));
   EXPECT_THAT(elements.mean_argument_of_periapsis_interval().midpoint(),
-              IsNear(90.03_(1) * Degree));
+              IsNear(90.00_(1) * Degree));
 
   // The nominal mean solar times of the nodes are 22:00 ascending, 10:00
   // descending.

astronomy/orbital_elements_body.hpp

@@ -357,11 +357,25 @@ OrbitalElements::MeanEquinoctialElements(
   };
 
   auto const tolerance_to_error_ratio =
-      [](Time const& step,
+      [period](Time const& step,
          ODE::State const& state,
          ODE::State::Error const& error) -> double {
+    // When the trajectory is very regular, the integrator is "too good" at
+    // approximating it, which causes the output of the integration to be very
+    // sparse, to the point where it confuses unwinding (because we have more
+    // than half a revolution between points).  To avoid this we reduce the
+    // tolerance-to-error ratio exponentially above 1/3 of the period.  For a
+    // step of |period / 2|, the reduction is e^-3.
+    double braking_factor = 1.0;
+    if (3 * step >= period) {
+      braking_factor = std::exp(6 - 18 * step / period);
+    }
+
+    // The braking factor can be very small (even 0) for large steps.  In that
+    // case we want to reject the step, but not drive it all the way to 0,
+    // hence the |std::max|.
     auto const& [Δa, Δh, Δk, Δλ, Δp, Δq, Δpʹ, Δqʹ] = error;
-    return eerk_a_tolerance / Abs(Δa);
+    return std::max(0.5, braking_factor * eerk_a_tolerance / Abs(Δa));
   };
 
   auto const initial_integration =
@@ -550,6 +564,10 @@ inline absl::Status OrbitalElements::ComputePeriodsAndPrecession() {
   anomalistic_period_ = 2 * π * Radian * Δt³ / (12 * ʃ_Mt_dt);
   nodal_period_ = 2 * π * Radian * Δt³ / (12 * ʃ_ut_dt);
   nodal_precession_ = 12 * ʃ_Ωt_dt / Δt³;
+
+  CHECK_LT(0 * Second, anomalistic_period_);
+  CHECK_LT(0 * Second, nodal_period_);
+
   return absl::OkStatus();
 }
 

astronomy/лидов_古在_test.cpp

@@ -181,7 +181,7 @@ TEST_F(Лидов古在Test, MercuryOrbiter) {
   // pumping energy into nor out of it.  The true values are 14'910.01 and
   // 14'910.28 km.
   EXPECT_THAT(elements.mean_semimajor_axis_interval().min,
-              IsNear(14'910.02_(1) * Kilo(Metre)));
+              IsNear(14'910.01_(1) * Kilo(Metre)));
   EXPECT_THAT(elements.mean_semimajor_axis_interval().max,
               IsNear(14'910.28_(1) * Kilo(Metre)));
 

journal/player_test.cpp

@@ -113,7 +113,7 @@ TEST_F(PlayerTest, DISABLED_SECULAR_Scan) {
 // |method_out_return| protocol buffers.
 TEST_F(PlayerTest, DISABLED_SECULAR_Debug) {
   std::string path =
-      R"(P:\Public Mockingbird\Principia\Issues\3782\JOURNAL.20231025-212010)";  // NOLINT
+      R"(P:\Public Mockingbird\Principia\Issues\3569\JOURNAL.20231026-191142)";  // NOLINT
   Player player(path);
   int count = 0;
   while (player.Play(count)) {

testing_utilities/matchers.hpp

@@ -5,6 +5,7 @@
 #include <string>
 
 #include "absl/status/status.h"
+#include "absl/status/statusor.h"
 #include "gmock/gmock.h"
 #include "google/protobuf/message.h"
 #include "google/protobuf/util/message_differencer.h"
@@ -21,6 +22,15 @@ namespace internal {
 #define EXPECT_OK(value) \
   EXPECT_THAT((value), ::principia::testing_utilities::_matchers::IsOk());
 
+inline absl::Status StatusOf(absl::Status const& s) {
+  return s;
+}
+
+template<typename T>
+absl::Status StatusOf(absl::StatusOr<T> const& s) {
+  return s.status();
+}
+
 MATCHER_P(EqualsProto,
           expected,
           std::string(negation ? "is not" : "is") + " equal to:\n" +
@@ -37,7 +47,11 @@ MATCHER_P(EqualsProto,
 
 MATCHER(IsOk,
         std::string(negation ? "is not" : "is") + " ok") {
-  return arg.ok();
+  if (!arg.ok()) {
+    *result_listener << "Status is " << StatusOf(arg);
+    return false;
+  }
+  return true;
 }
 
 MATCHER_P(IsOkAndHolds,