Merge pull request #4213 from pleroy/ANewHope

PlotMethod4: A New Hope

Author: pleroy

Principia.sln.DotSettings

@@ -221,6 +221,7 @@
 	<s:Boolean x:Key="/Default/UserDictionary/Words/=colour/@EntryIndexedValue">True</s:Boolean>
 	<s:Boolean x:Key="/Default/UserDictionary/Words/=colours/@EntryIndexedValue">True</s:Boolean>
 	<s:Boolean x:Key="/Default/UserDictionary/Words/=downsampling/@EntryIndexedValue">True</s:Boolean>
+	<s:Boolean x:Key="/Default/UserDictionary/Words/=equipotential/@EntryIndexedValue">True</s:Boolean>
 	<s:Boolean x:Key="/Default/UserDictionary/Words/=Frenet/@EntryIndexedValue">True</s:Boolean>
 	<s:Boolean x:Key="/Default/UserDictionary/Words/=geopotential/@EntryIndexedValue">True</s:Boolean>
 	<s:Boolean x:Key="/Default/UserDictionary/Words/=gipfeli/@EntryIndexedValue">True</s:Boolean>

benchmarks/planetarium_benchmark.cpp

@@ -220,12 +220,13 @@ class Satellites {
   }
 
   Planetarium MakePlanetarium(
+      Angle const& angular_resolution,
       Perspective<Navigation, Camera> const& perspective,
       not_null<PlottingFrame const*> plotting_frame) const {
     // No dark area, human visual acuity, wide field of view.
     Planetarium::Parameters parameters(
         /*sphere_radius_multiplier=*/1,
-        /*angular_resolution=*/0.4 * ArcMinute,
+        angular_resolution,
         /*field_of_view=*/90 * Degree);
     Instant const t = goes_8_trajectory().front().time;
     Similarity<Navigation, GCRS> const plotting_to_gcrs =
@@ -316,6 +317,7 @@ void BM_PlanetariumPlotMethod3(
   Instant const t = satellites.goes_8_trajectory().front().time;
   PlottingFrame const& plotting = (satellites.*plotting_frame)();
   Planetarium planetarium = satellites.MakePlanetarium(
+      0.4 * ArcMinute,
       perspective((satellites.gcrs().ToThisFrameAtTimeSimilarly(t) *
                    plotting.FromThisFrameAtTimeSimilarly(t)).similarity(),
                   distance_from_earth),
@@ -350,6 +352,52 @@ void BM_PlanetariumPlotMethod3(
                      .str());
 }
 
+void BM_PlanetariumPlotMethod4(
+    benchmark::State& state,
+    Perspective<Navigation, Camera> (*const perspective)(
+        Similarity<Navigation, GCRS> const& navigation_to_gcrs_at_epoch,
+        Length const distance_from_earth),
+    Length const distance_from_earth,
+    PlottingFrame const& (Satellites::*const plotting_frame)() const) {
+  static Satellites satellites;
+  Instant const t = satellites.goes_8_trajectory().front().time;
+  PlottingFrame const& plotting = (satellites.*plotting_frame)();
+  Planetarium planetarium = satellites.MakePlanetarium(
+      1 * ArcMinute,
+      perspective((satellites.gcrs().ToThisFrameAtTimeSimilarly(t) *
+                   plotting.FromThisFrameAtTimeSimilarly(t)).similarity(),
+                  distance_from_earth),
+      &plotting);
+  std::vector<ScaledSpacePoint> line;
+  int iterations = 0;
+  // This is the time of a lunar eclipse in January 2000.
+  constexpr Instant now = "2000-01-21T04:41:30,5"_TT;
+  for (auto _ : state) {
+    line.clear();
+    planetarium.PlotMethod4(
+        satellites.goes_8_trajectory(),
+        satellites.goes_8_trajectory().begin(),
+        satellites.goes_8_trajectory().end(),
+        /*t_max=*/InfiniteFuture,
+        /*reverse=*/false,
+        /*add_point=*/
+        [&line](ScaledSpacePoint const& point) { line.push_back(point); },
+        /*max_points=*/std::numeric_limits<int>::max());
+    ++iterations;
+  }
+  Interval<double> x;
+  Interval<double> y;
+  Interval<double> z;
+  for (auto const& point : line) {
+    x.Include(point.x);
+    y.Include(point.y);
+    z.Include(point.z);
+  }
+  state.SetLabel((std::stringstream() << line.size() << " points within " << x
+                                      << " × " << y << " × " << z)
+                     .str());
+}
+
 #define PRINCIPIA_BENCHMARK_PLANETARIUM_PLOT_METHODS_NEAR_AND_FAR( \
     name, perspective, plotting_frame)                             \
   BENCHMARK_CAPTURE(BM_PlanetariumPlotMethod3,                     \
@@ -358,11 +406,23 @@ void BM_PlanetariumPlotMethod3(
                     near,                                          \
                     (plotting_frame))                              \
       ->Unit(benchmark::kMillisecond);                             \
+  BENCHMARK_CAPTURE(BM_PlanetariumPlotMethod4,                     \
+                    Near##name,                                    \
+                    (perspective),                                 \
+                    near,                                          \
+                    (plotting_frame))                              \
+      ->Unit(benchmark::kMillisecond);                             \
   BENCHMARK_CAPTURE(BM_PlanetariumPlotMethod3,                     \
                     Far##name,                                     \
                     (perspective),                                 \
                     far,                                           \
                     (plotting_frame))                              \
+      ->Unit(benchmark::kMillisecond);                             \
+  BENCHMARK_CAPTURE(BM_PlanetariumPlotMethod4,                     \
+                    Far##name,                                     \
+                    (perspective),                                 \
+                    far,                                           \
+                    (plotting_frame))                              \
       ->Unit(benchmark::kMillisecond)
 
 #define PRINCIPIA_BENCHMARK_PLANETARIUM_PLOT_METHODS_POLAR_AND_EQUATORIAL( \

geometry/perspective.hpp

@@ -43,6 +43,7 @@ class Perspective final {
   Perspective(Similarity<FromFrame, ToFrame> const& to_camera,
               Length const& focal);
 
+  Position<FromFrame> const& camera() const;
   Length const& focal() const;
 
   // Returns the ℝP² element resulting from the projection of `point`.  This

geometry/perspective_body.hpp

@@ -42,6 +42,12 @@ Perspective<FromFrame, ToFrame>::Perspective(
       camera_(from_camera_(ToFrame::origin)),
       focal_(focal) {}
 
+template<typename FromFrame, typename ToFrame>
+Position<FromFrame> const&
+Perspective<FromFrame, ToFrame>::camera() const {
+  return camera_;
+}
+
 template<typename FromFrame, typename ToFrame>
 Length const& Perspective<FromFrame, ToFrame>::focal() const {
   return focal_;

ksp_plugin/interface_planetarium.cpp

@@ -50,6 +50,7 @@ Planetarium* __cdecl principia__PlanetariumCreate(
     double const focal,
     double const field_of_view,
     double const inverse_scale_factor,
+    double const angular_resolution,
     XYZ const scaled_space_origin) {
   journal::Method<journal::PlanetariumCreate> m({plugin,
                                                  sun_world_position,
@@ -60,6 +61,7 @@ Planetarium* __cdecl principia__PlanetariumCreate(
                                                  focal,
                                                  field_of_view,
                                                  inverse_scale_factor,
+                                                 angular_resolution,
                                                  scaled_space_origin});
   Renderer const& renderer = CHECK_NOTNULL(plugin)->renderer();
 
@@ -90,7 +92,7 @@ Planetarium* __cdecl principia__PlanetariumCreate(
 
   Planetarium::Parameters parameters(
       /*sphere_radius_multiplier=*/1.0,
-      /*angular_resolution=*/0.4 * ArcMinute,
+      angular_resolution * Radian,
       field_of_view * Radian);
   Perspective<Navigation, Camera> perspective(
       world_to_plotting_affine_map *
@@ -151,7 +153,7 @@ void __cdecl principia__PlanetariumPlotFlightPlanSegment(
   if (index % 2 == 0 ||
       segment->empty() ||
       segment->front().time >= plugin->renderer().GetPlottingFrame()->t_min()) {
-    planetarium->PlotMethod3(
+    planetarium->PlotMethod4(
         *segment,
         segment->begin(),
         segment->end(),
@@ -183,7 +185,7 @@ void __cdecl principia__PlanetariumPlotPrediction(
   *vertex_count = 0;
 
   auto const prediction = plugin->GetVessel(vessel_guid)->prediction();
-  planetarium->PlotMethod3(
+  planetarium->PlotMethod4(
       *prediction,
       prediction->begin(),
       prediction->end(),
@@ -239,7 +241,7 @@ void __cdecl principia__PlanetariumPlotPsychohistory(
     // time the history will be shorter than desired.
     vessel->RequestReanimation(desired_first_time);
 
-    planetarium->PlotMethod3(
+    planetarium->PlotMethod4(
         trajectory,
         trajectory.lower_bound(desired_first_time),
         psychohistory->end(),
@@ -296,7 +298,7 @@ void __cdecl principia__PlanetariumPlotCelestialPastTrajectory(
     Instant const first_time =
         std::max(desired_first_time, celestial_trajectory.t_min());
     Length minimal_distance;
-    planetarium->PlotMethod3(
+    planetarium->PlotMethod4(
         celestial_trajectory,
         first_time,
         /*last_time=*/plugin->CurrentTime(),
@@ -354,7 +356,7 @@ void __cdecl principia__PlanetariumPlotCelestialFutureTrajectory(
     // No need to request reanimation here because the current time of the
     // plugin is necessarily covered.
     Length minimal_distance;
-    planetarium->PlotMethod3(
+    planetarium->PlotMethod4(
         celestial_trajectory,
         /*first_time=*/plugin->CurrentTime(),
         /*last_time=*/final_time,
@@ -392,7 +394,7 @@ void __cdecl principia__PlanetariumPlotEquipotential(
   DiscreteTrajectory<Navigation> const& equipotential =
       equipotentials.lines[index];
 
-  planetarium->PlotMethod3(
+  planetarium->PlotMethod4(
       equipotential,
       equipotential.front().time,
       equipotential.back().time,

ksp_plugin/planetarium.cpp

@@ -30,6 +30,7 @@ Planetarium::Parameters::Parameters(double const sphere_radius_multiplier,
                                     Angle const& angular_resolution,
                                     Angle const& field_of_view)
     : sphere_radius_multiplier_(sphere_radius_multiplier),
+      angular_resolution_(angular_resolution),
       sin²_angular_resolution_(Pow<2>(Sin(angular_resolution))),
       tan_angular_resolution_(Tan(angular_resolution)),
       tan_field_of_view_(Tan(field_of_view)) {}
@@ -284,6 +285,25 @@ void Planetarium::PlotMethod3(
       trajectory, begin_time, last_time, reverse, add_point, max_points);
 }
 
+void Planetarium::PlotMethod4(
+    Trajectory<Barycentric> const& trajectory,
+    DiscreteTrajectory<Barycentric>::iterator begin,
+    DiscreteTrajectory<Barycentric>::iterator end,
+    Instant const& t_max,
+    bool const reverse,
+    std::function<void(ScaledSpacePoint const&)> const& add_point,
+    int max_points) const {
+  if (begin == end) {
+    return;
+  }
+  auto last = std::prev(end);
+  auto const begin_time = std::max(begin->time, plotting_frame_->t_min());
+  auto const last_time =
+      std::min({last->time, plotting_frame_->t_max(), t_max});
+  PlotMethod4(
+      trajectory, begin_time, last_time, reverse, add_point, max_points);
+}
+
 std::vector<Sphere<Navigation>> Planetarium::ComputePlottableSpheres(
     Instant const& now) const {
   SimilarMotion<Barycentric, Navigation> const similar_motion_at_now =

ksp_plugin/planetarium.hpp

@@ -4,7 +4,6 @@
 
 #include "base/not_null.hpp"
 #include "geometry/instant.hpp"
-#include "geometry/orthogonal_map.hpp"
 #include "geometry/perspective.hpp"
 #include "geometry/r3_element.hpp"
 #include "geometry/rp2_point.hpp"
@@ -14,7 +13,6 @@
 #include "physics/degrees_of_freedom.hpp"
 #include "physics/discrete_trajectory.hpp"
 #include "physics/ephemeris.hpp"
-#include "physics/rigid_motion.hpp"
 #include "physics/trajectory.hpp"
 #include "quantities/quantities.hpp"
 
@@ -25,7 +23,6 @@ namespace internal {
 
 using namespace principia::base::_not_null;
 using namespace principia::geometry::_instant;
-using namespace principia::geometry::_orthogonal_map;
 using namespace principia::geometry::_perspective;
 using namespace principia::geometry::_r3_element;
 using namespace principia::geometry::_rp2_point;
@@ -35,7 +32,6 @@ using namespace principia::ksp_plugin::_frames;
 using namespace principia::physics::_degrees_of_freedom;
 using namespace principia::physics::_discrete_trajectory;
 using namespace principia::physics::_ephemeris;
-using namespace principia::physics::_rigid_motion;
 using namespace principia::physics::_trajectory;
 using namespace principia::quantities::_quantities;
 
@@ -69,6 +65,7 @@ class Planetarium {
 
    private:
     double const sphere_radius_multiplier_;
+    Angle const angular_resolution_;
     double const sin²_angular_resolution_;
     double const tan_angular_resolution_;
     double const tan_field_of_view_;
@@ -125,7 +122,8 @@ class Planetarium {
       Length* minimal_distance = nullptr) const;
 
   // A method similar to PlotMethod2, but which produces a three-dimensional
-  // trajectory in scaled space instead of projecting and hiding.
+  // trajectory in scaled space instead of projecting and hiding.  It uses the
+  // apparent angle of the sagitta as the metric to analyse curvature.
   void PlotMethod3(
       Trajectory<Barycentric> const& trajectory,
       DiscreteTrajectory<Barycentric>::iterator begin,
@@ -147,6 +145,29 @@ class Planetarium {
       int max_points,
       Length* minimal_distance = nullptr) const;
 
+  // A method similar to PlotMethod4, but which uses the RMS of the apparent
+  // distance between the trajectory and line segments.
+  void PlotMethod4(
+      Trajectory<Barycentric> const& trajectory,
+      DiscreteTrajectory<Barycentric>::iterator begin,
+      DiscreteTrajectory<Barycentric>::iterator end,
+      Instant const& t_max,
+      bool reverse,
+      std::function<void(ScaledSpacePoint const&)> const& add_point,
+      int max_points) const;
+
+  // The same method, operating on the `Trajectory` interface for any frame that
+  // can be converted to `Navigation`.
+  template<typename Frame>
+  void PlotMethod4(
+      Trajectory<Frame> const& trajectory,
+      Instant const& first_time,
+      Instant const& last_time,
+      bool reverse,
+      std::function<void(ScaledSpacePoint const&)> const& add_point,
+      int max_points,
+      Length* minimal_distance = nullptr) const;
+
  private:
   // Computes the coordinates of the spheres that represent the `ephemeris_`
   // bodies.  These coordinates are in the `plotting_frame_` at time `now`.
@@ -160,6 +181,12 @@ class Planetarium {
       DiscreteTrajectory<Barycentric>::iterator begin,
       DiscreteTrajectory<Barycentric>::iterator end) const;
 
+  // Computes the proper motion (in the astronomical sense) of the given point
+  // and velocity seen from a frame centered at the origin of `Camera`.  This
+  // is a uniform rotation on a great circle.
+  AngularVelocity<Navigation> ProperMotion(
+      DegreesOfFreedom<Navigation> const& degrees_of_freedom) const;
+
   Parameters const parameters_;
   Perspective<Navigation, Camera> const perspective_;
   not_null<Ephemeris<Barycentric> const*> const ephemeris_;