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: e379d152b037
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 6a7e2698470b
Choose a head ref

Commits on Nov 6, 2019

  1. Rotation computation.

    pleroy committed Nov 6, 2019
    Copy the full SHA
    312f4a5 View commit details

Commits on Nov 7, 2019

  1. Copy the full SHA
    b97950d View commit details

Commits on Nov 9, 2019

  1. Flipping m.

    pleroy committed Nov 9, 2019
    Copy the full SHA
    9d444f7 View commit details
  2. Case (i).

    pleroy committed Nov 9, 2019
    Copy the full SHA
    ad414c1 View commit details

Commits on Nov 10, 2019

  1. Copy the full SHA
    c1c0b0b View commit details
  2. Copy the full SHA
    68dffd7 View commit details
  3. Copy the full SHA
    a29a60a View commit details
  4. Copy the full SHA
    8f5eb2f View commit details
  5. One test passing.

    pleroy committed Nov 10, 2019
    Copy the full SHA
    a219f84 View commit details
  6. Some more tests passing.

    pleroy committed Nov 10, 2019
    Copy the full SHA
    d223491 View commit details
  7. Case of a motionless body.

    pleroy committed Nov 10, 2019
    Copy the full SHA
    5af2c82 View commit details
  8. Cleanup.

    pleroy committed Nov 10, 2019
    Copy the full SHA
    b1bb7d8 View commit details
  9. Wrong sign in the offset.

    pleroy committed Nov 10, 2019
    Copy the full SHA
    8dfb397 View commit details
  10. Copy the full SHA
    a33988e View commit details

Commits on Nov 11, 2019

  1. Copy the full SHA
    5482c91 View commit details
  2. All tests passing.

    pleroy committed Nov 11, 2019
    Copy the full SHA
    f292446 View commit details
  3. Comments.

    pleroy committed Nov 11, 2019
    Copy the full SHA
    590c3b7 View commit details
  4. Whatever.

    pleroy committed Nov 11, 2019
    Copy the full SHA
    3b12300 View commit details
  5. Removing the σ's.

    pleroy committed Nov 11, 2019
    Copy the full SHA
    0acd4ab View commit details

Commits on Nov 12, 2019

  1. Tolerances.

    pleroy committed Nov 12, 2019
    Copy the full SHA
    79763a7 View commit details
  2. Copy the full SHA
    74b1814 View commit details
  3. After egg's review.

    pleroy committed Nov 12, 2019
    Copy the full SHA
    6e05753 View commit details

Commits on Nov 13, 2019

  1. Merge pull request #2371 from pleroy/Quaternions

    Reimplement the Arnold solver using quaternions
    pleroy authored Nov 13, 2019
    Copy the full SHA
    6a7e269 View commit details
Showing with 351 additions and 191 deletions.
  1. +39 −21 physics/euler_solver.hpp
  2. +277 −127 physics/euler_solver_body.hpp
  3. +35 −43 physics/euler_solver_test.cpp
60 changes: 39 additions & 21 deletions physics/euler_solver.hpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

#pragma once

#include <optional>

