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

Commits on Nov 23, 2019

  1. Copy the full SHA
    c57cf2c View commit details
  2. Merge pull request #2377 from pleroy/Fixes

    Separate out the code changes from #2375
    pleroy authored Nov 23, 2019

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    2dd0116 View commit details
Showing with 34 additions and 37 deletions.
  1. +1 −1 physics/euler_solver.hpp
  2. +23 −26 physics/euler_solver_body.hpp
  3. +10 −10 physics/euler_solver_test.cpp
2 changes: 1 addition & 1 deletion physics/euler_solver.hpp
Original file line number Diff line number Diff line change
@@ -127,7 +127,7 @@ class EulerSolver {
MomentOfInertia ψ_sn_multiplier_ = NaN<MomentOfInertia>();
AngularMomentum ψ_cosh_multiplier_ = NaN<AngularMomentum>();
AngularMomentum ψ_sinh_multiplier_ = NaN<AngularMomentum>();
double ψ_integral_multiplier_ = NaN<double>();
double ψ_elliptic_pi_multiplier_ = NaN<double>();
AngularFrequency ψ_t_multiplier_ = NaN<AngularFrequency>();
};

49 changes: 23 additions & 26 deletions physics/euler_solver_body.hpp
Original file line number Diff line number Diff line change
@@ -219,12 +219,12 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
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₃);
ψ_arctan_multiplier_ = -1;
ψ_elliptic_pi_multiplier_ = G_ * I₁₃ / (λ_ * I₁ * I₃);
ψ_offset_ = ψ_elliptic_pi_multiplier_ *
EllipticΠ(JacobiAmplitude(-ν_, mc_), n_, mc_) +
ψ_arctan_multiplier_ *
ArcTan(ψ_sn_multiplier_ * sn, ψ_cn_multiplier_ * cn);
ψ_t_multiplier_ = G_ / I₁;

break;
@@ -244,12 +244,12 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
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₃);
ψ_arctan_multiplier_ = 1;
ψ_elliptic_pi_multiplier_ = G_ * I₃₁ / (λ_ * I₁ * I₃);
ψ_offset_ = ψ_elliptic_pi_multiplier_ *
EllipticΠ(JacobiAmplitude(-ν_, mc_), n_, mc_) +
ψ_arctan_multiplier_ *
ArcTan(ψ_sn_multiplier_ * sn, ψ_cn_multiplier_ * cn);
ψ_t_multiplier_ = G_ / I₃;

