mockingbirdnest · Apr 1, 2018 · Apr 1, 2018 · Apr 1, 2018 · Apr 1, 2018 · Apr 1, 2018
Showing with 251 additions and 0 deletions.

+2 −0 benchmarks/benchmarks.vcxproj

+6 −0 benchmarks/benchmarks.vcxproj.filters

+243 −0 benchmarks/planetarium_plot_methods.cpp
diff --git a/benchmarks/benchmarks.vcxproj b/benchmarks/benchmarks.vcxproj
@@ -9,13 +9,15 @@
   <Import Project="$(SolutionDir)principia.props" />
   <ItemGroup>
     <ClCompile Include="..\base\status.cpp" />
+    <ClCompile Include="..\ksp_plugin\planetarium.cpp" />
     <ClCompile Include="dynamic_frame.cpp" />
     <ClCompile Include="embedded_explicit_runge_kutta_nyström_integrator.cpp" />
     <ClCompile Include="ephemeris.cpp" />
     <ClCompile Include="hexadecimal.cpp" />
     <ClCompile Include="main.cpp" />
     <ClCompile Include="newhall.cpp" />
     <ClCompile Include="perspective.cpp" />
+    <ClCompile Include="planetarium_plot_methods.cpp" />
     <ClCompile Include="polynomial.cpp" />
     <ClCompile Include="quantities.cpp" />
     <ClCompile Include="symplectic_runge_kutta_nyström_integrator.cpp" />

diff --git a/benchmarks/benchmarks.vcxproj.filters b/benchmarks/benchmarks.vcxproj.filters
@@ -47,6 +47,12 @@
     <ClCompile Include="newhall.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="..\ksp_plugin\planetarium.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="planetarium_plot_methods.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="quantities.hpp">