#include "geometry/grassmann.hpp"
#include "geometry/named_quantities.hpp"
#include "geometry/r3_element.hpp"
@@ -61,55 +63,71 @@ class EulerSolver {
struct ℬₜ;
struct ℬʹ;

// A frame which is rotated from PrincipalAxesFrame such that the coordinates
// of m along which we project is positive. Used for all internal
// computations.
struct PreferredPrincipalAxesFrame;

using PreferredAngularMomentumBivector =
Bivector<AngularMomentum, PreferredPrincipalAxesFrame>;

// The formula to use, following [CFSZ07], Section 2.2. They don't have a
// formula for the spherical case.
// formula for the spherical case. Also note our singular case for case (iii)
// which arises when the coordinate along which we don't project is 0.
enum class Formula {
i,
ii,
iii,
Sphere
Sphere,
};

Rotation<PrincipalAxesFrame, ℬₜ> Compute𝒫ₜ(
AngularMomentumBivector const& angular_momentum,
bool& ṁ_is_zero) const;
// For case (iii) we use project either on e₁ (as we do in case (i)) or on
// e₃ (as we do in case (ii)) depending on which of the x and z coordinates of
// m is larger (in absolute value).
enum class Region {
e₁,
e₃,
Motionless,
};

// If m is constant in the principal axes frames, we cannot construct ℬₜ using
// ṁ as specified after the demonstration of proposition 2.2 in [CFSZ07].
// Instead, we use a constant v orthogonal to m. This member is set when
// initializating ℛ_.
bool ṁ_is_zero_ = false;
Rotation<PreferredPrincipalAxesFrame, ℬₜ> Compute𝒫ₜ(
PreferredAngularMomentumBivector const& angular_momentum) const;

// Construction parameters.
R3Element<MomentOfInertia> const moments_of_inertia_;
AngularMomentumBivector const initial_angular_momentum_;
Instant const initial_time_;
Rotation<ℬʹ, InertialFrame> const ℛ_;
AngularMomentum const G_;
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> 𝒮_;

// Amusingly, the formula to use is a constant of motion.
// Importantly, the formula and the region to use are constants of motion.
Formula formula_;
Region region_;

// Only the parameters needed for the selected formula are non-NaN after
// construction.

AngularFrequency λ_ = NaN<AngularFrequency>();

AngularMomentum G_ = NaN<AngularMomentum>();
AngularMomentum B₂₃_ = NaN<AngularMomentum>();
AngularMomentum B₁₃_ = NaN<AngularMomentum>();
AngularMomentum B₃₁_ = NaN<AngularMomentum>();
AngularMomentum B₂₁_ = NaN<AngularMomentum>();

AngularMomentum σB₁₃_ = NaN<AngularMomentum>();
AngularMomentum σB₃₁_ = NaN<AngularMomentum>();
AngularMomentum σʹB₁₃_ = NaN<AngularMomentum>();
AngularMomentum σʺB₃₁_ = NaN<AngularMomentum>();

double n_ = NaN<double>();
double mc_ = NaN<double>();
Angle ν_ = NaN<Angle>();
Angle ψ_Π_offset_ = NaN<Angle>();
double ψ_Π_multiplier_ = NaN<double>();
Angle ψ_offset_ = NaN<Angle>();
double ψ_arctan_multiplier_ = NaN<double>();
MomentOfInertia ψ_cn_multiplier_ = NaN<MomentOfInertia>();
MomentOfInertia ψ_sn_multiplier_ = NaN<MomentOfInertia>();
AngularMomentum ψ_cosh_multiplier_ = NaN<AngularMomentum>();
AngularMomentum ψ_sinh_multiplier_ = NaN<AngularMomentum>();
double ψ_integral_multiplier_ = NaN<double>();
AngularFrequency ψ_t_multiplier_ = NaN<AngularFrequency>();
};

404 changes: 277 additions & 127 deletions physics/euler_solver_body.hpp
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@
#include <algorithm>

