text: ﻿
  #include "physics/dynamic_frame.hpp"
  
  #include "geometry/frame.hpp"
  #include "gmock/gmock.h"
  #include "gtest/gtest.h"
  #include "quantities/named_quantities.hpp"
  #include "testing_utilities/almost_equals.hpp"
  #include "testing_utilities/componentwise.hpp"
  
  namespace principia {
  namespace physics {
  namespace internal_dynamic_frame {
  
  using astronomy::InfiniteFuture;
  using astronomy::InfinitePast;
  using geometry::AngularVelocity;
  using geometry::Displacement;
  using geometry::Frame;
  using geometry::Inertial;
  using geometry::InnerProduct;
  using geometry::OrthogonalMap;
  using geometry::Position;
  using geometry::Velocity;
  using quantities::GravitationalParameter;
  using quantities::SIUnit;
  using quantities::Sqrt;
  using quantities::si::Metre;
  using quantities::si::Second;
  using testing_utilities::AlmostEquals;
  using testing_utilities::Componentwise;
  
  namespace {
  
  using Circular = Frame<enum class CircularTag, Inertial>;
  using Helical = Frame<enum class HelicalTag, Inertial>;
  
  Vector<Acceleration, Circular> Gravity(Instant const& t,
                                         Position<Circular> const& q) {
    Displacement<Circular> const r = q - Circular::origin;
    auto const r² = r.Norm²();
    return -SIUnit<GravitationalParameter>() * r / (Sqrt(r²) * r²);
  }
  
  // An inertial frame.
  template<typename OtherFrame, typename ThisFrame>
  class InertialFrame : public DynamicFrame<OtherFrame, ThisFrame> {
   public:
    InertialFrame(
        DegreesOfFreedom<OtherFrame> const& origin_degrees_of_freedom_at_epoch,
        Instant const& epoch,
        OrthogonalMap<OtherFrame, ThisFrame> const& orthogonal_map,
        std::function<Vector<Acceleration, OtherFrame>(
            Instant const& t,
            Position<OtherFrame> const& q)> gravity);
  
    Instant t_min() const override {
      return InfinitePast;
    }
  
    Instant t_max() const override {
      return InfiniteFuture;
    }
  
    RigidMotion<OtherFrame, ThisFrame> ToThisFrameAtTime(
        Instant const& t) const override;
  
    void WriteToMessage(
        not_null<serialization::DynamicFrame*> message) const override;
  
   private:
    Vector<Acceleration, OtherFrame> GravitationalAcceleration(
        Instant const& t,
        Position<OtherFrame> const& q) const override;
    AcceleratedRigidMotion<OtherFrame, ThisFrame> MotionOfThisFrame(
        Instant const& t) const override;
  
    DegreesOfFreedom<OtherFrame> const origin_degrees_of_freedom_at_epoch_;
    Instant const epoch_;
    OrthogonalMap<OtherFrame, ThisFrame> const orthogonal_map_;
    std::function<Vector<Acceleration, OtherFrame>(
        Instant const& t,
        Position<OtherFrame> const& q)> gravity_;
  };
  
  template<typename OtherFrame, typename ThisFrame>
  InertialFrame<OtherFrame, ThisFrame>::InertialFrame(
      DegreesOfFreedom<OtherFrame> const& origin_degrees_of_freedom_at_epoch,
      Instant const& epoch,
      OrthogonalMap<OtherFrame, ThisFrame> const& orthogonal_map,
      std::function<Vector<Acceleration, OtherFrame>(
          Instant const& t,
          Position<OtherFrame> const& q)> gravity)
      : origin_degrees_of_freedom_at_epoch_(origin_degrees_of_freedom_at_epoch),
        epoch_(epoch),
        orthogonal_map_(orthogonal_map),
        gravity_(std::move(gravity)) {}
  
  template<typename OtherFrame, typename ThisFrame>
  RigidMotion<OtherFrame, ThisFrame>
  InertialFrame<OtherFrame, ThisFrame>::ToThisFrameAtTime(
      Instant const& t) const {
    return RigidMotion<OtherFrame, ThisFrame>(
               RigidTransformation<OtherFrame, ThisFrame>(
                   origin_degrees_of_freedom_at_epoch_.position() +
                       (t - epoch_) *
                           origin_degrees_of_freedom_at_epoch_.velocity(),
                   ThisFrame::origin,
                   orthogonal_map_),
               AngularVelocity<OtherFrame>(),
               origin_degrees_of_freedom_at_epoch_.velocity());
  }
  
  template<typename OtherFrame, typename ThisFrame>
  void InertialFrame<OtherFrame, ThisFrame>::WriteToMessage(
      not_null<serialization::DynamicFrame*> message) const {}
  
  template<typename OtherFrame, typename ThisFrame>
  Vector<Acceleration, OtherFrame>
  InertialFrame<OtherFrame, ThisFrame>::GravitationalAcceleration(
      Instant const& t,
      Position<OtherFrame> const& q) const {
    return gravity_(t, q);
  }
  
  template<typename OtherFrame, typename ThisFrame>
  AcceleratedRigidMotion<OtherFrame, ThisFrame>
  InertialFrame<OtherFrame, ThisFrame>::MotionOfThisFrame(
      Instant const& t) const {
    return AcceleratedRigidMotion<OtherFrame, ThisFrame>(
        ToThisFrameAtTime(t),
        /*angular_acceleration_of_to_frame=*/{},
        /*acceleration_of_to_frame_origin=*/{});
  }
  
  }  // namespace
  
  class DynamicFrameTest : public testing::Test {
   protected:
    DegreesOfFreedom<Circular> circular_degrees_of_freedom_ = {
        Circular::origin +
            Displacement<Circular>({1 * Metre, 0 * Metre, 0 * Metre}),
        Velocity<Circular>(
            {0 * Metre / Second, 1 * Metre / Second, 0 * Metre / Second})};
  
    InertialFrame<Circular, Helical> helix_frame_ =
        InertialFrame<Circular, Helical>(
            {Circular::origin,
             Velocity<Circular>(
                 {0 * Metre / Second, 0 * Metre / Second, 1 * Metre / Second})},
            /*epoch=*/Instant() ,
            OrthogonalMap<Circular, Helical>::Identity(),
            &Gravity);
  };
  
  TEST_F(DynamicFrameTest, Helix) {
    auto helix_frenet_frame = helix_frame_.FrenetFrame(
        Instant(),
        helix_frame_.ToThisFrameAtTime(Instant())(circular_degrees_of_freedom_));
    Vector<double, Helical> tangent =
        helix_frenet_frame(Vector<double, Frenet<Helical>>({1, 0, 0}));
    Vector<double, Helical> normal =
        helix_frenet_frame(Vector<double, Frenet<Helical>>({0, 1, 0}));
    EXPECT_THAT(normal, Componentwise(-1, 0, 0));
    EXPECT_THAT(tangent,
                Componentwise(0,
                              AlmostEquals(Sqrt(0.5), 1),
                              AlmostEquals(-Sqrt(0.5), 1)));
  }
  
  }  // namespace internal_dynamic_frame
  }  // namespace physics
  }  // namespace principia
  