better accuracy parameters

physics/ephemeris.hpp

@@ -75,7 +75,7 @@ class Ephemeris {
     // Implicit for compatibility.
     AccuracyParameters(Length const& fitting_tolerance);  // NOLINT
     AccuracyParameters(Length const& fitting_tolerance,
-                       serialization::Numerics::Mode geopotential_mode);
+                       double geopotential_tolerance);
 
     void WriteToMessage(
         not_null<serialization::Ephemeris::AccuracyParameters*> const
@@ -85,8 +85,7 @@ class Ephemeris {
 
    private:
     Length fitting_tolerance_;
-    serialization::Numerics::Mode geopotential_mode_ =
-        serialization::Numerics::PRECISE;
+    double geopotential_tolerance_ = 0;
     friend class Ephemeris<Frame>;
   };
 

physics/ephemeris_body.hpp

@@ -121,16 +121,16 @@ Ephemeris<Frame>::AccuracyParameters::AccuracyParameters(
 template<typename Frame>
 Ephemeris<Frame>::AccuracyParameters::AccuracyParameters(
     Length const& fitting_tolerance,
-    serialization::Numerics::Mode const geopotential_mode)
+    double const geopotential_tolerance)
     : fitting_tolerance_(fitting_tolerance),
-      geopotential_mode_(geopotential_mode) {}
+      geopotential_tolerance_(geopotential_tolerance) {}
 
 template<typename Frame>
 void Ephemeris<Frame>::AccuracyParameters::WriteToMessage(
     not_null<serialization::Ephemeris::AccuracyParameters*> const message)
     const {
   fitting_tolerance_.WriteToMessage(message->mutable_fitting_tolerance());
-  message->set_geopotential_mode(geopotential_mode_);
+  message->set_geopotential_tolerance(geopotential_tolerance_);
 }
 
 template<typename Frame>
@@ -139,7 +139,7 @@ Ephemeris<Frame>::AccuracyParameters::ReadFromMessage(
     serialization::Ephemeris::AccuracyParameters const& message) {
   return AccuracyParameters(
       Length::ReadFromMessage(message.fitting_tolerance()),
-      message.geopotential_mode());
+      message.geopotential_tolerance());
 }
 
 template<typename Frame>