#include "geometry/grassmann.hpp"
#include "geometry/quaternion.hpp"
#include "numerics/elliptic_functions.hpp"
#include "numerics/elliptic_integrals.hpp"
#include "quantities/elementary_functions.hpp"
@@ -18,6 +19,7 @@ namespace internal_euler_solver {
using geometry::Commutator;
using geometry::DefinesFrame;
using geometry::Normalize;
using geometry::Quaternion;
using geometry::Vector;
using numerics::EllipticF;
using numerics::EllipticΠ;
@@ -29,6 +31,7 @@ using quantities::ArcTanh;
using quantities::Cosh;
using quantities::Energy;
using quantities::Inverse;
using quantities::IsFinite;
using quantities::Pow;
using quantities::Quotient;
using quantities::Sinh;
@@ -49,30 +52,20 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
AttitudeRotation const& initial_attitude,
Instant const& initial_time)
: moments_of_inertia_(moments_of_inertia),
initial_angular_momentum_(initial_angular_momentum),
initial_time_(initial_time),
ℛ_([this, initial_attitude]() -> Rotation<ℬʹ, InertialFrame> {
auto const 𝒴ₜ₀⁻¹ = Rotation<ℬʹ, ℬₜ>::Identity();
auto const 𝒫ₜ₀⁻¹ = Compute𝒫ₜ(initial_angular_momentum_,
ṁ_is_zero_).Inverse();

// This ℛ follows the assumptions in the third paragraph of section 2.3
// of [CFSZ07], that is, the inertial frame is identified with the
// (initial) principal axes frame.
Rotation<ℬʹ, PrincipalAxesFrame> const ℛ = 𝒫ₜ₀⁻¹ * 𝒴ₜ₀⁻¹;

// The multiplication by initial_attitude makes up for the loss of
// generality due to the assumptions in the third paragraph of section
// 2.3 of [CFSZ07].
return initial_attitude * ℛ;
}()) {
G_(initial_angular_momentum.Norm()),
ℛ_(Rotation<ℬʹ, InertialFrame>::Identity()),
𝒮_(Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity()) {
auto const& I₁ = moments_of_inertia_.x;
auto const& I₂ = moments_of_inertia_.y;
auto const& I₃ = moments_of_inertia_.z;
CHECK_LE(I₁, I₂);
CHECK_LE(I₂, I₃);

auto const& m = initial_angular_momentum.coordinates();
// 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();

// These computations are such that if, say I₁ == I₂, I₂₁ is +0.0 and I₁₂ is
// -0.0.
@@ -91,90 +84,208 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
CHECK_LE(Square<AngularMomentum>(), Δ₁);
CHECK_LE(Δ₃, Square<AngularMomentum>());

// These quantities are NaN in the spherical case, so they be used with care
// before we have checked for this case.
auto const B₂₃² = I * Δ / I₂₃;
// These quantities are NaN in the spherical case, so they must be used with
// care before we have checked for this case.
auto const B₃₁² = I * Δ / I₃₁;
auto const B₂₁² = I₂ * Δ₁ / I₂₁;
B₁₃_ = Sqrt(I₁ * Δ₃ / I₁₃);
B₃₁_ = Sqrt(I₃ * Δ₁ / I₃₁);

auto const G² = initial_angular_momentum_.Norm²();
G_ = Sqrt(G²);
auto const two_T = m.x * m.x / I₁ + m.y * m.y / I₂ + m.z * m.z / I₃;
ψ_t_multiplier_ = two_T / G_;
auto const B₂₃² = I₂ * Δ₃ / I₂₃;
auto const B₁₃² = I₁ * Δ₃ / I₁₃;
B₁₃_ = Sqrt(B₁₃²);
B₃₁_ = Sqrt(B₃₁²);

// Note that [CFSZ07] gives k, but we need mc = 1 - k^2. We write mc in a way
// that reduces cancellations when k is close to 1.
// Determine the formula and region to use.
if (Δ₂ < Square<AngularMomentum>()) {
CHECK_LE(Square<AngularMomentum>(), B₂₃²);
CHECK_LE(Square<AngularMomentum>(), B₂₁²);
B₂₁_ = Sqrt(B₂₁²);
mc_ = std::min(Δ₂ * I₃₁ / (Δ₃ * I₂₁), 1.0);
ν_ = EllipticF(ArcTan(m.y * B₃₁_, m.z * B₂₁_), mc_);
auto const λ₃ = Sqrt(Δ₃ * I₁₂ / (I₁ * I₂ * I₃));
// TODO(phl): These tests on the signs of coordinates should probably handle
// -0.0 correctly.
if (m.x < AngularMomentum()) {
σB₁₃_ = -B₁₃_;
λ_ = λ₃;
} else {
σB₁₃_ = B₁₃_;
λ_ = -λ₃;
}
n_ = std::min(G² / B₂₃², 1.0);
ψ_Π_offset_ = EllipticΠ(JacobiAmplitude(-ν_, mc_), n_, mc_);
ψ_Π_multiplier_ = ṁ_is_zero_ ? 0 : Δ₂ / (λ_ * I₂ * G_);
formula_ = Formula::i;
region_ = Region::e₁;
} else if (Square<AngularMomentum>() < Δ₂) {
CHECK_LE(Square<AngularMomentum>(), B₂₃²);
CHECK_LE(Square<AngularMomentum>(), B₂₁²);
B₂₃_ = Sqrt(B₂₃²);
mc_ = std::min(Δ₂ * I₃₁ / (Δ₁ * I₃₂), 1.0);
ν_ = EllipticF(ArcTan(m.y * B₁₃_, m.x * B₂₃_), mc_);
auto const λ₁ = Sqrt(Δ₁ * I₃₂ / (I₁ * I₂ * I₃));
if (m.z < AngularMomentum()) {
σB₃₁_ = -B₃₁_;
λ_ = λ₁;
} else {
σB₃₁_ = B₃₁_;
λ_ = -λ₁;
}
n_ = std::min(G² / B₂₁², 1.0);
ψ_Π_offset_ = EllipticΠ(JacobiAmplitude(-ν_, mc_), n_, mc_);
ψ_Π_multiplier_ = ṁ_is_zero_ ? 0 : Δ₂ / (λ_ * I₂ * G_);
formula_ = Formula::ii;
region_ = Region::e₃;
} else {
CHECK_EQ(Square<AngularMomentum>(), Δ₂);
if (G_ == AngularMomentum()) {
// No rotation. Might as well handle it as a sphere.
ψ_t_multiplier_ = AngularFrequency();
formula_ = Formula::Sphere;
region_ = Region::Motionless;
} else if (I₃₁ == MomentOfInertia()) {
// The degenerate case of a sphere. It would create NaNs.
// The degenerate case of a sphere. It would create NaNs. Pick the
// region that corresponds to the smallest coordinate.
CHECK_EQ(MomentOfInertia(), I₂₁);
CHECK_EQ(MomentOfInertia(), I₃₂);
formula_ = Formula::Sphere;
if (Abs(m.x) < Abs(m.z)) {
region_ = Region::e₁;
} else {
region_ = Region::e₃;
}
} else {
ν_ = -ArcTanh(m.y / G_);
auto const λ₂ = Sqrt(-Δ₁ * Δ₃ / (I₁ * I₃)) / G_;
λ_ = λ₂;
if (m.x < AngularMomentum()) {
σʹB₁₃_ = -B₁₃_;
λ_ = -λ_;
formula_ = Formula::iii;
// Project along the smallest coordinate of x and z in absolute value.
if (B₁₃_ < B₃₁_) {
region_ = Region::e₁;
} else {
σʹB₁₃_ = B₁₃_;
region_ = Region::e₃;
}
if (m.z < AngularMomentum()) {
σʺB₃₁_ = -B₃₁_;
λ_ = -λ_;
}
}

// Compute the rotation 𝒮_ that adjusts the signs of the coordinates of m in a
// way that ensures that the quaternions are well-conditioned and that the σ's
// disappear.
Bivector<double, PreferredPrincipalAxesFrame> e₁({1, 0, 0});
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 {
σʺB₃₁_ = B₃₁_;
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(e₁, -e₂, -e₃);
}
// Δ₂ shows up in the multiplier, and λ_ is finite in the non-spherical
// case, so things simplify tremendously.
ψ_Π_offset_ = Angle();
ψ_Π_multiplier_ = 0;
formula_ = Formula::iii;
} else {
if (m.z >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, -e₂, e₃);
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, e₂, -e₃);
}
}
} else {
switch (region_) {
case Region::e₁: {
if (m.x >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, e₂, -e₃);
}
break;
}
case Region::e₃: {
if (m.z >= AngularMomentum()) {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
} else {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>(-e₁, e₂, -e₃);
}
break;
}
case Region::Motionless: {
𝒮_ = Rotation<PrincipalAxesFrame,
PreferredPrincipalAxesFrame>::Identity();
break;
}
default:
LOG(FATAL) << "Unexpected region " << static_cast<int>(region_);
}
}

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

