mockingbirdnest · Jan 17, 2018 · Jan 17, 2018 · Jan 18, 2018 · Jan 18, 2018 · Jan 18, 2018
diff --git a/benchmarks/newhall.cpp b/benchmarks/newhall.cpp
@@ -34,6 +34,7 @@ void BM_NewhallApproximation(benchmark::State& state) {
   Instant const t_min = t0 + static_cast<double>(random()) * Second;
   Instant const t_max = t_min + static_cast<double>(random()) * Second;
 
+  double error_estimate;
   while (state.KeepRunning()) {
     state.PauseTiming();
     p.clear();
@@ -44,7 +45,10 @@ void BM_NewhallApproximation(benchmark::State& state) {
     }
     state.ResumeTiming();
     auto const series =
-        NewhallApproximationInЧебышёвBasis<double>(degree, p, v, t_min, t_max);
+        NewhallApproximationInЧебышёвBasis<double>(degree,
+                                                   p, v,
+                                                   t_min, t_max,
+                                                   error_estimate);
   }
 }
 

diff --git a/numerics/newhall.hpp b/numerics/newhall.hpp
@@ -15,14 +15,17 @@ using quantities::Variation;
 
 // Computes a Newhall approximation of the given |degree| in the Чебышёв basis.
 // |q| and |v| are the positions and velocities over a constant division of
-// [t_min, t_max].
+// [t_min, t_max].  |error_estimate| gives an estimate of the error between the
+// approximation the input data.  The client probably wants to compute some
+// norm of that estimate.
 template<typename Vector>
 ЧебышёвSeries<Vector> NewhallApproximationInЧебышёвBasis(
     int degree,
     std::vector<Vector> const& q,
     std::vector<Variation<Vector>> const& v,
     Instant const& t_min,
-    Instant const& t_max);
+    Instant const& t_max,
+    Vector& error_estimate);
 
 }  // namespace internal_newhall
 

diff --git a/numerics/newhall_body.hpp b/numerics/newhall_body.hpp
@@ -19,7 +19,8 @@ template<typename Vector>
     std::vector<Vector> const& q,
     std::vector<Variation<Vector>> const& v,
     Instant const& t_min,
-    Instant const& t_max) {
+    Instant const& t_max,
+    Vector& error_estimate) {
   // Only supports 8 divisions for now.
   int const divisions = 8;
   CHECK_EQ(divisions + 1, q.size());
@@ -90,6 +91,7 @@ template<typename Vector>
       break;
   }
   CHECK_EQ(degree + 1, coefficients.size());
+  error_estimate = coefficients[degree];
   return ЧебышёвSeries<Vector>(coefficients, t_min, t_max);
 }
 

