Skip to content

Commit

Permalink
Merge pull request #2377 from pleroy/Fixes
Browse files Browse the repository at this point in the history
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.
2 parents 6a04769 + c57cf2c commit 2dd0116
Showing 3 changed files with 34 additions and 37 deletions.
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 = {

0 comments on commit 2dd0116

Please sign in to comment.