// This ℛ follows the assumptions in the third paragraph of section 2.3
// of [CFSZ07], that is, the inertial frame is identified with the
// (initial) principal axes frame.
Rotation<ℬʹ, PrincipalAxesFrame> const ℛ = 𝒮⁻¹ * 𝒫ₜ₀⁻¹ * 𝒴ₜ₀⁻¹;

// The multiplication by initial_attitude makes up for the loss of
// generality due to the assumptions in the third paragraph of section
// 2.3 of [CFSZ07].
return initial_attitude * ℛ;
}();

switch (formula_) {
case Formula::i: {
CHECK_LE(Square<AngularMomentum>(), B₂₁²);
B₂₁_ = Sqrt(B₂₁²);
mc_ = std::min(Δ₂ * I₃₁ / (Δ₃ * I₂₁), 1.0);
ν_ = EllipticF(ArcTan(m.y * B₃₁_, m.z * B₂₁_), mc_);
auto const λ₃ = Sqrt(Δ₃ * I₁₂ / (I₁ * I₂ * I₃));
λ_ = -λ₃;

double sn;
double cn;
double dn;
JacobiSNCNDN(-ν_, mc_, sn, cn, dn);
n_ = I₁ * I₃₂ / (I₃ * I₁₂);
ψ_cn_multiplier_ = Sqrt(I₃ * I₂₁);
ψ_sn_multiplier_ = Sqrt(I₂ * I₃₁);
ψ_arctan_multiplier_ = -B₁₃_ * ψ_cn_multiplier_ /
(ψ_sn_multiplier_ * G_);
ψ_offset_ = EllipticΠ(JacobiAmplitude(-ν_, mc_), n_, mc_) +
ψ_arctan_multiplier_ * ArcTan(ψ_sn_multiplier_ * sn,
ψ_cn_multiplier_ * cn);
ψ_integral_multiplier_ = G_ * I₁₃ / (λ_ * I₁ * I₃);
ψ_t_multiplier_ = G_ / I₁;

break;
}
case Formula::ii: {
CHECK_LE(Square<AngularMomentum>(), B₂₃²);
B₂₃_ = Sqrt(B₂₃²);
mc_ = std::min(Δ₂ * I₃₁ / (Δ₁ * I₃₂), 1.0);
ν_ = EllipticF(ArcTan(m.y * B₁₃_, m.x * B₂₃_), mc_);
auto const λ₁ = Sqrt(Δ₁ * I₃₂ / (I₁ * I₂ * I₃));
λ_ = -λ₁;

double sn;
double cn;
double dn;
JacobiSNCNDN(-ν_, mc_, sn, cn, dn);
n_ = I₃ * I₂₁ / (I₁ * I₂₃);
ψ_cn_multiplier_ = Sqrt(I₁ * I₃₂);
ψ_sn_multiplier_ = Sqrt(I₂ * I₃₁);
ψ_arctan_multiplier_ = -B₃₁_ * ψ_cn_multiplier_ /
(ψ_sn_multiplier_ * G_);
ψ_offset_ = EllipticΠ(JacobiAmplitude(-ν_, mc_), n_, mc_) +
ψ_arctan_multiplier_ * ArcTan(ψ_sn_multiplier_ * sn,
ψ_cn_multiplier_ * cn);
ψ_integral_multiplier_ = G_ * I₃₁ / (λ_ * I₁ * I₃);
ψ_t_multiplier_ = G_ / I₃;

break;
}
case Formula::iii: {
ν_ = -ArcTanh(m.y / G_);
auto const λ₂ = Sqrt(-Δ₁ * Δ₃ / (I₁ * I₃)) / G_;
λ_ = λ₂;

switch (region_) {
case Region::e₁: {
ψ_arctan_multiplier_ = 2 * B₁₃_ / B₃₁_;
ψ_cosh_multiplier_ = B₃₁_;
ψ_sinh_multiplier_ = B₁₃_ - G_;
break;
}
case Region::e₃: {
ψ_arctan_multiplier_ = 2 * B₃₁_ / B₁₃_;
ψ_cosh_multiplier_ = B₁₃_;
ψ_sinh_multiplier_ = B₃₁_ - G_;
break;
}
case Region::Motionless:
default:
LOG(FATAL) << "Unexpected region " << static_cast<int>(region_);
}
ψ_offset_ = ArcTan(ψ_sinh_multiplier_ * Tanh(-0.5 * ν_),
ψ_cosh_multiplier_);
auto const two_T = m.x * m.x / I₁ + m.y * m.y / I₂ + m.z * m.z / I₃;
ψ_t_multiplier_ = two_T / G_;

break;
}
case Formula::Sphere: {
ψ_t_multiplier_ = G_ / I₂;
break;
}
}
}
@@ -184,36 +295,42 @@ typename EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularMomentumBivector
EulerSolver<InertialFrame, PrincipalAxesFrame>::AngularMomentumAt(
Instant const& time) const {
Time const Δt = time - initial_time_;
PreferredAngularMomentumBivector m;
switch (formula_) {
case Formula::i: {
double sn;
double cn;
double dn;
JacobiSNCNDN(λ_ * Δt - ν_, mc_, sn, cn, dn);
return AngularMomentumBivector({σB₁₃_ * dn, -B₂₁_ * sn, B₃₁_ * cn});
m = PreferredAngularMomentumBivector({B₁₃_ * dn, -B₂₁_ * sn, B₃₁_ * cn});
break;
}
case Formula::ii: {
double sn;
double cn;
double dn;
JacobiSNCNDN(λ_ * Δt - ν_, mc_, sn, cn, dn);
return AngularMomentumBivector({B₁₃_ * cn, -B₂₃_ * sn, σB₃₁_ * dn});
m = PreferredAngularMomentumBivector({B₁₃_ * cn, -B₂₃_ * sn, B₃₁_ * dn});
break;
}
case Formula::iii: {
Angle const angle = λ_ * Δt - ν_;
double const sech = 1.0 / Cosh(angle);
return AngularMomentumBivector(
{σʹB₁₃_ * sech, G_ * Tanh(angle), σʺB₃₁_ * sech});
m = PreferredAngularMomentumBivector(
{B₁₃_ * sech, G_ * Tanh(angle), B₃₁_ * sech});
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.
return initial_angular_momentum_;
m = initial_angular_momentum_;
break;
}
default:
LOG(FATAL) << "Unexpected formula " << static_cast<int>(formula_);
};
return 𝒮_.Inverse()(m);
}

