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

Commits on Oct 2, 2019

  1. Copy the full SHA
    f7cb12e View commit details
  2. After egg's review.

    pleroy committed Oct 2, 2019
    Copy the full SHA
    e68f3a4 View commit details
  3. Merge pull request #2349 from pleroy/Conformity2

    Change the code to conform to the documentation in Celledoni.pdf
    pleroy authored Oct 2, 2019

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    a425a94 View commit details
Showing with 10 additions and 10 deletions.
  1. +10 −10 physics/euler_solver_body.hpp
20 changes: 10 additions & 10 deletions physics/euler_solver_body.hpp
Original file line number Diff line number Diff line change
@@ -58,31 +58,31 @@ EulerSolver<PrincipalAxesFrame>::EulerSolver(
// The formulæ for the Δs in Celledoni cannot be used directly because of
// cancellations.
auto const Δ₁ = m.y * m.y * I₂₁ / I₂ + m.z * m.z * I₃₁ / I₃;
auto const Δ₂ = m.x * m.x * I₂ / I + m.z * m.z * I₂ / I;
auto const Δ₃ = m.x * m.x * I₃₁ / I₁ + m.y * m.y * I₃₂ / I₂;
auto const Δ₂ = m.z * m.z * I₂ / I + m.x * m.x * I₂ / I;
auto const Δ₃ = m.x * m.x * I₁₃ / I₁ + m.y * m.y * I₂₃ / I₂;
DCHECK_LE(Square<AngularMomentum>(), Δ₁);
DCHECK_LE(Square<AngularMomentum>(), Δ₃);
DCHECK_LE(Δ₃, Square<AngularMomentum>());

B₁₃_ = Sqrt(I₁ * Δ₃ / I₃₁);
B₁₃_ = Sqrt(I₁ * Δ₃ / I₁₃);
B₃₁_ = Sqrt(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₃));
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>() < Δ₂) {
B₂₃_ = Sqrt(I₂ * Δ₃ / I₃₂);
B₂₃_ = Sqrt(I₂ * Δ₃ / I₂₃);
mc_ = Δ₂ * I₃₁ / (Δ₁ * I₃₂);
ν_ = EllipticF(ArcTan(m.y / B₂₃_, m.x / B₁₃_), mc_);
ν_ = EllipticF(ArcTan(m.y * B₁₃_, m.x * B₂₃_), mc_);
λ₁_ = Sqrt(Δ₁ * I₃₂ / (I₁ * I₂ * I₃));
if (m.z < AngularMomentum()) {
B₃₁_ = -B₃₁_;
@@ -100,7 +100,7 @@ EulerSolver<PrincipalAxesFrame>::EulerSolver(
} else {
G_ = initial_angular_momentum_.Norm();
ν_ = -ArcTanh(m.y / G_);
λ₂_ = Sqrt(Δ₁ * Δ₃ / (I₁ * I₃)) / G_;
λ₂_ = Sqrt(-Δ₁ * Δ₃ / (I₁ * I₃)) / G_;
if (m.x < AngularMomentum()) {
B₁₃_ = -B₁₃_;
λ₂_ = -λ₂_;