diff --git a/numerics/newhall_test.cpp b/numerics/newhall_test.cpp
@@ -18,6 +18,7 @@ namespace principia {
 namespace numerics {
 
 using geometry::Instant;
+using quantities::Abs;
 using quantities::Length;
 using quantities::Speed;
 using quantities::si::Metre;
@@ -38,6 +39,7 @@ class NewhallTest : public ::testing::Test {
   void NewhallApproximationErrors(
       std::function<Length(Instant const)> length_function,
       std::function<Speed(Instant const)> speed_function,
+      std::vector<Length>& length_error_estimates,
       std::vector<Length>& length_absolute_errors,
       std::vector<Speed>& speed_absolute_errors) {
     std::vector<Length> lengths;
@@ -47,12 +49,15 @@ class NewhallTest : public ::testing::Test {
       speeds.push_back(speed_function(t));
     }
 
+    length_error_estimates.clear();
     length_absolute_errors.clear();
     speed_absolute_errors.clear();
     for (int degree = 3; degree <= 17; ++degree) {
+      Length length_error_estimate;
       ЧебышёвSeries<Length> const approximation =
           NewhallApproximationInЧебышёвBasis<Length>(
-              degree, lengths, speeds, t_min_, t_max_);
+              degree, lengths, speeds, t_min_, t_max_, length_error_estimate);
+      length_error_estimates.push_back(Abs(length_error_estimate));
 
       // Compute the absolute error of both functions throughout the interval.
       Length length_absolute_error;
@@ -105,6 +110,7 @@ class NewhallTest : public ::testing::Test {
 };
 
 TEST_F(NewhallTest, ApproximationInЧебышёвBasis) {
+  std::vector<Length> length_error_estimates;
   std::vector<Length> length_absolute_errors;
   std::vector<Speed> speed_absolute_errors;
 
@@ -131,10 +137,27 @@ TEST_F(NewhallTest, ApproximationInЧебышёвBasis) {
 
     NewhallApproximationErrors(length_function,
                                speed_function,
+                               length_error_estimates,
                                length_absolute_errors,
                                speed_absolute_errors);
   }
 
+  ExpectMultipart(length_error_estimates,
+                  ElementsAre(near_length(3.9e1 * Metre),
+                              near_length(1.8e1 * Metre),
+                              near_length(8.3e0 * Metre),
+                              near_length(1.4e0 * Metre),
+                              near_length(7.5e0 * Metre),
+                              near_length(4.2e0 * Metre),
+                              near_length(3.1e-1 * Metre),
+                              near_length(7.5e-1 * Metre),
+                              near_length(1.8e-1 * Metre),
+                              near_length(4.9e-2 * Metre)),
+                  ElementsAre(near_length(2.6e-2 * Metre),
+                              near_length(8.9e-4 * Metre),
+                              near_length(1.6e-3 * Metre),
+                              near_length(2.9e-4 * Metre),
+                              near_length(3.7e-5 * Metre)));
   ExpectMultipart(length_absolute_errors,
                   ElementsAre(near_length(1.7e2 * Metre),
                               near_length(4.7e1 * Metre),
@@ -181,10 +204,27 @@ TEST_F(NewhallTest, ApproximationInЧебышёвBasis) {
 
     NewhallApproximationErrors(length_function,
                                speed_function,
+                               length_error_estimates,
                                length_absolute_errors,
                                speed_absolute_errors);
   }
 
+  ExpectMultipart(length_error_estimates,
+                  ElementsAre(near_length(7.9e-1 * Metre),
+                              near_length(2.5e-1 * Metre),
+                              near_length(5.7e-2 * Metre),
+                              near_length(8.6e-3 * Metre),
+                              near_length(6.2e-4 * Metre),
+                              near_length(4.5e-16 * Metre),
+                              near_length(4.8e-16 * Metre),
+                              near_length(4.8e-16 * Metre),
+                              near_length(1.2e-16 * Metre),
+                              near_length(7.3e-16 * Metre)),
+                  ElementsAre(near_length(1.4e-15 * Metre),
+                              near_length(2.1e-16 * Metre),
+                              near_length(3.2e-16 * Metre),
+                              near_length(1.2e-15 * Metre),
+                              near_length(4.8e-15 * Metre)));
   ExpectMultipart(length_absolute_errors,
                   ElementsAre(near_length(2.0 * Metre),
                               near_length(2.9e-1 * Metre),

diff --git a/physics/continuous_trajectory.hpp b/physics/continuous_trajectory.hpp
@@ -145,7 +145,8 @@ class ContinuousTrajectory : public Trajectory<Frame> {
           std::vector<Displacement<Frame>> const& q,
           std::vector<Velocity<Frame>> const& v,
           Instant const& t_min,
-          Instant const& t_max));
+          Instant const& t_max,
+          Displacement<Frame>& displacement_error_estimate));
 
   // Returns an iterator to the series applicable for the given |time|, or
   // |begin()| if |time| is before the first series or |end()| if |time| is

diff --git a/physics/continuous_trajectory_body.hpp b/physics/continuous_trajectory_body.hpp
@@ -295,7 +295,8 @@ Status ContinuousTrajectory<Frame>::ComputeBestNewhallApproximation(
         std::vector<Displacement<Frame>> const& q,
         std::vector<Velocity<Frame>> const& v,
         Instant const& t_min,
-        Instant const& t_max)) {
+        Instant const& t_max,
+        Displacement<Frame>& error_estimate)) {
   Length const previous_adjusted_tolerance = adjusted_tolerance_;
 
   // If the degree is too old, restart from the lowest degree.  This ensures
@@ -310,12 +311,15 @@ Status ContinuousTrajectory<Frame>::ComputeBestNewhallApproximation(
   }
 
   // Compute the approximation with the current degree.
-  series_.push_back(
-      newhall_approximation(degree_, q, v, last_points_.cbegin()->first, time));
+  Displacement<Frame> displacement_error_estimate;
+  series_.push_back(newhall_approximation(degree_,
+                                          q, v,
+                                          last_points_.cbegin()->first, time,
+                                          displacement_error_estimate));
 
   // Estimate the error.  For initializing |previous_error_estimate|, any value
   // greater than |error_estimate| will do.
-  Length error_estimate = series_.back().last_coefficient().Norm();
+  Length error_estimate = displacement_error_estimate.Norm();
   Length previous_error_estimate = error_estimate + error_estimate;
 
   // If we are in the zone of numerical instabilities and we exceeded the
@@ -343,11 +347,12 @@ Status ContinuousTrajectory<Frame>::ComputeBestNewhallApproximation(
     ++degree_;
     VLOG(1) << "Increasing degree for " << this << " to " <<degree_
             << " because error estimate was " << error_estimate;
-    series_.back() =
-        newhall_approximation(
-            degree_, q, v, last_points_.cbegin()->first, time);
+    series_.back() = newhall_approximation(degree_,
+                                           q, v,
+                                           last_points_.cbegin()->first, time,
+                                           displacement_error_estimate);
     previous_error_estimate = error_estimate;
-    error_estimate = series_.back().last_coefficient().Norm();
+    error_estimate = displacement_error_estimate.Norm();
   }
 
   // If we have entered the zone of numerical instability, go back to the

diff --git a/physics/continuous_trajectory_test.cpp b/physics/continuous_trajectory_test.cpp
@@ -56,8 +56,9 @@ class ContinuousTrajectoryTest : public testing::Test {
       std::vector<Displacement<World>> const& q,
       std::vector<Velocity<World>> const& v,
       Instant const& t_min,
-      Instant const& t_max) {
-    Displacement<World> const error_estimate = error_estimates_->front();
+      Instant const& t_max,
+      Displacement<World>& error_estimate) {
+    error_estimate = error_estimates_->front();
     error_estimates_->pop_front();
     return ЧебышёвSeries<Displacement<World>>({error_estimate}, t_min, t_max);
   }