template<typename InertialFrame, typename PrincipalAxesFrame>
@@ -237,25 +354,43 @@ typename EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeRotation
EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeAt(
AngularMomentumBivector const& angular_momentum,
Instant const& time) const {
bool ṁ_is_zero;
Rotation<PrincipalAxesFrame, ℬₜ> const 𝒫ₜ = Compute𝒫ₜ(angular_momentum,
ṁ_is_zero);
CHECK_EQ(ṁ_is_zero_, ṁ_is_zero);
Rotation<PreferredPrincipalAxesFrame, ℬₜ> const 𝒫ₜ =
Compute𝒫ₜ(𝒮_(angular_momentum));

Time const Δt = time - initial_time_;
Angle ψ = ψ_t_multiplier_ * Δt;
switch (formula_) {
case Formula::i: {
double sn;
double cn;
double dn;
JacobiSNCNDN(λ_ * Δt - ν_, mc_, sn, cn, dn);
Angle const φ = JacobiAmplitude(λ_ * Δt - ν_, mc_);
ψ += ψ_Π_multiplier_ * (EllipticΠ(φ, n_, mc_) - ψ_Π_offset_);
ψ += ψ_integral_multiplier_ *
(EllipticΠ(φ, n_, mc_) +
ψ_arctan_multiplier_ * ArcTan(ψ_sn_multiplier_ * sn,
ψ_cn_multiplier_ * cn) -
ψ_offset_);
break;
}
case Formula::ii: {
double sn;
double cn;
double dn;
JacobiSNCNDN(λ_ * Δt - ν_, mc_, sn, cn, dn);
Angle const φ = JacobiAmplitude(λ_ * Δt - ν_, mc_);
ψ += ψ_Π_multiplier_ * (EllipticΠ(φ, n_, mc_) - ψ_Π_offset_);
ψ += ψ_integral_multiplier_ *
(EllipticΠ(φ, n_, mc_) +
ψ_arctan_multiplier_ * ArcTan(ψ_sn_multiplier_ * sn,
ψ_cn_multiplier_ * cn) -
ψ_offset_);
break;
}
case Formula::iii: {
ψ += ψ_arctan_multiplier_ *
(ArcTan(ψ_sinh_multiplier_ * Tanh(0.5 * (λ_ * Δt - ν_)),
ψ_cosh_multiplier_) -
ψ_offset_);
break;
}
case Formula::Sphere: {
@@ -264,51 +399,66 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeAt(
default:
LOG(FATAL) << "Unexpected formula " << static_cast<int>(formula_);
};
Bivector<double, ℬʹ> const e₃({0, 0, 1});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, e₃, DefinesFrame<ℬₜ>{});

return ℛ_ * 𝒴ₜ * 𝒫ₜ;
switch (region_) {
case Region::e₁: {
Bivector<double, ℬʹ> const e₁({1, 0, 0});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, e₁, DefinesFrame<ℬₜ>{});
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_;
}
case Region::e₃: {
Bivector<double, ℬʹ> const e₃({0, 0, 1});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, e₃, DefinesFrame<ℬₜ>{});
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_;
}
case Region::Motionless: {
Bivector<double, ℬʹ> const unused({0, 1, 0});
Rotation<ℬₜ, ℬʹ> const 𝒴ₜ(ψ, unused, DefinesFrame<ℬₜ>{});
return ℛ_ * 𝒴ₜ * 𝒫ₜ * 𝒮_;
}
default:
LOG(FATAL) << "Unexpected region " << static_cast<int>(region_);
}
}

