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: fde2f2d8ff98
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 315852203622
Choose a head ref
  • 5 commits
  • 2 files changed
  • 1 contributor

Commits on Dec 28, 2019

  1. Copy the full SHA
    7078624 View commit details
  2. Copy the full SHA
    776afe6 View commit details
  3. construction

    eggrobin committed Dec 28, 2019
    Copy the full SHA
    f6ca628 View commit details
  4. after pleroy's review

    eggrobin committed Dec 28, 2019
    Copy the full SHA
    3e128cf View commit details
  5. Merge pull request #2416 from eggrobin/signature

    Use Signature in the Euler-Arnold solver
    eggrobin authored Dec 28, 2019

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    3158522 View commit details
Showing with 21 additions and 42 deletions.
  1. +5 −3 physics/euler_solver.hpp
  2. +16 −39 physics/euler_solver_body.hpp
8 changes: 5 additions & 3 deletions physics/euler_solver.hpp
Original file line number Diff line number Diff line change
@@ -8,6 +8,7 @@
#include "geometry/named_quantities.hpp"
#include "geometry/r3_element.hpp"
#include "geometry/rotation.hpp"
#include "geometry/signature.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/quantities.hpp"

@@ -21,6 +22,7 @@ using geometry::Frame;
using geometry::Instant;
using geometry::R3Element;
using geometry::Rotation;
using geometry::Signature;
using quantities::Angle;
using quantities::AngularFrequency;
using quantities::AngularMomentum;
@@ -103,9 +105,9 @@ class EulerSolver {
PreferredAngularMomentumBivector initial_angular_momentum_;
Rotation<ℬʹ, InertialFrame> ℛ_;

// A rotation that describes which axes are flipped to adjust the signs of the
// coordinates of m. It incorporates σ, σʹ and σʺ from [CFSZ07].
Rotation<PrincipalAxesFrame, PreferredPrincipalAxesFrame> 𝒮_;
// A signature that describes which axes are flipped to adjust the signs of
// the coordinates of m. It incorporates σ, σʹ and σʺ from [CFSZ07].
Signature<PrincipalAxesFrame, PreferredPrincipalAxesFrame> 𝒮_;

// Importantly, the formula and the region to use are constants of motion.
Formula formula_;
55 changes: 16 additions & 39 deletions physics/euler_solver_body.hpp
Original file line number Diff line number Diff line change
@@ -17,9 +17,11 @@ namespace physics {
namespace internal_euler_solver {

using geometry::Commutator;
using geometry::DeduceSignPreservingOrientation;
using geometry::DefinesFrame;
using geometry::Normalize;
using geometry::Quaternion;
using geometry::Sign;
using geometry::Vector;
using numerics::EllipticF;
using numerics::EllipticΠ;
@@ -55,8 +57,8 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
initial_time_(initial_time),
G_(initial_angular_momentum.Norm()),
ℛ_(Rotation<ℬʹ, InertialFrame>::Identity()),
𝒮_(Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity()) {
𝒮_(Signature<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity()) {
auto const& I₁ = moments_of_inertia_.x;
auto const& I₂ = moments_of_inertia_.y;
auto const& I₃ = moments_of_inertia_.z;
@@ -135,48 +137,23 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
Bivector<double, PreferredPrincipalAxesFrame> e₂({0, 1, 0});
Bivector<double, PreferredPrincipalAxesFrame> e₃({0, 0, 1});
if (formula_ == Formula::iii) {
if (m.x >= AngularMomentum()) {
if (m.z >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(e₁, -e₂, -e₃);
}
} else {
if (m.z >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, -e₂, e₃);
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, e₂, -e₃);
}
}
𝒮_ = Signature<PrincipalAxesFrame, PreferredPrincipalAxesFrame>(
Sign(m.x), DeduceSignPreservingOrientation{}, Sign(m.z));
} else {
switch (region_) {
case Region::e₁: {
if (m.x >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, e₂, -e₃);
}
𝒮_ = Signature<PrincipalAxesFrame, PreferredPrincipalAxesFrame>(
Sign(m.x), Sign::Positive(), DeduceSignPreservingOrientation{});
break;
}
case Region::e₃: {
if (m.z >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, e₂, -e₃);
}
𝒮_ = Signature<PrincipalAxesFrame, PreferredPrincipalAxesFrame>(
DeduceSignPreservingOrientation{}, Sign::Positive(), Sign(m.z));
break;
}
case Region::Motionless: {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
𝒮_ = Signature<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
break;
}
default:
@@ -190,7 +167,7 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
ℛ_ = [this, initial_attitude]() -> Rotation<ℬʹ, InertialFrame> {
auto const 𝒴ₜ₀⁻¹ = Rotation<ℬʹ, ℬₜ>::Identity();
auto const 𝒫ₜ₀⁻¹ = Compute𝒫ₜ(initial_angular_momentum_).Inverse();
auto const 𝒮⁻¹ = 𝒮_.Inverse();
auto const 𝒮⁻¹ = 𝒮_.Inverse().Forget().AsRotation();

// This ℛ follows the assumptions in the third paragraph of section 2.3
// of [CFSZ07], that is, the inertial frame is identified with the
@@ -401,17 +378,17 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeAt(
case Region::e₁: {
Bivector<double, ℬʹ> const e₁({1, 0, 0});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, e₁, DefinesFrame<ℬₜ>{});
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_;
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_.Forget().AsRotation();
}
case Region::e₃: {
Bivector<double, ℬʹ> const e₃({0, 0, 1});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, e₃, DefinesFrame<ℬₜ>{});
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_;
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_.Forget().AsRotation();
}
case Region::Motionless: {
Bivector<double, ℬʹ> const unused({0, 1, 0});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, unused, DefinesFrame<ℬₜ>{});
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_;
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_.Forget().AsRotation();
}
default:
LOG(FATAL) << "Unexpected region " << static_cast<int>(region_);