Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 315852203622
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: e4d629ece056
Choose a head ref
  • 4 commits
  • 3 files changed
  • 2 contributors

Commits on Dec 28, 2019

  1. Copy the full SHA
    17c21dd View commit details
  2. Lint.

    pleroy committed Dec 28, 2019
    Copy the full SHA
    dcbb7fe View commit details
  3. Merge.

    pleroy committed Dec 28, 2019
    Copy the full SHA
    b39e2fc View commit details
  4. Merge pull request #2418 from pleroy/SolverConstruction

    Change the construction of the solver
    pleroy authored Dec 28, 2019
    Copy the full SHA
    e4d629e View commit details
Showing with 87 additions and 79 deletions.
  1. +14 −10 physics/euler_solver.hpp
  2. +11 −9 physics/euler_solver_body.hpp
  3. +62 −60 physics/euler_solver_test.cpp
24 changes: 14 additions & 10 deletions physics/euler_solver.hpp
Original file line number Diff line number Diff line change
@@ -40,28 +40,32 @@ class EulerSolver {
static_assert(InertialFrame::is_inertial);

public:
using AngularMomentumBivector = Bivector<AngularMomentum, PrincipalAxesFrame>;
using AttitudeRotation = Rotation<PrincipalAxesFrame, InertialFrame>;

// Constructs a solver for a body with the given moments_of_inertia in its
// principal axes frame. The moments must be in increasing order. At
// initial_time the angular momentum is initial_angular_momentum and the
// attitude initial_attitude.
EulerSolver(R3Element<MomentOfInertia> const& moments_of_inertia,
AngularMomentumBivector const& initial_angular_momentum,
AttitudeRotation const& initial_attitude,
Instant const& initial_time);
EulerSolver(
R3Element<MomentOfInertia> const& moments_of_inertia,
Bivector<AngularMomentum, InertialFrame> const& initial_angular_momentum,
AttitudeRotation const& initial_attitude,
Instant const& initial_time);

// Computes the angular momentum at the given time.
AngularMomentumBivector AngularMomentumAt(Instant const& time) const;
// Computes the angular momentum at the given time in the principal axes.
// This is mostly useful as input to the following two functions.
Bivector<AngularMomentum, PrincipalAxesFrame> AngularMomentumAt(
Instant const& time) const;

AngularVelocity<PrincipalAxesFrame> AngularVelocityFor(
AngularMomentumBivector const& angular_momentum) const;
Bivector<AngularMomentum, PrincipalAxesFrame> const& angular_momentum)
const;

// Computes the attitude at the given time, using the angular momentum
// computed by the previous function.
AttitudeRotation AttitudeAt(AngularMomentumBivector const& angular_momentum,
Instant const& time) const;
AttitudeRotation AttitudeAt(
Bivector<AngularMomentum, PrincipalAxesFrame> const& angular_momentum,
Instant const& time) const;

private:
using ℬₜ = Frame<enum class ℬₜTag>;
20 changes: 11 additions & 9 deletions physics/euler_solver_body.hpp
Original file line number Diff line number Diff line change
@@ -50,7 +50,7 @@ using quantities::si::Radian;
template<typename InertialFrame, typename PrincipalAxesFrame>
EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
R3Element<MomentOfInertia> const& moments_of_inertia,
AngularMomentumBivector const& initial_angular_momentum,
Bivector<AngularMomentum, InertialFrame> const& initial_angular_momentum,
AttitudeRotation const& initial_attitude,
Instant const& initial_time)
: moments_of_inertia_(moments_of_inertia),
@@ -59,6 +59,10 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
ℛ_(Rotation<ℬʹ, InertialFrame>::Identity()),
𝒮_(Signature<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity()) {
// Do not use initial_angular_momentum after this point.
auto const initial_angular_momentum_in_principal_axes =
initial_attitude.Inverse()(initial_angular_momentum);

auto const& I₁ = moments_of_inertia_.x;
auto const& I₂ = moments_of_inertia_.y;
auto const& I₃ = moments_of_inertia_.z;
@@ -67,7 +71,7 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(

// The usages of this variable prior to the computation of 𝒮_ must not depend
// on the signs of its coordinates since we may flip it.
auto m = initial_angular_momentum.coordinates();
auto m = initial_angular_momentum_in_principal_axes.coordinates();

// These computations are such that if, say I₁ == I₂, I₂₁ is +0.0 and I₁₂ is
// -0.0.
@@ -162,7 +166,7 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
}

// Now that 𝒮_ has been computed we can use it to adjust m and to compute ℛ_.
initial_angular_momentum_ = 𝒮_(initial_angular_momentum);
initial_angular_momentum_ = 𝒮_(initial_angular_momentum_in_principal_axes);
m = initial_angular_momentum_.coordinates();
ℛ_ = [this, initial_attitude]() -> Rotation<ℬʹ, InertialFrame> {
auto const 𝒴ₜ₀⁻¹ = Rotation<ℬʹ, ℬₜ>::Identity();
@@ -267,7 +271,7 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
}

template<typename InertialFrame, typename PrincipalAxesFrame>
typename EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularMomentumBivector
Bivector<AngularMomentum, PrincipalAxesFrame>
EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularMomentumAt(
Instant const& time) const {
Time const Δt = time - initial_time_;
@@ -297,9 +301,6 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularMomentumAt(
break;
}
case Formula::Sphere : {
// NOTE(phl): It's unclear how the formulæ degenerate in this case, but
// surely λ₃_ becomes 0, so the dependency in time disappears, so this is
// my best guess.
m = initial_angular_momentum_;
break;
}
@@ -312,7 +313,8 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularMomentumAt(
template<typename InertialFrame, typename PrincipalAxesFrame>
AngularVelocity<PrincipalAxesFrame>
EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularVelocityFor(
AngularMomentumBivector const& angular_momentum) const {
Bivector<AngularMomentum, PrincipalAxesFrame> const& angular_momentum)
const {
auto const& m = angular_momentum;
auto const& m_coordinates = m.coordinates();

@@ -328,7 +330,7 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularVelocityFor(
template<typename InertialFrame, typename PrincipalAxesFrame>
typename EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeRotation
EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeAt(
AngularMomentumBivector const& angular_momentum,
Bivector<AngularMomentum, PrincipalAxesFrame> const& angular_momentum,
Instant const& time) const {
Rotation<PreferredPrincipalAxesFrame, ℬₜ> const 𝒫ₜ =
Compute𝒫ₜ(𝒮_(angular_momentum));
Loading