break;
@@ -261,13 +261,13 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(

switch (region_) {
case Region::e₁: {
ψ_arctan_multiplier_ = 2 * B₁₃_ / B₃₁_;
ψ_arctan_multiplier_ = -2;
ψ_cosh_multiplier_ = B₃₁_;
ψ_sinh_multiplier_ = B₁₃_ - G_;
break;
}
case Region::e₃: {
ψ_arctan_multiplier_ = 2 * B₃₁_ / B₁₃_;
ψ_arctan_multiplier_ = 2;
ψ_cosh_multiplier_ = B₁₃_;
ψ_sinh_multiplier_ = B₃₁_ - G_;
break;
@@ -278,8 +278,7 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::EulerSolver(
}
ψ_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_;
ψ_t_multiplier_ = G_ / I₂;

break;
}
@@ -366,11 +365,10 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeAt(
double dn;
JacobiSNCNDN(λ_ * Δt - ν_, mc_, sn, cn, dn);
Angle const φ = JacobiAmplitude(λ_ * Δt - ν_, mc_);
ψ += ψ_integral_multiplier_ *
(EllipticΠ(φ, n_, mc_) +
ψ_arctan_multiplier_ * ArcTan(ψ_sn_multiplier_ * sn,
ψ_cn_multiplier_ * cn) -
ψ_offset_);
ψ += ψ_elliptic_pi_multiplier_ * EllipticΠ(φ, n_, mc_) +
ψ_arctan_multiplier_ *
ArcTan(ψ_sn_multiplier_ * sn, ψ_cn_multiplier_ * cn) -
ψ_offset_;
break;
}
case Formula::ii: {
@@ -379,11 +377,10 @@ EulerSolver<InertialFrame, PrincipalAxesFrame>::AttitudeAt(
double dn;
JacobiSNCNDN(λ_ * Δt - ν_, mc_, sn, cn, dn);
Angle const φ = JacobiAmplitude(λ_ * Δt - ν_, mc_);
ψ += ψ_integral_multiplier_ *
(EllipticΠ(φ, n_, mc_) +
ψ_arctan_multiplier_ * ArcTan(ψ_sn_multiplier_ * sn,
ψ_cn_multiplier_ * cn) -
ψ_offset_);
ψ += ψ_elliptic_pi_multiplier_ * EllipticΠ(φ, n_, mc_) +
ψ_arctan_multiplier_ *
ArcTan(ψ_sn_multiplier_ * sn, ψ_cn_multiplier_ * cn) -
ψ_offset_;
break;
}
case Formula::iii: {
20 changes: 10 additions & 10 deletions physics/euler_solver_test.cpp
Original file line number Diff line number Diff line change
@@ -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=*/3);
CheckPoinsotConstruction(solver, angular_momenta, attitudes, /*ulps=*/2);
}

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

// A general body that doesn't rotate.
@@ -898,14 +898,14 @@ TEST_F(EulerSolverTest, SeparatrixConstantMomentum) {
}

// The data in this test are from Takahashi, Busch and Scheeres, Spin state and
// moment of inertia characterization of 4179 Toutatis, 2013.
// moment of inertia characterization of 4179 Toutatis, 2013 [TBS13].
TEST_F(EulerSolverTest, Toutatis) {
Instant const epoch = "1992-11-09T17:49:47"_UTC;

// Takahashi et al. adopt a bizarre convention where their x axis is our I₂,
// their y axis is our I₃ and their z axis is our I₁, see Table 2. This
// appears to contradict Figure 1, but it is consistent with ω₁, ω₂, ω₃ being
// along their x, y, z axes respectively.
// [TBS13] adopt a bizarre convention where their x axis is our I₂, their y
// axis is our I₃ and their z axis is our I₁, see Table 2. This appears to
// contradict Figure 1, but it is consistent with ω₁, ω₂, ω₃ being along their
// x, y, z axes respectively.
struct TakahashiPrincipalAxes;
using TakahashiAttitudeRotation = Rotation<TakahashiPrincipalAxes, ICRS>;
using TakahashiPermutation = Permutation<TakahashiPrincipalAxes,
@@ -939,7 +939,7 @@ TEST_F(EulerSolverTest, Toutatis) {
takahashi_moments_of_inertia.z});

// From Zhao et al., Orientation and rotational parameters of asteroid 4179
// Toutatis: new insights from Change'e-2's close flyby.
// Toutatis: new insights from Change'e-2's close flyby [Zha+15].
auto const angular_momentum_orientation_in_inertial = Bivector<double, ICRS>(
RadiusLatitudeLongitude(1.0, -54.75 * Degree, 180.2 * Degree)
.ToCartesian());
@@ -1002,12 +1002,12 @@ TEST_F(EulerSolverTest, Toutatis) {
ApproximateQuantity<Angle> e3_direction_error;
};

// The errors for 1992 seem consistent with the residuals from Takahashi (3 σ
// The errors for 1992 seem consistent with the residuals from [TBS13] (3 σ
// for the attitude on 1992-12-04, 1.5 σ on 1992-12-07, for instance, figures
// 7 and 8) and they are generally smaller for the angular velocity than for
// the attitude, in agreement with that paper.
// The errors for 2008 are caused by the fact that we ignore the torque
// exerted by the Earth and the Sun. As shown in figure 6 of Takahashi, this
// exerted by the Earth and the Sun. As shown in figure 6 of [TBS13], this
// results in a completely different orientation, estimated in the text to be
// around 100°, which is again generally consistent with our results.
std::vector<Observation> const observations = {