A function to compute a collision between a vessel and a celestial

geometry/r3_element.hpp

@@ -94,6 +94,11 @@ SphericalCoordinates<Scalar> RadiusLatitudeLongitude(Scalar const& radius,
                                                      Angle const& latitude,
                                                      Angle const& longitude);
 
+template<typename Scalar>
+std::ostream& operator<<(
+    std::ostream& out,
+    SphericalCoordinates<Scalar> const& spherical_coordinates);
+
 template<typename Scalar>
 R3Element<Scalar> operator+(R3Element<Scalar> const& right);
 template<typename Scalar>

geometry/r3_element_body.hpp

@@ -206,6 +206,16 @@ SphericalCoordinates<Scalar> RadiusLatitudeLongitude(Scalar const& radius,
   return result;
 }
 
+template<typename Scalar>
+std::ostream& operator<<(
+    std::ostream& out,
+    SphericalCoordinates<Scalar> const& spherical_coordinates) {
+  out << spherical_coordinates.radius
+      << ", " << spherical_coordinates.latitude
+      << ", " << spherical_coordinates.longitude;
+  return out;
+}
+
 template<typename Scalar>
 R3Element<Scalar> operator+(R3Element<Scalar> const& right) {
   return R3Element<Scalar>(+right.x, +right.y, +right.z);

physics/apsides.hpp

@@ -6,7 +6,9 @@
 #include "base/constant_function.hpp"
 #include "geometry/grassmann.hpp"
 #include "physics/discrete_trajectory.hpp"
+#include "physics/rotating_body.hpp"
 #include "physics/trajectory.hpp"
+#include "quantities/quantities.hpp"
 
 namespace principia {
 namespace physics {
@@ -16,7 +18,9 @@ namespace internal {
 using namespace principia::base::_constant_function;
 using namespace principia::geometry::_grassmann;
 using namespace principia::physics::_discrete_trajectory;
+using namespace principia::physics::_rotating_body;
 using namespace principia::physics::_trajectory;
+using namespace principia::quantities::_quantities;
 
 // Computes the apsides with respect to |reference| for the section given by
 // |begin| and |end| of |trajectory|.  Appends to the given output trajectories
@@ -30,6 +34,21 @@ void ComputeApsides(Trajectory<Frame> const& reference,
                     DiscreteTrajectory<Frame>& apoapsides,
                     DiscreteTrajectory<Frame>& periapsides);
 
+// Computes a collision between a vessel and a rotating body.  |begin| and |end|
+// must be on opposite sides of the surface of the body (in particular, a
+// collision must exist).  |radius| gives the radius of the celestial at a
+// particular position given by its latitude and longitude.  It must never
+// exceed the |max_radius| of the body.
+template<typename Frame>
+typename DiscreteTrajectory<Frame>::value_type ComputeCollision(
+    RotatingBody<Frame> const& reference_body,
+    Trajectory<Frame> const& reference,
+    Trajectory<Frame> const& trajectory,
+    typename DiscreteTrajectory<Frame>::iterator begin,
+    typename DiscreteTrajectory<Frame>::iterator end,
+    std::function<Length(Angle const& latitude,
+                         Angle const& longitude)> const& radius);
+
 // Computes the crossings of the section given by |begin| and |end| of
 // |trajectory| with the xy plane.  Appends the crossings that go towards the
 // |north| side of the xy plane to |ascending|, and those that go away from the
@@ -60,6 +79,7 @@ void ComputeApsides(Trajectory<Frame> const& trajectory1,
 }  // namespace internal
 
 using internal::ComputeApsides;
+using internal::ComputeCollision;
 using internal::ComputeNodes;
 
 }  // namespace _apsides

physics/apsides_body.hpp

@@ -3,17 +3,17 @@
 #include "physics/apsides.hpp"
 
 #include <optional>
-#include <vector>
 
 #include "base/array.hpp"
 #include "geometry/barycentre_calculator.hpp"
 #include "geometry/instant.hpp"
+#include "geometry/r3_element.hpp"
 #include "geometry/sign.hpp"
+#include "geometry/space.hpp"
 #include "numerics/hermite3.hpp"
 #include "numerics/root_finders.hpp"
 #include "physics/degrees_of_freedom.hpp"
 #include "quantities/named_quantities.hpp"
-#include "quantities/quantities.hpp"
 
 namespace principia {
 namespace physics {
@@ -23,12 +23,13 @@ namespace internal {
 using namespace principia::base::_array;
 using namespace principia::geometry::_barycentre_calculator;
 using namespace principia::geometry::_instant;
+using namespace principia::geometry::_r3_element;
 using namespace principia::geometry::_sign;
+using namespace principia::geometry::_space;
 using namespace principia::numerics::_hermite3;
 using namespace principia::numerics::_root_finders;
 using namespace principia::physics::_degrees_of_freedom;
 using namespace principia::quantities::_named_quantities;
-using namespace principia::quantities::_quantities;
 
 template<typename Frame>
 void ComputeApsides(Trajectory<Frame> const& reference,
@@ -134,6 +135,72 @@ void ComputeApsides(Trajectory<Frame> const& reference,
   }
 }
 
+template<typename Frame>
+typename DiscreteTrajectory<Frame>::value_type ComputeCollision(
+    RotatingBody<Frame> const& reference_body,
+    Trajectory<Frame> const& reference,
+    Trajectory<Frame> const& trajectory,
+    typename DiscreteTrajectory<Frame>::iterator const begin,
+    typename DiscreteTrajectory<Frame>::iterator const end,
+    std::function<Length(Angle const& latitude,
+                         Angle const& longitude)> const& radius) {
+  // The frame of the surface of the celestial.
+  using SurfaceFrame = geometry::_frame::Frame<struct SurfaceFrameTag>;
+
+  auto const last = std::prev(end);
+  Square<Length> max_radius² = Pow<2>(reference_body.max_radius());
+
+  std::optional<typename DiscreteTrajectory<Frame>::iterator>
+      last_above_max_radius;
+  for (auto it = begin; it != end; ++it) {
+    auto const& [time, degrees_of_freedom] = *it;
+    DegreesOfFreedom<Frame> const body_degrees_of_freedom =
+        reference.EvaluateDegreesOfFreedom(time);
+    RelativeDegreesOfFreedom<Frame> const relative =
+        degrees_of_freedom - body_degrees_of_freedom;
+    Square<Length> const squared_distance = relative.displacement().Norm²();
+    if (squared_distance <= max_radius²) {
+      break;
+    } else {
+      last_above_max_radius = it;
+    }
+  }
+
+  auto radius_at_time = [&radius,
+                         &reference,
+                         &reference_body,
+                         &trajectory](Instant const& t) {
+    auto const reference_position = reference.EvaluatePosition(t);
+    auto const trajectory_position = trajectory.EvaluatePosition(t);
+    Displacement<Frame> const displacement_in_frame =
+        trajectory_position - reference_position;
+
+    auto const to_surface_frame =
+        reference_body.ToSurfaceFrame<SurfaceFrame>(t);
+    Displacement<SurfaceFrame> const displacement_in_surface =
+        to_surface_frame(displacement_in_frame);
+
+    SphericalCoordinates<Length> const spherical_coordinates =
+        displacement_in_surface.coordinates().ToSpherical();
+
+    return spherical_coordinates.radius -
+           radius(spherical_coordinates.latitude,
+                  spherical_coordinates.longitude);
+  };
+
+  CHECK_LT(Length{}, radius_at_time(begin->time));
+  CHECK_LT(radius_at_time(last->time), Length{});
+
+  Instant const collision_time =
+      Brent(radius_at_time,
+            last_above_max_radius.value_or(begin)->time,
+            last->time);
+
+  return typename DiscreteTrajectory<Frame>::value_type(
+      collision_time,
+      trajectory.EvaluateDegreesOfFreedom(collision_time));
+}
+
 template<typename Frame, typename Predicate>
 absl::Status ComputeNodes(
     Trajectory<Frame> const& trajectory,

physics/apsides_test.cpp

@@ -23,12 +23,17 @@
 #include "physics/kepler_orbit.hpp"
 #include "physics/massive_body.hpp"
 #include "physics/massless_body.hpp"
+#include "physics/rotating_body.hpp"
 #include "quantities/astronomy.hpp"
 #include "quantities/elementary_functions.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/quantities.hpp"
 #include "quantities/si.hpp"
 #include "testing_utilities/almost_equals.hpp"
+#include "testing_utilities/approximate_quantity.hpp"
+#include "testing_utilities/componentwise.hpp"
+#include "testing_utilities/discrete_trajectory_factories.hpp"
+#include "testing_utilities/is_near.hpp"
 #include "testing_utilities/matchers.hpp"  // 🧙 For EXPECT_OK.
 
 namespace principia {
@@ -51,12 +56,17 @@ using namespace principia::physics::_ephemeris;
 using namespace principia::physics::_kepler_orbit;
 using namespace principia::physics::_massive_body;
 using namespace principia::physics::_massless_body;
+using namespace principia::physics::_rotating_body;
 using namespace principia::quantities::_astronomy;
 using namespace principia::quantities::_elementary_functions;
 using namespace principia::quantities::_named_quantities;
 using namespace principia::quantities::_quantities;
 using namespace principia::quantities::_si;
+using namespace principia::testing_utilities::_approximate_quantity;
 using namespace principia::testing_utilities::_almost_equals;
+using namespace principia::testing_utilities::_componentwise;
+using namespace principia::testing_utilities::_discrete_trajectory_factories;
+using namespace principia::testing_utilities::_is_near;
 
 class ApsidesTest : public ::testing::Test {
  protected:
@@ -162,6 +172,78 @@ TEST_F(ApsidesTest, ComputeApsidesDiscreteTrajectory) {
   }
 }
 
+TEST_F(ApsidesTest, ComputeCollision) {
+  Instant const t0;
+  DiscreteTrajectory<World> reference_trajectory;
+  DiscreteTrajectory<World> vessel_trajectory;
+
+  // At |t0| the vessel is inside the celestial, so we expect the collision at a
+  // negative time.
+  AppendTrajectoryTimeline(
+      NewLinearTrajectoryTimeline(
+          DegreesOfFreedom<World>(
+              World::origin +
+                  Displacement<World>({0 * Metre, -1 * Metre, 0 * Metre}),
+              Velocity<World>({1 * Metre / Second,
+                               0 * Metre / Second,
+                               0 * Metre / Second})),
+          /*Δt=*/1 * Second,
+          t0,
+          /*t1=*/t0 - 10 * Second,
+          /*t2=*/t0 + 10 * Second),
+      reference_trajectory);
+  AppendTrajectoryTimeline(
+      NewLinearTrajectoryTimeline(
+          DegreesOfFreedom<World>(
+              World::origin +
+                  Displacement<World>({1 * Metre, -1 * Metre, 0 * Metre}),
+              Velocity<World>({0 * Metre / Second,
+                               -1 * Metre / Second,
+                               0 * Metre / Second})),
+          /*Δt=*/1 * Second,
+          t0,
+          /*t1=*/t0 - 10 * Second,
+          /*t2=*/t0 + 1 * Second),
+      vessel_trajectory);
+
+  RotatingBody<World> const body(
+      1 * Kilogram,
+      RotatingBody<World>::Parameters(
+          /*min_radius=*/1 * Metre,
+          /*mean_radius=*/2 * Metre,
+          /*max_radius=*/3 * Metre,
+          /*reference_angle=*/0 * Radian,
+          /*reference_instant=*/t0,
+          /*angular_frequency=*/2 * π * Radian / Second,
+          /*right_ascension_of_pole=*/0 * Radian,
+          /*declination_of_pole=*/π / 2 * Radian));
+
+  // The celestial is infinite in the z direction and has four lobes in the x-y
+  // plane.  Think of a LEGO® axle.
+  auto radius = [](Angle const& latitude, Angle const& longitude) {
+    return (Cos(4 * longitude) + 2) * Metre;
+  };
+
+  auto const collision = ComputeCollision(body,
+                                          reference_trajectory,
+                                          vessel_trajectory,
+                                          vessel_trajectory.begin(),
+                                          vessel_trajectory.end(),
+                                          radius);
+
+  // The collision was verified with Mathematica to the given accuracy.
+  EXPECT_THAT(collision.time - t0, IsNear(-1.43861971643135_(1) * Second));
+  EXPECT_THAT(collision.degrees_of_freedom.position() - World::origin,
+              Componentwise(1 * Metre,
+                            IsNear(0.43861971643135_(1) * Metre),
+                            0 * Metre));
+  EXPECT_THAT(
+      collision.degrees_of_freedom.velocity(),
+      AlmostEquals(Velocity<World>({0 * Metre / Second,
+                                    -1 * Metre / Second,
+                                    0 * Metre / Second}), 0));
+}
+
 TEST_F(ApsidesTest, ComputeNodes) {
   Instant const t0;
   GravitationalParameter const μ = SolarGravitationalParameter;

physics/rotating_body.hpp

@@ -139,9 +139,6 @@ class RotatingBody : public MassiveBody {
   template<typename CelestialFrame>
   Rotation<Frame, CelestialFrame> ToCelestialFrame() const;
 
-  // Returns the rotation at time |t|.
-  Rotation<Frame, Frame> RotationAt(Instant const& t) const;
-
   // Returns false.
   bool is_massless() const override;
 

physics/rotating_body_body.hpp

@@ -143,11 +143,6 @@ Angle RotatingBody<Frame>::AngleAt(Instant const& t) const {
          (t - parameters_.reference_instant_) * parameters_.angular_frequency_;
 }
 
-template<typename Frame>
-Rotation<Frame, Frame> RotatingBody<Frame>::RotationAt(Instant const& t) const {
-  return Exp((t - parameters_.reference_instant_) * angular_velocity_);
-}
-
 template<typename Frame>
 bool RotatingBody<Frame>::is_massless() const {
   return false;