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: e077ecb33d5d
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 1839adf92b62
Choose a head ref
  • 6 commits
  • 3 files changed
  • 2 contributors

Commits on Jul 19, 2019

  1. Start a test of precession.

    pleroy committed Jul 19, 2019
    Copy the full SHA
    08bdc7a View commit details
  2. A test that fails.

    pleroy committed Jul 19, 2019

    Unverified

    This commit is not signed, but one or more authors requires that any commit attributed to them is signed.
    Copy the full SHA
    0dfb70a View commit details
  3. Copy the full SHA
    1893ae2 View commit details

Commits on Jul 25, 2019

  1. Copy the full SHA
    addb124 View commit details

Commits on Jul 26, 2019

  1. Fix incorrect assertion.

    pleroy committed Jul 26, 2019
    Copy the full SHA
    3b1172a View commit details
  2. Merge pull request #2250 from pleroy/Tests

    A test of the precession of a symmetric top
    pleroy authored Jul 26, 2019
    Copy the full SHA
    1839adf View commit details
Showing with 90 additions and 7 deletions.
  1. +11 −7 physics/euler_solver.cpp
  2. +7 −0 physics/euler_solver.hpp
  3. +72 −0 physics/euler_solver_test.cpp
18 changes: 11 additions & 7 deletions physics/euler_solver.cpp
Original file line number Diff line number Diff line change
@@ -61,17 +61,18 @@ EulerSolver::EulerSolver(

B₁₃_ = Sqrt(I₁ * Δ₃ / I₃₁);
B₃₁_ = Sqrt(I₃ * Δ₁ / I₃₁);
λ₃_ = Sqrt(Δ₃ * I₂₁ / (I₁ * I₂ * I₃));

// Note that Celledoni et al. give k, but we need mc = 1 - k^2. We write mc
// in a way that reduces cancellations when k is close to 1.
if (Δ₂ < Square<AngularMomentum>()) {
B₂₁_ = Sqrt(I₂ * Δ₁ / I₂₁);
mc_ = -Δ₂ * I₃₁ / (Δ₃ * I₂₁);
ν_ = EllipticF(ArcTan(m.y / B₂₁_, m.z / B₃₁_), mc_);
λ₃_ = Sqrt(Δ₃ * I₂₁ / (I₁ * I₂ * I₃));
if (m.x < AngularMomentum()) {
λ₃_ = -λ₃_;
B₁₃_ = -B₁₃_;
} else {
λ₃_ = -λ₃_;
}
formula_ = Formula::i;
} else if (Square<AngularMomentum>() < Δ₂) {
@@ -80,26 +81,29 @@ EulerSolver::EulerSolver(
ν_ = EllipticF(ArcTan(m.y / B₂₃_, m.x / B₁₃_), mc_);
λ₁_ = Sqrt(Δ₁ * I₃₂ / (I₁ * I₂ * I₃));
if (m.z < AngularMomentum()) {
λ₁_ = -λ₁_;
B₃₁_ = -B₃₁_;
} else {
λ₁_ = -λ₁_;
}
formula_ = Formula::ii;
} else {
// Δ₂ == Square<AngularMomentum>()
CHECK_EQ(Square<AngularMomentum>(), Δ₂);
if (I₃₁ == MomentOfInertia()) {
// The degenerate case of a sphere. It would create NaNs.
DCHECK_EQ(MomentOfInertia(), I₂₁);
DCHECK_EQ(AngularFrequency(), λ₃_);
DCHECK_EQ(MomentOfInertia(), I₃₂);
formula_ = Formula::Sphere;
} else {
G_ = initial_angular_momentum_.Norm();
ν_ = -ArcTanh(m.y / G_);
// NOTE(phl): The sign adjustments on this path are unclear.
λ₂_ = Sqrt(Δ₁ * Δ₃ / (I₁ * I₃)) / G_;
if (m.x < AngularMomentum()) {
B₁₃_ = -B₁₃_;
λ₂_ = -λ₂_;
}
if (m.z < AngularMomentum()) {
B₃₁_ = -B₃₁_;
λ₂_ = -λ₂_;
}
formula_ = Formula::iii;
}
@@ -125,7 +129,7 @@ EulerSolver::AngularMomentumBivector EulerSolver::AngularMomentumAt(
return AngularMomentumBivector({B₁₃_ * cn, -B₂₃_ * sn, B₃₁_ * dn});
}
case Formula::iii: {
Angle const angle = λ_ * Δt - ν_;
Angle const angle = λ_ * Δt - ν_;
double const sech = 1.0 / Cosh(angle);
return AngularMomentumBivector(
{B₁₃_ * sech, G_ * Tanh(angle), B₃₁_ * sech});
7 changes: 7 additions & 0 deletions physics/euler_solver.hpp
Original file line number Diff line number Diff line change
@@ -25,6 +25,12 @@ using quantities::NaN;
// A solver for Euler's rotation equations. It follows Celledoni, Fassò,
// Säfström and Zanna, 2007, The exact computation of the free rigid body motion
// and its use in splitting method.
// NOTE(phl): There are a number of errors in the formulæ in Proposition 2.1, as
// can be seen by differentiation:
// In case (i) λ should be defined as -σ λ₃.
// In case (ii) λ should be defined as -σ λ₁.
// In case (iii) the first coordinate should include a factor σ (not σʹ) and
// λ should be defined as σ σʹ λ₂ (where λ₂ is the common value of λ₁ and λ₃).
class EulerSolver {
public:
using PrincipalAxesFrame = Frame<serialization::Frame::PhysicsTag,
@@ -67,6 +73,7 @@ class EulerSolver {
AngularMomentum B₂₃_ = NaN<AngularMomentum>();
AngularMomentum G_ = NaN<AngularMomentum>();
AngularFrequency λ₁_ = NaN<AngularFrequency>();
AngularFrequency λ₂_ = NaN<AngularFrequency>();
AngularFrequency λ₃_ = NaN<AngularFrequency>();
double mc_ = NaN<double>();
Angle ν_ = NaN<Angle>();
72 changes: 72 additions & 0 deletions physics/euler_solver_test.cpp
Original file line number Diff line number Diff line change
@@ -11,17 +11,24 @@
#include "gtest/gtest.h"
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"

namespace principia {
namespace physics {

using geometry::Instant;
using geometry::R3Element;
using quantities::AngularFrequency;
using quantities::AngularMomentum;
using quantities::MomentOfInertia;
using quantities::Cos;
using quantities::Sin;
using quantities::SIUnit;
using quantities::Sqrt;
using quantities::Time;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;

class EulerSolverTest : public ::testing::Test {
@@ -201,5 +208,70 @@ TEST_F(EulerSolverTest, InitialStateFormulæ) {
}
}

// This test and the next come from
// http://n.ethz.ch/~stiegerc/HS09/Mechanik/Unterlagen/Lecture19.pdf.
TEST_F(EulerSolverTest, ShortFatSymmetricTopPrecession) {
R3Element<MomentOfInertia> const moments_of_inertia{
3.0 * SIUnit<MomentOfInertia>(),
3.0 * SIUnit<MomentOfInertia>(),
9.0 * SIUnit<MomentOfInertia>()};

EulerSolver::AngularMomentumBivector const initial_angular_momentum(
{0.0 * SIUnit<AngularMomentum>(),
5.0 * SIUnit<AngularMomentum>(),
7.0 * SIUnit<AngularMomentum>()});

// Correspondence with the referential of lecture 19: x = e1, y = e2, z = e3.
AngularFrequency Ω = initial_angular_momentum.coordinates().z *
(moments_of_inertia[0] - moments_of_inertia[2]) /
(moments_of_inertia[0] * moments_of_inertia[2]);

EulerSolver const solver(
moments_of_inertia, initial_angular_momentum, Instant());
for (Time t = 0 * Second; t < 5.0 * Second; t += 0.1 * Second) {
auto const angular_momentum_at_t = solver.AngularMomentumAt(Instant() + t);
EXPECT_THAT(angular_momentum_at_t,
AlmostEquals(EulerSolver::AngularMomentumBivector(
{5.0 * Sin(Ω * t) * SIUnit<AngularMomentum>(),
5.0 * Cos(Ω * t) * SIUnit<AngularMomentum>(),
7.0 * SIUnit<AngularMomentum>()}),
0,
102))
<< t;
}
}

TEST_F(EulerSolverTest, TallSkinnySymmetricTopPrecession) {
R3Element<MomentOfInertia> const moments_of_inertia{
3.0 * SIUnit<MomentOfInertia>(),
9.0 * SIUnit<MomentOfInertia>(),
9.0 * SIUnit<MomentOfInertia>()};

EulerSolver::AngularMomentumBivector const initial_angular_momentum(
{7.0 * SIUnit<AngularMomentum>(),
0.0 * SIUnit<AngularMomentum>(),
5.0 * SIUnit<AngularMomentum>()});

// Correspondence with the referential of lecture 19: x = e3, y = e1, z = e2.
AngularFrequency Ω = initial_angular_momentum.coordinates().x *
(moments_of_inertia[1] - moments_of_inertia[0]) /
(moments_of_inertia[1] * moments_of_inertia[0]);

EulerSolver const solver(
moments_of_inertia, initial_angular_momentum, Instant());
for (Time t = 0 * Second; t < 5.0 * Second; t += 0.1 * Second) {
auto const angular_momentum_at_t = solver.AngularMomentumAt(Instant() + t);
EXPECT_THAT(angular_momentum_at_t,
AlmostEquals(EulerSolver::AngularMomentumBivector(
{7.0 * SIUnit<AngularMomentum>(),
5.0 * Sin(Ω * t) * SIUnit<AngularMomentum>(),
5.0 * Cos(Ω * t) * SIUnit<AngularMomentum>(),
}),
0,
34))
<< t;
}
}

} // namespace physics
} // namespace principia