Merge pull request #3170 from pleroy/Redo

Convert mathematica, physics, astronomy and testing_utilities to DiscreteTraject0ry

Author: pleroy

astronomy/geodesy_test.cpp

@@ -7,6 +7,7 @@
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 #include "physics/body_surface_dynamic_frame.hpp"
+#include "physics/discrete_traject0ry.hpp"
 #include "physics/solar_system.hpp"
 #include "quantities/si.hpp"
 #include "testing_utilities/approximate_quantity.hpp"
@@ -32,7 +33,7 @@ using integrators::methods::QuinlanTremaine1990Order12;
 using physics::BodySurfaceDynamicFrame;
 using physics::ContinuousTrajectory;
 using physics::DegreesOfFreedom;
-using physics::DiscreteTrajectory;
+using physics::DiscreteTraject0ry;
 using physics::Ephemeris;
 using physics::KeplerianElements;
 using physics::KeplerOrbit;
@@ -125,15 +126,15 @@ TEST_F(GeodesyTest, DISABLED_LAGEOS2) {
 
   ephemeris_->Prolong(final_time);
 
-  DiscreteTrajectory<ICRS> primary_lageos2_trajectory;
+  DiscreteTraject0ry<ICRS> primary_lageos2_trajectory;
   primary_lageos2_trajectory.Append(
       initial_time, itrs_.FromThisFrameAtTime(initial_time)(initial_dof_ilrsa));
-  DiscreteTrajectory<ICRS> secondary_lageos2_trajectory;
+  DiscreteTraject0ry<ICRS> secondary_lageos2_trajectory;
   secondary_lageos2_trajectory.Append(
       initial_time, itrs_.FromThisFrameAtTime(initial_time)(initial_dof_ilrsb));
   auto flow_lageos2 =
       [this, final_time](
-          DiscreteTrajectory<ICRS>& lageos2_trajectory) -> absl::Status {
+          DiscreteTraject0ry<ICRS>& lageos2_trajectory) -> absl::Status {
         return ephemeris_->FlowWithAdaptiveStep(
             &lageos2_trajectory,
             Ephemeris<ICRS>::NoIntrinsicAcceleration,

astronomy/lunar_orbit_test.cpp

@@ -18,7 +18,7 @@
 #include "mathematica/mathematica.hpp"
 #include "physics/apsides.hpp"
 #include "physics/body_surface_dynamic_frame.hpp"
-#include "physics/discrete_trajectory.hpp"
+#include "physics/discrete_traject0ry.hpp"
 #include "physics/kepler_orbit.hpp"
 #include "physics/massless_body.hpp"
 #include "physics/oblate_body.hpp"
@@ -57,7 +57,7 @@ using physics::BodySurfaceDynamicFrame;
 using physics::ComputeApsides;
 using physics::ComputeNodes;
 using physics::DegreesOfFreedom;
-using physics::DiscreteTrajectory;
+using physics::DiscreteTraject0ry;
 using physics::Ephemeris;
 using physics::KeplerianElements;
 using physics::KeplerOrbit;
@@ -363,7 +363,7 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
                 IsNear(4.7e-13_⑴));
   }
 
-  DiscreteTrajectory<ICRS> trajectory;
+  DiscreteTraject0ry<ICRS> trajectory;
   trajectory.Append(J2000, initial_state);
   auto const instance = ephemeris_->NewInstance(
       {&trajectory},
@@ -377,7 +377,7 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
 
   // To find the nodes, we need to convert the trajectory to a reference frame
   // whose xy plane is the Moon's equator.
-  DiscreteTrajectory<LunarSurface> surface_trajectory;
+  DiscreteTraject0ry<LunarSurface> surface_trajectory;
   for (auto const& [time, degrees_of_freedom] : trajectory) {
     surface_trajectory.Append(
         time, lunar_frame_.ToThisFrameAtTime(time)(degrees_of_freedom));
@@ -413,18 +413,20 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
                   mathematica::ExpressIn(Metre));
   }
 
-  DiscreteTrajectory<LunarSurface> ascending_nodes;
-  DiscreteTrajectory<LunarSurface> descending_nodes;
-  ComputeNodes(surface_trajectory.begin(),
+  DiscreteTraject0ry<LunarSurface> ascending_nodes;
+  DiscreteTraject0ry<LunarSurface> descending_nodes;
+  ComputeNodes(surface_trajectory,
+               surface_trajectory.begin(),
                surface_trajectory.end(),
                /*north=*/Vector<double, LunarSurface>({0, 0, 1}),
                /*max_points=*/std::numeric_limits<int>::max(),
                ascending_nodes,
                descending_nodes);
 
-  DiscreteTrajectory<ICRS> apoapsides;
-  DiscreteTrajectory<ICRS> periapsides;
+  DiscreteTraject0ry<ICRS> apoapsides;
+  DiscreteTraject0ry<ICRS> periapsides;
   ComputeApsides(*ephemeris_->trajectory(moon_),
+                 trajectory,
                  trajectory.begin(),
                  trajectory.end(),
                  /*max_points=*/std::numeric_limits<int>::max(),
@@ -433,12 +435,12 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
 
   struct Nodes {
     std::string_view const name;
-    DiscreteTrajectory<LunarSurface> const& trajectory;
+    DiscreteTraject0ry<LunarSurface> const& trajectory;
   };
 
   struct Apsides {
     std::string_view const name;
-    DiscreteTrajectory<ICRS> const& trajectory;
+    DiscreteTraject0ry<ICRS> const& trajectory;
   };
 
   std::vector<double> descending_node_eccentricities;
@@ -530,7 +532,7 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
   {
     EccentricityVectorRange actual_period_ends;
     for (int orbit = 0;
-         orbit < descending_nodes.Size();
+         orbit < descending_nodes.size();
          orbit += orbits_per_period) {
       auto& actual = actual_period_ends;
       auto const e = descending_node_eccentricities[orbit];

astronomy/mercury_perihelion_test.cpp

@@ -12,7 +12,6 @@
 #include "integrators/symmetric_linear_multistep_integrator.hpp"
 #include "mathematica/mathematica.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"
@@ -33,7 +32,6 @@ using geometry::Position;
 using integrators::SymmetricLinearMultistepIntegrator;
 using integrators::methods::QuinlanTremaine1990Order12;
 using physics::ContinuousTrajectory;
-using physics::DiscreteTrajectory;
 using physics::Ephemeris;
 using physics::KeplerianElements;
 using physics::KeplerOrbit;

astronomy/orbit_analysis_test.cpp

@@ -14,6 +14,7 @@
 #include "mathematica/mathematica.hpp"
 #include "numerics/polynomial.hpp"
 #include "physics/body_centred_non_rotating_dynamic_frame.hpp"
+#include "physics/discrete_traject0ry.hpp"
 #include "physics/ephemeris.hpp"
 #include "physics/solar_system.hpp"
 #include "testing_utilities/approximate_quantity.hpp"
@@ -35,7 +36,7 @@ using numerics::EstrinEvaluator;
 using numerics::PolynomialInMonomialBasis;
 using physics::BodyCentredNonRotatingDynamicFrame;
 using physics::BodySurfaceDynamicFrame;
-using physics::DiscreteTrajectory;
+using physics::DiscreteTraject0ry;
 using physics::Ephemeris;
 using physics::MasslessBody;
 using physics::RotatingBody;
@@ -143,18 +144,18 @@ class OrbitAnalysisTest : public ::testing::Test {
 
   // Returns a GCRS trajectory obtained by stitching together the trajectories
   // of |sp3_orbit.satellites| in |sp3_orbit.files|.
-  not_null<std::unique_ptr<DiscreteTrajectory<GCRS>>> EarthCentredTrajectory(
+  not_null<std::unique_ptr<DiscreteTraject0ry<GCRS>>> EarthCentredTrajectory(
       SP3Orbit const& sp3_orbit) {
     BodyCentredNonRotatingDynamicFrame<ICRS, GCRS> gcrs{ephemeris_.get(),
                                                         &earth_};
     BodySurfaceDynamicFrame<ICRS, ITRS> itrs{ephemeris_.get(), &earth_};
 
-    auto result = make_not_null_unique<DiscreteTrajectory<GCRS>>();
+    auto result = make_not_null_unique<DiscreteTraject0ry<GCRS>>();
     for (auto const& file : sp3_orbit.files.names) {
       StandardProduct3 sp3(
           SOLUTION_DIR / "astronomy" / "standard_product_3" / file,
           sp3_orbit.files.dialect);
-      std::vector<not_null<DiscreteTrajectory<ITRS> const*>> const& orbit =
+      std::vector<not_null<DiscreteTraject0ry<ITRS> const*>> const& orbit =
           sp3.orbit(sp3_orbit.satellite);
       CHECK_EQ(orbit.size(), 1);
       auto const& arc = *orbit.front();

astronomy/orbit_ground_track.hpp

@@ -7,7 +7,7 @@
 #include "astronomy/orbit_recurrence.hpp"
 #include "astronomy/orbital_elements.hpp"
 #include "geometry/interval.hpp"
-#include "physics/discrete_trajectory.hpp"
+#include "physics/discrete_traject0ry.hpp"
 #include "physics/rotating_body.hpp"
 
 namespace principia {
@@ -16,7 +16,7 @@ namespace internal_orbit_ground_track {
 
 using geometry::Instant;
 using geometry::Interval;
-using physics::DiscreteTrajectory;
+using physics::DiscreteTraject0ry;
 using physics::RotatingBody;
 using quantities::Angle;
 using quantities::AngularFrequency;
@@ -67,7 +67,7 @@ class OrbitGroundTrack {
   // sun-synchronicity is analysed.
   template<typename PrimaryCentred, typename Inertial>
   static absl::StatusOr<OrbitGroundTrack> ForTrajectory(
-      DiscreteTrajectory<PrimaryCentred> const& trajectory,
+      DiscreteTraject0ry<PrimaryCentred> const& trajectory,
       RotatingBody<Inertial> const& primary,
       std::optional<MeanSun> const& mean_sun);
 

astronomy/orbit_ground_track_body.hpp

@@ -31,14 +31,14 @@ Angle CelestialLongitude(Position<PrimaryCentred> const& q) {
 // The resulting angles are neither normalized nor unwound.
 template<typename PrimaryCentred, typename Inertial>
 std::vector<Angle> PlanetocentricLongitudes(
-    DiscreteTrajectory<PrimaryCentred> const& nodes,
+    DiscreteTraject0ry<PrimaryCentred> const& nodes,
     RotatingBody<Inertial> const& primary) {
   std::vector<Angle> longitudes;
-  longitudes.reserve(nodes.Size());
-  for (auto const& node : nodes) {
+  longitudes.reserve(nodes.size());
+  for (auto const& [time, degrees_of_freedom] : nodes) {
     longitudes.push_back(
-        CelestialLongitude(node.degrees_of_freedom.position()) -
-        primary.AngleAt(node.time) - π / 2 * Radian);
+        CelestialLongitude(degrees_of_freedom.position()) -
+        primary.AngleAt(time) - π / 2 * Radian);
   }
   return longitudes;
 }
@@ -55,7 +55,7 @@ Angle MeanSolarTime(Iterator const& it,
 
 template<typename PrimaryCentred>
 Interval<Angle> MeanSolarTimesOfNodes(
-    DiscreteTrajectory<PrimaryCentred> const& nodes,
+    DiscreteTraject0ry<PrimaryCentred> const& nodes,
     OrbitGroundTrack::MeanSun const& mean_sun) {
   Interval<Angle> mean_solar_times;
   std::optional<Angle> mean_solar_time;
@@ -120,24 +120,25 @@ inline OrbitGroundTrack::EquatorCrossingLongitudes::EquatorCrossingLongitudes(
 
 template<typename PrimaryCentred, typename Inertial>
 absl::StatusOr<OrbitGroundTrack> OrbitGroundTrack::ForTrajectory(
-    DiscreteTrajectory<PrimaryCentred> const& trajectory,
+    DiscreteTraject0ry<PrimaryCentred> const& trajectory,
     RotatingBody<Inertial> const& primary,
     std::optional<MeanSun> const& mean_sun) {
-  DiscreteTrajectory<PrimaryCentred> ascending_nodes;
-  DiscreteTrajectory<PrimaryCentred> descending_nodes;
+  DiscreteTraject0ry<PrimaryCentred> ascending_nodes;
+  DiscreteTraject0ry<PrimaryCentred> descending_nodes;
   OrbitGroundTrack ground_track;
-  RETURN_IF_ERROR(ComputeNodes(trajectory.begin(),
+  RETURN_IF_ERROR(ComputeNodes(trajectory,
+                               trajectory.begin(),
                                trajectory.end(),
                                Vector<double, PrimaryCentred>({0, 0, 1}),
                                /*max_points=*/std::numeric_limits<int>::max(),
                                ascending_nodes,
                                descending_nodes));
   if (mean_sun.has_value()) {
-    if (!ascending_nodes.Empty()) {
+    if (!ascending_nodes.empty()) {
       ground_track.mean_solar_times_of_ascending_nodes_ =
           MeanSolarTimesOfNodes(ascending_nodes, *mean_sun);
     }
-    if (!descending_nodes.Empty()) {
+    if (!descending_nodes.empty()) {
       ground_track.mean_solar_times_of_descending_nodes_ =
           MeanSolarTimesOfNodes(descending_nodes, *mean_sun);
     }
@@ -147,7 +148,7 @@ absl::StatusOr<OrbitGroundTrack> OrbitGroundTrack::ForTrajectory(
   ground_track.longitudes_of_equator_crossings_of_descending_passes_ =
       PlanetocentricLongitudes(descending_nodes, primary);
   ground_track.first_descending_pass_before_first_ascending_pass_ =
-      !ascending_nodes.Empty() && !descending_nodes.Empty() &&
+      !ascending_nodes.empty() && !descending_nodes.empty() &&
       descending_nodes.front().time < ascending_nodes.front().time;
   return ground_track;
 }

astronomy/orbital_elements.hpp

@@ -6,7 +6,7 @@
 #include "geometry/interval.hpp"
 #include "geometry/named_quantities.hpp"
 #include "physics/body.hpp"
-#include "physics/discrete_trajectory.hpp"
+#include "physics/discrete_traject0ry.hpp"
 #include "physics/massive_body.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/quantities.hpp"
@@ -18,7 +18,7 @@ namespace internal_orbital_elements {
 using geometry::Instant;
 using geometry::Interval;
 using physics::Body;
-using physics::DiscreteTrajectory;
+using physics::DiscreteTraject0ry;
 using physics::MassiveBody;
 using quantities::Angle;
 using quantities::AngularFrequency;
@@ -37,7 +37,7 @@ class OrbitalElements {
 
   template<typename PrimaryCentred>
   static absl::StatusOr<OrbitalElements> ForTrajectory(
-      DiscreteTrajectory<PrimaryCentred> const& trajectory,
+      DiscreteTraject0ry<PrimaryCentred> const& trajectory,
       MassiveBody const& primary,
       Body const& secondary);
 
@@ -146,13 +146,13 @@ class OrbitalElements {
 
   template<typename PrimaryCentred>
   static std::vector<EquinoctialElements> OsculatingEquinoctialElements(
-      DiscreteTrajectory<PrimaryCentred> const& trajectory,
+      DiscreteTraject0ry<PrimaryCentred> const& trajectory,
       MassiveBody const& primary,
       Body const& secondary);
 
   template<typename PrimaryCentred>
   static std::vector<Length> RadialDistances(
-      DiscreteTrajectory<PrimaryCentred> const& trajectory);
+      DiscreteTraject0ry<PrimaryCentred> const& trajectory);
 
   // |equinoctial_elements| must contain at least 2 elements.
   static absl::StatusOr<Time> SiderealPeriod(

astronomy/orbital_elements_body.hpp

@@ -5,6 +5,7 @@
 #include <algorithm>
 #include <vector>
 
+#include "absl/strings/str_cat.h"
 #include "base/jthread.hpp"
 #include "base/status_utilities.hpp"
 #include "physics/kepler_orbit.hpp"
@@ -32,13 +33,13 @@ using quantities::si::Radian;
 
 template<typename PrimaryCentred>
 absl::StatusOr<OrbitalElements> OrbitalElements::ForTrajectory(
-    DiscreteTrajectory<PrimaryCentred> const& trajectory,
+    DiscreteTraject0ry<PrimaryCentred> const& trajectory,
     MassiveBody const& primary,
     Body const& secondary) {
   OrbitalElements orbital_elements;
-  if (trajectory.Size() < 2) {
+  if (trajectory.size() < 2) {
     return absl::InvalidArgumentError(
-        "trajectory.Size() is " + std::to_string(trajectory.Size()));
+        absl::StrCat("trajectory.Size() is ", trajectory.size()));
   }
   orbital_elements.osculating_equinoctial_elements_ =
       OsculatingEquinoctialElements(trajectory, primary, secondary);
@@ -138,13 +139,13 @@ inline Interval<Length> OrbitalElements::radial_distance_interval() const {
 template<typename PrimaryCentred>
 std::vector<OrbitalElements::EquinoctialElements>
 OrbitalElements::OsculatingEquinoctialElements(
-    DiscreteTrajectory<PrimaryCentred> const& trajectory,
+    DiscreteTraject0ry<PrimaryCentred> const& trajectory,
     MassiveBody const& primary,
     Body const& secondary) {
   DegreesOfFreedom<PrimaryCentred> const primary_dof{
       PrimaryCentred::origin, PrimaryCentred::unmoving};
   std::vector<EquinoctialElements> result;
-  result.reserve(trajectory.Size());
+  result.reserve(trajectory.size());
   for (auto const& [time, degrees_of_freedom] : trajectory) {
     auto const osculating_elements =
         KeplerOrbit<PrimaryCentred>(primary,
@@ -175,9 +176,9 @@ OrbitalElements::OsculatingEquinoctialElements(
 
 template<typename PrimaryCentred>
 std::vector<Length> OrbitalElements::RadialDistances(
-    DiscreteTrajectory<PrimaryCentred> const& trajectory) {
+    DiscreteTraject0ry<PrimaryCentred> const& trajectory) {
   std::vector<Length> radial_distances;
-  radial_distances.reserve(trajectory.Size());
+  radial_distances.reserve(trajectory.size());
   DegreesOfFreedom<PrimaryCentred> const primary_dof{PrimaryCentred::origin,
                                                      PrimaryCentred::unmoving};
   for (auto const& [time, degrees_of_freedom] : trajectory) {