@@ -302,7 +302,8 @@ Ephemeris<Frame>::Ephemeris(
     if (body->is_oblate()) {
       geopotentials_.emplace(
           geopotentials_.cbegin(),
-          dynamic_cast_not_null<OblateBody<Frame> const*>(body.get()));
+          dynamic_cast_not_null<OblateBody<Frame> const*>(body.get()),
+          accuracy_parameters_.geopotential_tolerance_);
       // Inserting at the beginning of the vectors is O(N).
       bodies_.insert(bodies_.begin(), std::move(body));
       trajectories_.insert(trajectories_.begin(), trajectory);

physics/geopotential.hpp

@@ -25,6 +25,37 @@ using quantities::Length;
 using quantities::Quotient;
 using quantities::Square;
 
+// Specification of the damping of a spherical harmonic, acting as a radial
+// multiplier on the potential:
+//   V_damped = σ(|r|) V(r).
+struct HarmonicDamping final {
+  explicit HarmonicDamping() = default;
+  explicit HarmonicDamping(Length const& inner_threshold);
+
+  // Above this threshold, the contribution to the potential from this
+  // harmonic is 0, i.e., σ = 0.
+  Length outer_threshold = Infinity<Length>();
+  // Below this threshold, the contribution to the potential from this
+  // harmonic is undamped, σ = 1.
+  // inner_threshold = outer_threshold / 3.
+  Length inner_threshold = Infinity<Length>();
+  // For r in [outer_threshold, inner_threshold], σ is a polynomial with the
+  // following coefficients in monomial basis.
+  // The constant term is always 0, and is thus ignored in the evaluation.
+  Derivatives<double, Length, 4> sigmoid_coefficients;
+
+  // Sets σℜ_over_r and grad_σℜ according to σ as defined by |*this|.
+  template<typename Frame>
+  void ComputeDampedRadialQuantities(
+      Length const& r_norm,
+      Square<Length> const& r²,
+      Vector<double, Frame> const& r_normalized,
+      Inverse<Square<Length>> const& ℜ_over_r,
+      Inverse<Square<Length>> const& ℜʹ,
+      Inverse<Square<Length>>& σℜ_over_r,
+      Vector<Inverse<Square<Length>>, Frame>& grad_σℜ) const;
+};
+
 // Representation of the geopotential model of an oblate body.
 template<typename Frame>
 class Geopotential {
@@ -77,36 +108,6 @@ class Geopotential {
   template<typename>
   struct AllDegrees;
 
-  // Specification of the damping of a spherical harmonic, acting as
-  // a radial multiplier on the potential:
-  //   V_damped = σ(|r|) V(r).
-  struct HarmonicDamping {
-    explicit HarmonicDamping() = default;
-    explicit HarmonicDamping(Length const& inner_threshold);
-
-    // Above this threshold, the contribution to the potential from this
-    // harmonic is 0, i.e., σ = 0.
-    Length const outer_threshold = Infinity<Length>();
-    // Below this threshold, the contribution to the potential from this
-    // harmonic is undamped, σ = 1.
-    // inner_threshold = outer_threshold / 3.
-    Length const inner_threshold = Infinity<Length>();
-    // For r in [outer_threshold, inner_threshold], σ is a polynomial with the
-    // following coefficients in monomial basis.
-    // The constant term is always 0, and is thus ignored in the evaluation.
-    Derivatives<double, Length, 4> const sigmoid_coefficients;
-
-    // Sets σℜ_over_r and grad_σℜ according to σ as defined by |*this|.
-    void ComputeDampedRadialQuantities(
-        Length const& r_norm,
-        Square<Length> const& r²,
-        Vector<double, Frame> const& r_normalized,
-        Inverse<Square<Length>> const& ℜ_over_r,
-        Inverse<Square<Length>> const& ℜʹ,
-        Inverse<Square<Length>>& σℜ_over_r,
-        Vector<Inverse<Square<Length>>, Frame>& grad_σℜ) const;
-  };
-
   // If z is a unit vector along the axis of rotation, and r a vector from the
   // center of |body_| to some point in space, the acceleration computed here
   // is:

physics/geopotential_body.hpp

@@ -38,6 +38,50 @@ using quantities::Sin;
 using quantities::SIUnit;
 
 // The notation in this file follows documentation/Geopotential.pdf.
+
+inline HarmonicDamping::HarmonicDamping(
+    Length const& inner_threshold)
+    : outer_threshold(inner_threshold * 3),
+      inner_threshold(outer_threshold),
+      sigmoid_coefficients{0,
+                           9 / (4 * inner_threshold),
+                           -3 / (2 * Pow<2>(inner_threshold)),
+                           1 / (4 * Pow<3>(inner_threshold))} {}
+
+template<typename Frame>
+void HarmonicDamping::ComputeDampedRadialQuantities(
+      Length const& r_norm,
+      Square<Length> const& r²,
+      Vector<double, Frame> const& r_normalized,
+      Inverse<Square<Length>> const& ℜ_over_r,
+      Inverse<Square<Length>> const& ℜʹ,
+      Inverse<Square<Length>>& σℜ_over_r,
+      Vector<Inverse<Square<Length>>, Frame>& grad_σℜ) const {
+  Length const& s1 = outer_threshold;
+  Length const& s0 = inner_threshold;
+  if (r_norm <= s0) {
+    // Below the inner threshold, σ = 1.
+    σℜ_over_r = ℜ_over_r;
+    grad_σℜ = ℜʹ * r_normalized;
+  } else {
+    auto const& c = sigmoid_coefficients;
+    Derivative<double, Length> const c1 = std::get<1>(c);
+    Derivative<double, Length, 2> const c2 = std::get<2>(c);
+    Derivative<double, Length, 3> const c3 = std::get<3>(c);
+    auto const r³ = r² * r_norm;
+    double const c3r³ = c3 * r³;
+    double const c2r² = c2 * r²;
+    double const c1r = c1 * r_norm;
+    double const σ = c3r³ + c2r² + c1r;
+    double const σʹr = 3 * c3r³ + 2 * c2r² + c1r;
+
+    σℜ_over_r = σ * ℜ_over_r;
+    // Writing this as σ′ℜ + ℜ′σ rather than ℜ∇σ + σ∇ℜ turns some vector
+    // operations into scalar ones.
+    grad_σℜ = (σʹr * ℜ_over_r + ℜʹ * σ) * r_normalized;
+  }
+}
+
 template<typename Frame>
 template<int size>
 struct Geopotential<Frame>::Precomputations {
@@ -243,50 +287,6 @@ auto Geopotential<Frame>::DegreeNOrderM<size, degree, order>::Acceleration(
   }
 }
 
-template<typename Frame>
-Geopotential<Frame>::HarmonicDamping::HarmonicDamping(
-    Length const& inner_threshold)
-    : outer_threshold(inner_threshold * 3),
-      inner_threshold(outer_threshold),
-      sigmoid_coefficients{0,
-                           9 / (4 * inner_threshold),
-                           -3 / (2 * Pow<2>(inner_threshold)),
-                           1 / (4 * Pow<3>(inner_threshold))} {}
-
-template<typename Frame>
-void Geopotential<Frame>::HarmonicDamping::ComputeDampedRadialQuantities(
-      Length const& r_norm,
-      Square<Length> const& r²,
-      Vector<double, Frame> const& r_normalized,
-      Inverse<Square<Length>> const& ℜ_over_r,
-      Inverse<Square<Length>> const& ℜʹ,
-      Inverse<Square<Length>>& σℜ_over_r,
-      Vector<Inverse<Square<Length>>, Frame>& grad_σℜ) const {
-  Length const& s1 = outer_threshold;
-  Length const& s0 = inner_threshold;
-  if (r_norm <= s0) {
-    // Below the inner threshold, σ = 1.
-    σℜ_over_r = ℜ_over_r;
-    grad_σℜ = ℜʹ * r_normalized;
-  } else {
-    auto const& c = sigmoid_coefficients;
-    Derivative<double, Length> const c1 = std::get<1>(c);
-    Derivative<double, Length, 2> const c2 = std::get<2>(c);
-    Derivative<double, Length, 3> const c3 = std::get<3>(c);
-    auto const r³ = r² * r_norm;
-    double const c3r³ = c3 * r³;
-    double const c2r² = c2 * r²;
-    double const c1r = c1 * r_norm;
-    double const σ = c3r³ + c2r² + c1r;
-    double const σʹr = 3 * c3r³ + 2 * c2r² + c1r;
-
-    σℜ_over_r = σ * ℜ_over_r;
-    // Writing this as σ′ℜ + ℜ′σ rather than ℜ∇σ + σ∇ℜ turns some vector
-    // operations into scalar ones.
-    grad_σℜ = (σʹr * ℜ_over_r + ℜʹ * σ) * r_normalized;
-  }
-}
-
 template<typename Frame>
 template<int size, int degree, int... orders>
 auto Geopotential<Frame>::
@@ -492,20 +492,20 @@ Geopotential<Frame>::Geopotential(not_null<OblateBody<Frame> const*> body,
   tesseral_damping_ = HarmonicDamping(Length{});
   Length tesseral_threshold;
   bool is_tesseral = false;
-  degree_damping_.push_back();
-  degree_damping_.push_back();
+  degree_damping_.emplace_back();
+  degree_damping_.emplace_back();
   for (int n = 2; n <= body_->geopotential_degree(); ++n) {
     Length degree_n_threshold = degree_damping_[n - 1].inner_threshold;
-    for (m = 0; m <= n; ++m) {
+    for (int m = 0; m <= n; ++m) {
       double const max_abs_Pnm =
           MaxAbsNormalizedAssociatedLegendreFunction[n][m];
       double const Cnm = body->cos()[n][m];
       double const Snm = body->sin()[n][m];
       // TODO(egg): write a rootn.
       Length const r =
           body->reference_radius() *
-          std::exp(
-              (max_abs_Pnm * (1 + n) * Sqrt(Pow<2>(Cnm) + Pow<2>(Snm))) / ε,
+          std::pow(
+              (max_abs_Pnm * (n + 1) * Sqrt(Pow<2>(Cnm) + Pow<2>(Snm))) / ε,
               1.0 / n);
       if (m == 0 || is_tesseral) {
         degree_n_threshold = std::max(r, degree_n_threshold);

physics/geopotential_test.cpp

@@ -97,7 +97,7 @@ TEST_F(GeopotentialTest, J2) {
       OblateBody<World>(massive_body_parameters_,
                         rotating_body_parameters_,
                         OblateBody<World>::Parameters(/*j2=*/6, 1 * Metre));
-  Geopotential<World> const geopotential(&body);
+  Geopotential<World> const geopotential(&body, /*tolerance=*/0);
 
   // The acceleration at a point located on the axis is along the axis.
   {
@@ -169,7 +169,7 @@ TEST_F(GeopotentialTest, C22S22) {
                         rotating_body_parameters_,
                         OblateBody<World>::Parameters::ReadFromMessage(
                             message, 1 * Metre));
-  Geopotential<World> const geopotential(&body);
+  Geopotential<World> const geopotential(&body, /*tolerance=*/0);
 
   // The acceleration at a point located on the axis is along the axis for the
   // (2, 2) harmonics.
@@ -224,7 +224,7 @@ TEST_F(GeopotentialTest, J3) {
                         rotating_body_parameters_,
                         OblateBody<World>::Parameters::ReadFromMessage(
                             message, 1 * Metre));
-  Geopotential<World> const geopotential(&body);
+  Geopotential<World> const geopotential(&body, /*tolerance=*/0);
 
   // The acceleration at a point located on the axis is along the axis.
   {
@@ -279,7 +279,7 @@ TEST_F(GeopotentialTest, TestVector) {
       rotating_body_parameters,
       OblateBody<ICRS>::Parameters::ReadFromMessage(
           earth_message.geopotential(), earth_reference_radius));
-  Geopotential<ICRS> const geopotential(&earth);
+  Geopotential<ICRS> const geopotential(&earth, /*tolerance=*/0);
 
   // This test vector is from Kuga & Carrara, "Fortran- and C-codes for higher
   // order and degree geopotential computation",

serialization/physics.proto

@@ -98,7 +98,7 @@ message Ephemeris {
   // Added in Ἐρατοσθένης.
   message AccuracyParameters {
     required Quantity fitting_tolerance = 1;
-    required Numerics.Mode geopotential_mode = 2;
+    required double geopotential_tolerance = 2;
   }
   message AdaptiveStepParameters {
     required AdaptiveStepSizeIntegrator integrator = 1;

suppress_useless_warnings.props

@@ -5,7 +5,7 @@
   <PropertyGroup />
   <ItemDefinitionGroup>
     <ClCompile>
-      <DisableSpecificWarnings>4018;4244;4251;4267;4373;4503;4506;4554;4566;4624;4800;4996;%(DisableSpecificWarnings)</DisableSpecificWarnings>
+      <DisableSpecificWarnings>4018;4244;4251;4267;4373;4503;4506;4554;4566;4624;4723;4800;4996;%(DisableSpecificWarnings)</DisableSpecificWarnings>
     </ClCompile>
   </ItemDefinitionGroup>
   <ItemGroup />