template<typename InertialFrame, typename PrincipalAxesFrame>
Rotation<PrincipalAxesFrame,
Rotation<typename EulerSolver<InertialFrame,
PrincipalAxesFrame>::PreferredPrincipalAxesFrame,
typename EulerSolver<InertialFrame, PrincipalAxesFrame>::ℬₜ>
EulerSolver<InertialFrame, PrincipalAxesFrame>::Compute𝒫ₜ(
AngularMomentumBivector const& angular_momentum,
bool& ṁ_is_zero) const {
PreferredAngularMomentumBivector const& angular_momentum) const {
auto const& m = angular_momentum;
auto const& m_coordinates = m.coordinates();
auto m_coordinates = m.coordinates();

// Compute ṁ using the Euler equation.
AngularVelocity<PrincipalAxesFrame> const ω = AngularVelocityFor(m);
Bivector<Variation<AngularMomentum>, PrincipalAxesFrame> const ṁ =
Commutator(m, ω) / Radian;

// Construct the orthonormal frame ℬₜ. If m is constant in the principal axes
// frame, the choice is arbitrary.
static Bivector<double, PrincipalAxesFrame> const zero;
ṁ_is_zero = false;
auto const m_normalized = NormalizeOrZero(m);
if (m_normalized == zero) {
ṁ_is_zero = true;
return Rotation<PrincipalAxesFrame, ℬₜ>::Identity();
}
auto v = NormalizeOrZero(ṁ);
if (v == zero) {
ṁ_is_zero = true;
v = NormalizeOrZero(AngularMomentumBivector(
{m_coordinates.y, -m_coordinates.x, AngularMomentum()}));
}
if (v == zero) {
ṁ_is_zero = true;
v = NormalizeOrZero(AngularMomentumBivector(
{AngularMomentum(), -m_coordinates.z, m_coordinates.y}));
Quaternion pₜ;
switch (region_) {
case Region::e₁: {
double const real_part = Sqrt(0.5 * (1 + m_coordinates.x / G_));
AngularMomentum const denominator = 2 * G_ * real_part;
pₜ = Quaternion(real_part,
{0,
m_coordinates.z / denominator,
-m_coordinates.y / denominator});
break;
}
case Region::e₃: {
double const real_part = Sqrt(0.5 * (1 + m_coordinates.z / G_));
AngularMomentum const denominator = 2 * G_ * real_part;
pₜ = Quaternion(real_part,
{m_coordinates.y / denominator,
-m_coordinates.x / denominator,
0});
break;
}
case Region::Motionless: {
pₜ = Quaternion(1);
break;
}
default:
LOG(FATAL) << "Unexpected region " << static_cast<int>(region_);
}
DCHECK_NE(v, zero);
auto const w = Commutator(m_normalized, v);

// 𝒫ₜ(m_normalized).coordinates() = {0, 0, 1} , etc.
Rotation<PrincipalAxesFrame, ℬₜ> const 𝒫ₜ(v, w, m_normalized);
Rotation<PreferredPrincipalAxesFrame, ℬₜ> const 𝒫ₜ(pₜ);

return 𝒫ₜ;
}
78 changes: 35 additions & 43 deletions physics/euler_solver_test.cpp
Original file line number Diff line number Diff line change
@@ -172,7 +172,7 @@ TEST_F(EulerSolverTest, InitialStateRandom) {
solver.AngularMomentumAt(Instant());

EXPECT_THAT(computed_initial_angular_momentum,
AlmostEquals(initial_angular_momentum, 0, 209))
AlmostEquals(initial_angular_momentum, 0, 2136))
<< moments_of_inertia << " " << initial_angular_momentum;
}
}
@@ -383,7 +383,7 @@ TEST_F(EulerSolverTest, ShortFatSymmetricTopPrecession) {
Componentwise(VanishesBefore(1 * SIUnit<AngularMomentum>(), 0, 32),
AlmostEquals(reference_momentum.coordinates().y, 0, 4),
AlmostEquals(reference_momentum.coordinates().z, 0, 2)));
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/13);
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/8);
}

