Merge pull request #3495 from pleroy/3493

A benchmark for orbital elements

benchmarks/apsides.cpp

@@ -123,9 +123,8 @@ class ApsidesBenchmark : public benchmark::Fixture {
 
 SolarSystem<ICRS>* ApsidesBenchmark::solar_system_2010_ = nullptr;
 Ephemeris<ICRS>* ApsidesBenchmark::ephemeris_ = nullptr;
-OblateBody<ICRS> const* ApsidesBenchmark::ApsidesBenchmark::earth_ = nullptr;
-ContinuousTrajectory<ICRS> const*
-    ApsidesBenchmark::ApsidesBenchmark::earth_trajectory_ = nullptr;
+OblateBody<ICRS> const* ApsidesBenchmark::earth_ = nullptr;
+ContinuousTrajectory<ICRS> const* ApsidesBenchmark::earth_trajectory_ = nullptr;
 DiscreteTrajectory<ICRS>* ApsidesBenchmark::ilrsa_lageos2_trajectory_icrs_ =
     nullptr;
 DiscreteTrajectory<GCRS>* ApsidesBenchmark::ilrsa_lageos2_trajectory_gcrs_ =

benchmarks/benchmarks.vcxproj

@@ -40,6 +40,7 @@
     <ClCompile Include="main.cpp" />
     <ClCompile Include="nearest_neighbour.cpp" />
     <ClCompile Include="newhall.cpp" />
+    <ClCompile Include="orbital_elements.cpp" />
     <ClCompile Include="perspective.cpp" />
     <ClCompile Include="planetarium_plot_methods.cpp" />
     <ClCompile Include="polynomial.cpp" />

benchmarks/benchmarks.vcxproj.filters

@@ -107,6 +107,9 @@
     <ClCompile Include="..\testing_utilities\optimization_test_functions.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="orbital_elements.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="quantities.hpp">

benchmarks/orbital_elements.cpp

@@ -0,0 +1,152 @@
+// .\Release\x64\benchmarks.exe --benchmark_repetitions=3 --benchmark_filter=OrbitalElements --benchmark_min_time=30  // NOLINT(whitespace/line_length)
+
+#include "astronomy/orbital_elements.hpp"
+
+#include <limits>
+
+#include "astronomy/frames.hpp"
+#include "base/not_null.hpp"
+#include "benchmark/benchmark.h"
+#include "geometry/named_quantities.hpp"
+#include "integrators/embedded_explicit_runge_kutta_nyström_integrator.hpp"
+#include "integrators/symmetric_linear_multistep_integrator.hpp"
+#include "integrators/methods.hpp"
+#include "physics/body_centred_non_rotating_dynamic_frame.hpp"
+#include "physics/degrees_of_freedom.hpp"
+#include "physics/discrete_trajectory.hpp"
+#include "physics/ephemeris.hpp"
+#include "physics/kepler_orbit.hpp"
+#include "physics/massive_body.hpp"
+#include "physics/massless_body.hpp"
+#include "physics/oblate_body.hpp"
+#include "physics/solar_system.hpp"
+#include "quantities/si.hpp"
+
+namespace principia {
+namespace astronomy {
+
+using astronomy::GCRS;
+using base::dynamic_cast_not_null;
+using base::make_not_null_unique;
+using base::not_null;
+using geometry::Instant;
+using geometry::Position;
+using integrators::EmbeddedExplicitRungeKuttaNyströmIntegrator;
+using integrators::SymmetricLinearMultistepIntegrator;
+using integrators::methods::DormandالمكاوىPrince1986RKN434FM;
+using integrators::methods::QuinlanTremaine1990Order12;
+using physics::BodyCentredNonRotatingDynamicFrame;
+using physics::DegreesOfFreedom;
+using physics::DiscreteTrajectory;
+using physics::Ephemeris;
+using physics::KeplerianElements;
+using physics::KeplerOrbit;
+using physics::MassiveBody;
+using physics::MasslessBody;
+using physics::OblateBody;
+using physics::SolarSystem;
+using quantities::Time;
+using quantities::si::ArcSecond;
+using quantities::si::Day;
+using quantities::si::Degree;
+using quantities::si::Kilo;
+using quantities::si::Metre;
+using quantities::si::Milli;
+using quantities::si::Minute;
+using quantities::si::Second;
+
+class OrbitalElementsBenchmark : public benchmark::Fixture {
+ protected:
+  static void SetUpFixture() {
+    solar_system_ = std::make_unique<SolarSystem<ICRS>>(
+        SOLUTION_DIR / "astronomy" / "sol_gravity_model.proto.txt",
+        SOLUTION_DIR / "astronomy" /
+            "sol_initial_state_jd_2451545_000000000.proto.txt").release();
+    ephemeris_ = solar_system_->MakeEphemeris(
+        /*accuracy_parameters=*/{/*fitting_tolerance=*/1 * Milli(Metre),
+                                 /*geopotential_tolerance=*/0x1p-24},
+        Ephemeris<ICRS>::FixedStepParameters(
+            SymmetricLinearMultistepIntegrator<QuinlanTremaine1990Order12,
+                                               Position<ICRS>>(),
+            /*step=*/10 * Minute)).release();
+    earth_ = dynamic_cast_not_null<OblateBody<ICRS> const*>(
+        solar_system_->massive_body(*ephemeris_, "Earth"));
+  }
+
+  void SetUp(benchmark::State&) override {
+    static int const set_up_fixture = []() {
+      SetUpFixture();
+      return 0;
+    }();
+  }
+
+  static not_null<std::unique_ptr<DiscreteTrajectory<GCRS>>>
+  EarthCentredTrajectory(
+      KeplerianElements<GCRS> const& initial_osculating_elements,
+      Instant const& initial_time,
+      Instant const& final_time) {
+    BodyCentredNonRotatingDynamicFrame<ICRS, GCRS> gcrs{ephemeris_, earth_};
+    DiscreteTrajectory<ICRS> icrs_trajectory;
+    KeplerOrbit<GCRS> const initial_osculating_orbit{
+        *earth_,
+        MasslessBody{},
+        initial_osculating_elements,
+        initial_time};
+    CHECK_OK(icrs_trajectory.Append(
+        initial_time,
+        gcrs.FromThisFrameAtTime(initial_time)(
+            DegreesOfFreedom<GCRS>{GCRS::origin, GCRS::unmoving} +
+            initial_osculating_orbit.StateVectors(initial_time))));
+    CHECK_OK(ephemeris_->FlowWithAdaptiveStep(
+        &icrs_trajectory,
+        Ephemeris<ICRS>::NoIntrinsicAcceleration,
+        final_time,
+        Ephemeris<ICRS>::AdaptiveStepParameters{
+            EmbeddedExplicitRungeKuttaNyströmIntegrator<
+                DormandالمكاوىPrince1986RKN434FM,
+                Position<ICRS>>(),
+            /*max_steps=*/std::numeric_limits<std::int64_t>::max(),
+            /*length_integration_tolerance=*/1 * Milli(Metre),
+            /*speed_integration_tolerance=*/1 * Milli(Metre) / Second
+        },
+        /*max_ephemeris_steps=*/std::numeric_limits<std::int64_t>::max()));
+    auto result = make_not_null_unique<DiscreteTrajectory<GCRS>>();
+    for (auto const& [time, degrees_of_freedom] : icrs_trajectory) {
+      CHECK_OK(result->Append(
+          time, gcrs.ToThisFrameAtTime(time)(degrees_of_freedom)));
+    }
+    return result;
+  }
+
+  static SolarSystem<ICRS>* solar_system_;
+  static Ephemeris<ICRS>* ephemeris_;
+  static OblateBody<ICRS> const* earth_;
+};
+
+SolarSystem<ICRS>* OrbitalElementsBenchmark::solar_system_ = nullptr;
+Ephemeris<ICRS>* OrbitalElementsBenchmark::ephemeris_ = nullptr;
+OblateBody<ICRS> const* OrbitalElementsBenchmark::earth_ = nullptr;
+
+BENCHMARK_F(OrbitalElementsBenchmark, ComputeOrbitalElements)(
+    benchmark::State& state) {
+  Time const mission_duration = 180 * Day;
+  Instant const final_time = J2000 + mission_duration;
+  CHECK_OK(ephemeris_->Prolong(final_time));
+
+  KeplerianElements<GCRS> initial_osculating;
+  initial_osculating.semimajor_axis = 7000 * Kilo(Metre);
+  initial_osculating.eccentricity = 1e-6;
+  initial_osculating.inclination = 10 * Milli(ArcSecond);
+  initial_osculating.longitude_of_ascending_node = 10 * Degree;
+  initial_osculating.argument_of_periapsis = 20 * Degree;
+  initial_osculating.mean_anomaly = 30 * Degree;
+  for (auto _ : state) {
+    OrbitalElements::ForTrajectory(
+        *EarthCentredTrajectory(initial_osculating, J2000, final_time),
+        *earth_,
+        MasslessBody{}).IgnoreError();
+  }
+}
+
+}  // namespace astronomy
+}  // namespace principia

numerics/piecewise_poisson_series_body.hpp

@@ -689,7 +689,7 @@ InnerProduct(PiecewisePoissonSeries<LValue,
              std::optional<int> max_points) {
   AngularFrequency const max_ω = left.max_ω() + right.max_ω() + weight.max_ω();
   std::optional<int> const max_points_heuristic =
-      MaxPointsHeuristicsForAutomaticClenshawCurtis(
+      quadrature::MaxPointsHeuristicsForAutomaticClenshawCurtis(
           max_ω,
           t_max - t_min,
           clenshaw_curtis_min_points_overall,

numerics/poisson_series.hpp

@@ -387,18 +387,9 @@ InnerProduct(PoissonSeries<LValue,
              Instant const& t_min,
              Instant const& t_max);
 
-// Computes a heuristic for the maximum number of points for an oscillating
-// function.
-std::optional<int> MaxPointsHeuristicsForAutomaticClenshawCurtis(
-    AngularFrequency const& max_ω,
-    Time const& Δt,
-    int min_points_overall,
-    int points_per_period);
-
 }  // namespace internal_poisson_series
 
 using internal_poisson_series::InnerProduct;
-using internal_poisson_series::MaxPointsHeuristicsForAutomaticClenshawCurtis;
 using internal_poisson_series::PoissonSeries;
 
 }  // namespace numerics

numerics/poisson_series_body.hpp

@@ -370,7 +370,7 @@ Norm(PoissonSeries<double,
 
   AngularFrequency const max_ω = 2 * split.slow.max_ω() + weight.max_ω();
   std::optional<int> const max_points =
-      MaxPointsHeuristicsForAutomaticClenshawCurtis(
+      quadrature::MaxPointsHeuristicsForAutomaticClenshawCurtis(
           max_ω,
           t_max - t_min,
           clenshaw_curtis_min_points_overall,
@@ -953,7 +953,7 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
   AngularFrequency const max_ω =
       left_split.slow.max_ω() + right_split.slow.max_ω() + weight.max_ω();
   std::optional<int> const max_points =
-      MaxPointsHeuristicsForAutomaticClenshawCurtis(
+      quadrature::MaxPointsHeuristicsForAutomaticClenshawCurtis(
           max_ω,
           t_max - t_min,
           clenshaw_curtis_min_points_overall,
@@ -980,18 +980,6 @@ typename Hilbert<LValue, RValue>::InnerProductType InnerProduct(
   return (slow_quadrature + fast_quadrature) / (t_max - t_min);
 }
 
-inline std::optional<int> MaxPointsHeuristicsForAutomaticClenshawCurtis(
-    AngularFrequency const& max_ω,
-    Time const& Δt,
-    int min_points_overall,
-    int points_per_period) {
-  return max_ω == AngularFrequency()
-             ? std::optional<int>{}
-             : std::max(min_points_overall,
-                        static_cast<int>(points_per_period * Δt * max_ω /
-                                         (2 * π * Radian)));
-}
-
 }  // namespace internal_poisson_series
 }  // namespace numerics
 }  // namespace principia

numerics/quadrature.hpp

@@ -61,6 +61,7 @@ Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(
 
 using internal_quadrature::AutomaticClenshawCurtis;
 using internal_quadrature::GaussLegendre;
+using internal_quadrature::MaxPointsHeuristicsForAutomaticClenshawCurtis;
 using internal_quadrature::Midpoint;
 
 }  // namespace quadrature

numerics/quadrature_body.hpp

@@ -2,6 +2,7 @@
 
 #include "numerics/quadrature.hpp"
 
+#include <algorithm>
 #include <memory>
 #include <vector>
 
@@ -294,6 +295,18 @@ Primitive<std::invoke_result_t<Function, Argument>, Argument> ClenshawCurtis(
       f, lower_bound, upper_bound, f_cos_N⁻¹π_bit_reversed);
 }
 
+inline std::optional<int> MaxPointsHeuristicsForAutomaticClenshawCurtis(
+    AngularFrequency const& max_ω,
+    Time const& Δt,
+    int min_points_overall,
+    int points_per_period) {
+  return max_ω == AngularFrequency()
+             ? std::optional<int>{}
+             : std::max(min_points_overall,
+                        static_cast<int>(points_per_period * Δt * max_ω /
+                                         (2 * π * Radian)));
+}
+
 template<typename Argument, typename Function>
 Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(
     Function const& f,