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

Commits on May 30, 2020

  1. Copy the full SHA
    1e8ee98 View commit details
  2. Merge pull request #2597 from pleroy/Diagonalize

    Fix a bug where Diagonalize would produce non-unit quaternions
    pleroy authored May 30, 2020
    Copy the full SHA
    ab70029 View commit details
Showing with 52 additions and 8 deletions.
  1. +4 −2 geometry/symmetric_bilinear_form_body.hpp
  2. +47 −1 geometry/symmetric_bilinear_form_test.cpp
  3. +1 −5 ksp_plugin/pile_up.cpp
6 changes: 4 additions & 2 deletions geometry/symmetric_bilinear_form_body.hpp
Original file line number Diff line number Diff line change
@@ -172,14 +172,16 @@ typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
// Use the Cayley-Hamilton theorem to efficiently find the eigenvectors. The
// m matrices contain, in columns, eigenvectors for the corresponding α.
// However it's possible for a column to be identically 0. To deal with this
// the call to PickEigenvector extracts the column with the largest norm.
// the call to PickEigenvector extracts the column with the largest norm, and
// we make sure that the eigenvectors remain orthonormal.
R3x3Matrix<Scalar> const A_minus_α₀I = A - α₀ * I;
R3x3Matrix<Scalar> const A_minus_α₁I = A - α₁ * I;
R3x3Matrix<Scalar> const A_minus_α₂I = A - α₂ * I;
auto const m₀ = A_minus_α₁I * A_minus_α₂I;
auto const m₁ = A_minus_α₂I * A_minus_α₀I;
auto const v₀ = Vector<double, Frame>(PickEigenvector(m₀));
auto const v₁ = Vector<double, Frame>(PickEigenvector(m₁));
auto const v₁ = Vector<double, Frame>(PickEigenvector(m₁)).
OrthogonalizationAgainst(v₀);

Rotation<Eigenframe, Frame> const rotation{v₀, v₁, Wedge(v₀, v₁)};
return {form, rotation};
48 changes: 47 additions & 1 deletion geometry/symmetric_bilinear_form_test.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

#include "geometry/symmetric_bilinear_form.hpp"

#include <random>

#include "geometry/frame.hpp"
#include "geometry/grassmann.hpp"
#include "geometry/r3_element.hpp"
@@ -347,7 +349,7 @@ TEST_F(SymmetricBilinearFormTest, Diagonalize) {
EXPECT_THAT(f_eigensystem.rotation.Inverse()(w₀),
Componentwise(AlmostEquals(1, 0),
VanishesBefore(1, 0),
VanishesBefore(1, 0)));
VanishesBefore(1, 1)));
EXPECT_THAT(f_eigensystem.rotation.Inverse()(w₁),
Componentwise(VanishesBefore(1, 0),
AlmostEquals(1, 0),
@@ -379,6 +381,50 @@ TEST_F(SymmetricBilinearFormTest, Diagonalize) {
EXPECT_THAT(f_eigensystem.form(e₂, e₁), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₂), AlmostEquals(2 * Metre, 0));
}

// Random matrices.
#if !defined(_DEBUG)
{
std::mt19937_64 random(42);
std::uniform_real_distribution<> inertia_tensor_distribution(-1000.0,
1000.0);
for (int i = 0; i < 1'000'000; ++i) {
double const i00 = inertia_tensor_distribution(random);
double const i01 = inertia_tensor_distribution(random);
double const i02 = inertia_tensor_distribution(random);
double const i11 = inertia_tensor_distribution(random);
double const i12 = inertia_tensor_distribution(random);
double const i22 = inertia_tensor_distribution(random);
auto const f = MakeSymmetricBilinearForm<World>(
R3x3Matrix<double>({i00, i01, i02},
{i01, i11, i12},
{i02, i12, i22}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

EXPECT_THAT(f_eigensystem.rotation.quaternion().Norm(),
AlmostEquals(1.0, 0, 3)) << f;
}
}
#endif

// A poorly-conditioned case that used to yield a non-unit quaternion because
// of non-orthogonal eigenvectors.
{
auto const f = MakeSymmetricBilinearForm<World>(
R3x3Matrix<double>({{+2.86210785240963560e+02,
+4.96371874564241580e+02,
+5.70775558185911450e+02},
{+4.96371874564241580e+02,
+7.16992846801826317e+02,
+8.67180777241129135e+02},
{+5.70775558185911450e+02,
+8.67180777241129135e+02,
+9.56737022360950050e+02}}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

EXPECT_THAT(f_eigensystem.rotation.quaternion().Norm(),
AlmostEquals(1.0, 0)) << f;
}
}

} // namespace internal_symmetric_bilinear_form
6 changes: 1 addition & 5 deletions ksp_plugin/pile_up.cpp
Original file line number Diff line number Diff line change
@@ -430,12 +430,8 @@ void PileUp::DeformPileUpIfNeeded(Instant const& t) {
}
auto const apparent_angular_momentum = apparent_system.AngularMomentum();
auto const apparent_inertia_tensor = apparent_system.InertiaTensor();
auto apparent_inertia_eigensystem =
auto const apparent_inertia_eigensystem =
apparent_inertia_tensor.Diagonalize<PileUpPrincipalAxes>();
// TODO(phl): |Diagonalize| should produce unit quaternions.
apparent_inertia_eigensystem.rotation =
Rotation<PileUpPrincipalAxes, ApparentPileUp>(
Normalize(apparent_inertia_eigensystem.rotation.quaternion()));

Rotation<PileUpPrincipalAxes, ApparentPileUp> const apparent_attitude =
apparent_inertia_eigensystem.rotation;