TEST_F(EulerSolverTest, TallSkinnySymmetricTopPrecession) {
@@ -435,7 +435,7 @@ TEST_F(EulerSolverTest, TallSkinnySymmetricTopPrecession) {
Componentwise(AlmostEquals(reference_momentum.coordinates().x, 0, 3),
VanishesBefore(1 * SIUnit<AngularMomentum>(), 0, 24),
AlmostEquals(reference_momentum.coordinates().z, 0, 6)));
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/4);
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/3);
}

// This test demonstrates the Джанибеков effect, also known as tennis racket
@@ -536,9 +536,9 @@ TEST_F(EulerSolverTest, ДжанибековEffect) {
angular_momenta,
attitudes,
Componentwise(VanishesBefore(1 * SIUnit<AngularMomentum>(), 0, 10),
AlmostEquals(reference_momentum.coordinates().y, 0, 13),
AlmostEquals(reference_momentum.coordinates().y, 0, 16),
AlmostEquals(reference_momentum.coordinates().z, 0, 965)));
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/33);
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/45);
}

// A general body that doesn't rotate.
@@ -604,8 +604,8 @@ TEST_F(EulerSolverTest, GeneralBodyRotationAlongFirstAxis) {
actual_angular_velocity,
AlmostEquals(solver.AngularVelocityFor(initial_angular_momentum), 0));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 0));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 1));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 1));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 17));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 17));
}

// A general body that rotates along the second principal axis.
@@ -646,7 +646,7 @@ TEST_F(EulerSolverTest, GeneralBodyRotationAlongSecondAxis) {
EXPECT_THAT(
actual_attitude(e1_),
Componentwise(AlmostEquals(expected_attitude(e1_).coordinates().x, 0),
VanishesBefore(1, 1),
VanishesBefore(1, 0),
AlmostEquals(expected_attitude(e1_).coordinates().z, 1)));
EXPECT_THAT(
actual_attitude(e2_),
@@ -656,7 +656,7 @@ TEST_F(EulerSolverTest, GeneralBodyRotationAlongSecondAxis) {
EXPECT_THAT(
actual_attitude(e3_),
Componentwise(AlmostEquals(expected_attitude(e3_).coordinates().x, 1),
VanishesBefore(1, 0),
VanishesBefore(1, 2),
AlmostEquals(expected_attitude(e3_).coordinates().z, 0)));
}

@@ -695,8 +695,8 @@ TEST_F(EulerSolverTest, GeneralBodyRotationAlongThirdAxis) {
EXPECT_THAT(
actual_angular_velocity,
AlmostEquals(solver.AngularVelocityFor(initial_angular_momentum), 0));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 52));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 52));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 136));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 136));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 0));
}