diff --git a/benchmarks/planetarium_plot_methods.cpp b/benchmarks/planetarium_plot_methods.cpp
@@ -0,0 +1,243 @@
+
+// .\Release\x64\benchmarks.exe --benchmark_repetitions=3 --benchmark_filter=Planetarium  // NOLINT(whitespace/line_length)
+
+#include "ksp_plugin/planetarium.hpp"
+
+#include <algorithm>
+
+#include "astronomy/time_scales.hpp"
+#include "benchmark/benchmark.h"
+#include "physics/body_centred_non_rotating_dynamic_frame.hpp"
+#include "physics/solar_system.hpp"
+#include "testing_utilities/solar_system_factory.hpp"
+
+namespace principia {
+namespace geometry {
+
+using astronomy::operator""_UT1;
+using astronomy::operator""_TT;
+using base::make_not_null_unique;
+using base::not_null;
+using geometry::Bivector;
+using geometry::Perspective;
+using geometry::RigidTransformation;
+using geometry::RP2Lines;
+using geometry::Vector;
+using geometry::Velocity;
+using integrators::SymmetricLinearMultistepIntegrator;
+using integrators::methods::Quinlan1999Order8A;
+using integrators::methods::QuinlanTremaine1990Order12;
+using ksp_plugin::Camera;
+using ksp_plugin::Barycentric;
+using ksp_plugin::Navigation;
+using ksp_plugin::NavigationFrame;
+using ksp_plugin::Planetarium;
+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::SolarSystem;
+using quantities::Angle;
+using quantities::Cos;
+using quantities::Infinity;
+using quantities::Length;
+using quantities::Sin;
+using quantities::Time;
+using quantities::si::ArcMinute;
+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::Radian;
+using quantities::si::Second;
+using testing_utilities::SolarSystemFactory;
+
+namespace {
+
+constexpr Length near = 40'000 * Kilo(Metre);
+constexpr Length far = 400'000 * Kilo(Metre);
+constexpr Length focal = 1 * Metre;
+
+Perspective<Navigation, Camera> PolarPerspective(
+    Length const distance_from_earth) {
+  return {
+      RigidTransformation<Navigation, Camera>(
+          Navigation::origin + Displacement<Navigation>(
+                                   {0 * Metre, 0 * Metre, distance_from_earth}),
+          Camera::origin,
+          Rotation<Navigation, Camera>(Vector<double, Navigation>({1, 0, 0}),
+                                       Vector<double, Navigation>({0, -1, 0}),
+                                       Bivector<double, Navigation>({0, 0, -1}))
+              .Forget()),
+      focal};
+}
+
+Perspective<Navigation, Camera> EquatorialPerspective(
+    Length const distance_from_earth) {
+  return {
+      RigidTransformation<Navigation, Camera>(
+          Navigation::origin + Displacement<Navigation>(
+                                   {0 * Metre, distance_from_earth, 0 * Metre}),
+          Camera::origin,
+          Rotation<Navigation, Camera>(Vector<double, Navigation>({1, 0, 0}),
+                                       Vector<double, Navigation>({0, 0, 1}),
+                                       Bivector<double, Navigation>({0, -1, 0}))
+              .Forget()),
+      focal};
+}
+
+class Satellites {
+ public:
+  Satellites()
+      : solar_system_(make_not_null_unique<SolarSystem<Barycentric>>(
+            SOLUTION_DIR / "astronomy" / "sol_gravity_model.proto.txt",
+            SOLUTION_DIR / "astronomy" /
+                "sol_initial_state_jd_2451545_000000000.proto.txt",
+            /*ignore_frame=*/true)),
+        ephemeris_(
+            solar_system_->MakeEphemeris(/*fitting_tolerance=*/1 * Milli(Metre),
+                                         EphemerisParameters())),
+        earth_(solar_system_->massive_body(
+            *ephemeris_,
+            SolarSystemFactory::name(SolarSystemFactory::Earth))),
+        earth_centred_inertial_(
+            make_not_null_unique<
+                BodyCentredNonRotatingDynamicFrame<Barycentric, Navigation>>(
+                ephemeris_.get(),
+                earth_)) {
+    // Two-line elements for GOES-8:
+    // 1 23051U 94022A   00004.06628221 -.00000243  00000-0  00000-0 0  9630
+    // 2 23051   0.4232  97.7420 0004776 192.8349 121.5613  1.00264613 28364
+    constexpr Instant goes_8_epoch = "JD2451548.56628221"_UT1;
+    KeplerianElements<Barycentric> goes_8_elements;
+    goes_8_elements.inclination = 0.4232 * Degree;
+    goes_8_elements.longitude_of_ascending_node = 97.7420 * Degree;
+    goes_8_elements.eccentricity = 0.0004776;
+    goes_8_elements.argument_of_periapsis = 192.8349 * Degree;
+    goes_8_elements.mean_anomaly = 121.5613 * Degree;
+    goes_8_elements.mean_motion = 1.00264613 * (2 * π * Radian / Day);
+
+    ephemeris_->Prolong(goes_8_epoch);
+    KeplerOrbit<Barycentric> const goes_8_orbit(
+        *earth_, MasslessBody{}, goes_8_elements, goes_8_epoch);
+    goes_8_trajectory_.Append(
+        goes_8_epoch,
+        ephemeris_->trajectory(earth_)->EvaluateDegreesOfFreedom(goes_8_epoch) +
+            goes_8_orbit.StateVectors(goes_8_epoch));
+    auto goes_8_instance = ephemeris_->NewInstance(
+        {&goes_8_trajectory_},
+        Ephemeris<Barycentric>::NoIntrinsicAccelerations,
+        HistoryParameters());
+    CHECK_OK(ephemeris_->FlowWithFixedStep(goes_8_epoch + 100 * Day,
+                                           *goes_8_instance));
+  }
+
+  DiscreteTrajectory<Barycentric> const& goes_8_trajectory() const {
+    return goes_8_trajectory_;
+  }
+
+  Planetarium MakePlanetarium(
+      Perspective<Navigation, Camera> const& perspective) const {
+    // No dark area, human visual acuity, wide field of view.
+    Planetarium::Parameters parameters(
+        /*sphere_radius_multiplier=*/1,
+        /*angular_resolution=*/0.4 * ArcMinute,
+        /*field_of_view=*/90 * Degree);
+    return Planetarium(parameters,
+                       perspective,
+                       ephemeris_.get(),
+                       earth_centred_inertial_.get());
+  }
+
+ private:
+  Ephemeris<Barycentric>::FixedStepParameters EphemerisParameters() {
+    return Ephemeris<Barycentric>::FixedStepParameters(
+        SymmetricLinearMultistepIntegrator<QuinlanTremaine1990Order12,
+                                           Position<Barycentric>>(),
+        /*step=*/10 * Minute);
+  }
+
+  Ephemeris<Barycentric>::FixedStepParameters HistoryParameters() {
+    return Ephemeris<Barycentric>::FixedStepParameters(
+        SymmetricLinearMultistepIntegrator<Quinlan1999Order8A,
+                                           Position<Barycentric>>(),
+        /*step=*/10 * Second);
+  }
+
+  not_null<std::unique_ptr<SolarSystem<Barycentric>>> const solar_system_;
+  not_null<std::unique_ptr<Ephemeris<Barycentric>>> const ephemeris_;
+  not_null<MassiveBody const*> const earth_;
+  not_null<std::unique_ptr<NavigationFrame>> const earth_centred_inertial_;
+  DiscreteTrajectory<Barycentric> goes_8_trajectory_;
+};
+
+}  // namespace
+
+void RunBenchmark(benchmark::State& state,
+                  Perspective<Navigation, Camera> const& perspective) {
+  Satellites satellites;
+  Planetarium planetarium = satellites.MakePlanetarium(perspective);
+  RP2Lines<Length, Camera> lines;
+  int total_lines = 0;
+  int iterations = 0;
+  // This is the time of a lunar eclipse in January 2000.
+  constexpr Instant now = "2000-01-21T04:41:30,5"_TT;
+  while (state.KeepRunning()) {
+    lines = planetarium.PlotMethod2(satellites.goes_8_trajectory().Begin(),
+                                    satellites.goes_8_trajectory().End(),
+                                    now,
+                                    /*reverse=*/false);
+    total_lines += lines.size();
+    ++iterations;
+  }
+  Length min_x = Infinity<Length>();
+  Length min_y = Infinity<Length>();
+  Length max_x = -Infinity<Length>();
+  Length max_y = -Infinity<Length>();
+  int points = 0;
+  for (auto const& line : lines) {
+    points += line.size();
+    for (auto const& point : line) {
+      min_x = std::min(min_x, point.x());
+      min_y = std::min(min_y, point.y());
+      max_x = std::max(max_x, point.x());
+      max_y = std::max(max_y, point.y());
+    }
+  }
+  state.SetLabel(std::to_string(points) + " points in " +
+                 std::to_string(total_lines / iterations) + " lines within [" +
+                 DebugString(min_x) + ", " + DebugString(max_x) + "] × [" +
+                 DebugString(min_y) + ", " + DebugString(max_y) + "]");
+}
+
+void BM_PlanetariumPlotMethod2NearPolarPerspective(benchmark::State& state) {
+  RunBenchmark(state, PolarPerspective(near));
+}
+
+void BM_PlanetariumPlotMethod2FarPolarPerspective(benchmark::State& state) {
+  RunBenchmark(state, PolarPerspective(far));
+}
+
+void BM_PlanetariumPlotMethod2NearEquatorialPerspective(
+    benchmark::State& state) {
+  RunBenchmark(state, EquatorialPerspective(near));
+}
+
+void BM_PlanetariumPlotMethod2FarEquatorialPerspective(
+    benchmark::State& state) {
+  RunBenchmark(state, EquatorialPerspective(far));
+}
+
+BENCHMARK(BM_PlanetariumPlotMethod2NearPolarPerspective);
+BENCHMARK(BM_PlanetariumPlotMethod2FarPolarPerspective);
+BENCHMARK(BM_PlanetariumPlotMethod2NearEquatorialPerspective);
+BENCHMARK(BM_PlanetariumPlotMethod2FarEquatorialPerspective);
+
+}  // namespace geometry
+}  // namespace principia