@@ -738,13 +738,13 @@ TEST_F(EulerSolverTest, GeneralBodyRotationCloseToThirdAxis) {

EXPECT_THAT(
actual_attitude(e1_),
Componentwise(AlmostEquals(expected_attitude(e1_).coordinates().x, 12),
AlmostEquals(expected_attitude(e1_).coordinates().y, 1),
VanishesBefore(std::numeric_limits<double>::epsilon(), 8)));
Componentwise(AlmostEquals(expected_attitude(e1_).coordinates().x, 8),
AlmostEquals(expected_attitude(e1_).coordinates().y, 0, 2),
VanishesBefore(std::numeric_limits<double>::epsilon(), 4)));
EXPECT_THAT(
actual_attitude(e2_),
Componentwise(AlmostEquals(expected_attitude(e2_).coordinates().x, 1),
AlmostEquals(expected_attitude(e2_).coordinates().y, 12),
Componentwise(AlmostEquals(expected_attitude(e2_).coordinates().x, 0, 2),
AlmostEquals(expected_attitude(e2_).coordinates().y, 8),
VanishesBefore(1, 2)));
EXPECT_THAT(
actual_attitude(e3_),
@@ -779,11 +779,8 @@ TEST_F(EulerSolverTest, GeneralBodyRotationVeryCloseToThirdAxis) {
auto const actual_angular_momentum = solver.AngularMomentumAt(t);
auto const actual_angular_velocity =
solver.AngularVelocityFor(actual_angular_momentum);
EXPECT_DEATH({
auto const actual_attitude = solver.AttitudeAt(actual_angular_momentum, t);
}, ".*_is_zero_ == .*_is_zero");
auto const actual_attitude = solver.AttitudeAt(actual_angular_momentum, t);

#if 0
// The expected attitude ignores any precession and is just a rotation
// around z.
auto const expected_angular_frequency = actual_angular_velocity.Norm();
@@ -797,19 +794,18 @@ TEST_F(EulerSolverTest, GeneralBodyRotationVeryCloseToThirdAxis) {
EXPECT_THAT(
actual_attitude(e1_),
Componentwise(AlmostEquals(expected_attitude(e1_).coordinates().x, 0, 8),
AlmostEquals(expected_attitude(e1_).coordinates().y, 0, 1),
Le(std::numeric_limits<double>::min())));
AlmostEquals(expected_attitude(e1_).coordinates().y, 0, 2),
VanishesBefore(1, 0)));
EXPECT_THAT(
actual_attitude(e2_),
Componentwise(AlmostEquals(expected_attitude(e2_).coordinates().x, 0, 1),
Componentwise(AlmostEquals(expected_attitude(e2_).coordinates().x, 0, 2),
AlmostEquals(expected_attitude(e2_).coordinates().y, 0, 8),
Le(std::numeric_limits<double>::min())));
EXPECT_THAT(
actual_attitude(e3_),
Componentwise(Le(std::numeric_limits<double>::min()),
Le(std::numeric_limits<double>::min()),
AlmostEquals(expected_attitude(e3_).coordinates().z, 0)));
#endif
VanishesBefore(1, 0)));
EXPECT_THAT(actual_attitude(e3_),
Componentwise(
VanishesBefore(1, 0),
VanishesBefore(1, 0),
AlmostEquals(expected_attitude(e3_).coordinates().z, 0, 0)));
}

// A sphere that rotates around an axis with a random orientation. Also, a
@@ -855,9 +851,9 @@ TEST_F(EulerSolverTest, SphereRotationAlongRandomAxis) {
AlmostEquals(initial_angular_momentum, 0));
EXPECT_THAT(actual_angular_velocity,
AlmostEquals(expected_angular_velocity, 0));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 29));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 68));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 6));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 8));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 12));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 12));
}

// Rotation on the separatrix with a constant momentum.
@@ -880,11 +876,8 @@ TEST_F(EulerSolverTest, SeparatrixConstantMomentum) {
auto const actual_angular_momentum = solver.AngularMomentumAt(t);
auto const actual_angular_velocity =
solver.AngularVelocityFor(actual_angular_momentum);
EXPECT_DEATH({
auto const actual_attitude = solver.AttitudeAt(actual_angular_momentum, t);
}, ".*_is_zero_ == .*_is_zero");
auto const actual_attitude = solver.AttitudeAt(actual_angular_momentum, t);

#if 0
auto const expected_angular_velocity =
AngularVelocity<PrincipalAxes>({0.0 * Radian / Second,
4.0 / 3.0 * Radian / Second,
@@ -896,13 +889,12 @@ TEST_F(EulerSolverTest, SeparatrixConstantMomentum) {
expected_angular_frequency * 10 * Second, initial_angular_momentum);

EXPECT_THAT(actual_angular_momentum,
AlmostEquals(initial_angular_momentum, 0));
AlmostEquals(initial_angular_momentum, 1));
EXPECT_THAT(actual_angular_velocity,
AlmostEquals(expected_angular_velocity, 0));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 0));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 0));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 0));
#endif
AlmostEquals(expected_angular_velocity, 2));
EXPECT_THAT(actual_attitude(e1_), AlmostEquals(expected_attitude(e1_), 8));
EXPECT_THAT(actual_attitude(e2_), AlmostEquals(expected_attitude(e2_), 4));
EXPECT_THAT(actual_attitude(e3_), AlmostEquals(expected_attitude(e3_), 7));
}

// The data in this test are from Takahashi, Busch and Scheeres, Spin state and