Skip to content

Commit

Permalink
Merge pull request #3347 from pleroy/Bisect
Browse files Browse the repository at this point in the history 
Use Brent's method instead of bisection for computing apsides and nodes
  • Loading branch information
@pleroy
pleroy committed Apr 22, 2022
2 parents 0983a2b + 8ea268d commit 08171b5
Show file tree
Hide file tree
Showing 4 changed files with 25 additions and 26 deletions.
16 changes: 8 additions & 8 deletions astronomy/ksp_resonance_test.cpp
View file
Expand Up @@ -287,13 +287,13 @@ TEST_F(KSPResonanceTest, MSVC_ONLY_TEST(Stock)) {
ephemeris->t_max() - 2 * longest_joolian_period_);
EXPECT_THAT(RelativeError(periods_at_mid_term.at(laythe_),
expected_periods_.at(laythe_)),
IsNear(0.874_⑴));
IsNear(0.121_⑴));
EXPECT_THAT(periods_at_mid_term.at(vall_), Eq(Infinity<Time>));
EXPECT_THAT(periods_at_mid_term.at(tylo_), Eq(Infinity<Time>));
EXPECT_THAT(RelativeError(periods_at_mid_term.at(bop_),
expected_periods_.at(bop_)), IsNear(0.38_⑴));
expected_periods_.at(bop_)), IsNear(0.19_⑴));
EXPECT_THAT(RelativeError(periods_at_mid_term.at(pol_),
expected_periods_.at(pol_)), IsNear(0.30_⑴));
expected_periods_.at(pol_)), IsNear(0.13_⑴));

LogEphemeris(*ephemeris,
ephemeris->t_max() - 5 * longest_joolian_period_,
Expand Down Expand Up @@ -357,15 +357,15 @@ TEST_F(KSPResonanceTest, MSVC_ONLY_TEST(Corrected)) {
ComputePeriods(*ephemeris,
ephemeris->t_max() - 2 * longest_joolian_period_);
EXPECT_THAT(RelativeError(periods_at_long_term.at(laythe_),
expected_periods_.at(laythe_)), IsNear(4.7e-3_⑴));
expected_periods_.at(laythe_)), IsNear(4.3e-3_⑴));
EXPECT_THAT(RelativeError(periods_at_long_term.at(vall_),
expected_periods_.at(vall_)), IsNear(5.7e-3_⑴));
expected_periods_.at(vall_)), IsNear(2.9e-3_⑴));
EXPECT_THAT(RelativeError(periods_at_long_term.at(tylo_),
expected_periods_.at(tylo_)), IsNear(0.75e-3_⑴));
expected_periods_.at(tylo_)), IsNear(0.87e-3_⑴));
EXPECT_THAT(RelativeError(periods_at_long_term.at(bop_),
expected_periods_.at(bop_)), IsNear(43e-3_⑴));
expected_periods_.at(bop_)), IsNear(1.3e-3_⑴));
EXPECT_THAT(RelativeError(periods_at_long_term.at(pol_),
expected_periods_.at(pol_)), IsNear(9.9e-3_⑴));
expected_periods_.at(pol_)), IsNear(12e-3_⑴));

LogEphemeris(*ephemeris,
ephemeris->t_max() - 5 * longest_joolian_period_,
Expand Down
9 changes: 4 additions & 5 deletions physics/apsides_body.hpp
View file
Expand Up @@ -19,7 +19,7 @@ using geometry::Barycentre;
using geometry::Instant;
using geometry::Position;
using geometry::Sign;
using numerics::Bisect;
using numerics::Brent;
using numerics::Hermite3;
using quantities::IsFinite;
using quantities::Length;
Expand Down Expand Up @@ -176,10 +176,9 @@ absl::Status ComputeNodes(
node_time = Barycentre<Instant, Length>({*previous_time, time},
{z, -*previous_z});
} else {
// The normal case, find the intersection with z = 0 using bisection.
// TODO(egg): Bisection on a polynomial seems daft; we should have
// Newton's method.
node_time = Bisect(
// The normal case, find the intersection with z = 0 using Brent's
// method.
node_time = Brent(
[&z_approximation](Instant const& t) {
return z_approximation.Evaluate(t);
},
Expand Down
12 changes: 6 additions & 6 deletions physics/ephemeris_body.hpp
View file
Expand Up @@ -49,7 +49,7 @@ using integrators::ExplicitSecondOrderOrdinaryDifferentialEquation;
using integrators::IntegrationProblem;
using integrators::Integrator;
using integrators::methods::Fine1987RKNG34;
using numerics::Bisect;
using numerics::Brent;
using numerics::DoublePrecision;
using numerics::Hermite3;
using quantities::Abs;
Expand Down Expand Up @@ -721,11 +721,11 @@ void Ephemeris<Frame>::ComputeApsides(not_null<MassiveBody const*> const body1,
CHECK(previous_time);

// The derivative of |squared_distance| changed sign. Find its zero by
// bisection, this is the time of the apsis. Then compute the apsis and
// append it to one of the output trajectories.
Instant const apsis_time = Bisect(evaluate_square_distance_derivative,
*previous_time,
time);
// Brent's method, this is the time of the apsis. Then compute the apsis
// and append it to one of the output trajectories.
Instant const apsis_time = Brent(evaluate_square_distance_derivative,
*previous_time,
time);
DegreesOfFreedom<Frame> const apsis1_degrees_of_freedom =
body1_trajectory->EvaluateDegreesOfFreedomLocked(apsis_time);
DegreesOfFreedom<Frame> const apsis2_degrees_of_freedom =
Expand Down
14 changes: 7 additions & 7 deletions physics/kepler_orbit_body.hpp
View file
Expand Up @@ -29,7 +29,7 @@ using geometry::Sign;
using geometry::Vector;
using geometry::Velocity;
using geometry::Wedge;
using numerics::Bisect;
using numerics::Brent;
using quantities::Abs;
using quantities::ArcCos;
using quantities::ArcCosh;
Expand Down Expand Up @@ -642,9 +642,9 @@ void KeplerOrbit<Frame>::CompleteAnomalies(KeplerianElements<Frame>& elements) {
(eccentric_anomaly - e * Sin(eccentric_anomaly) * Radian);
};
Angle const eccentric_anomaly = e == 0 ? *mean_anomaly
: Bisect(kepler_equation,
*mean_anomaly - e * Radian,
*mean_anomaly + e * Radian);
: Brent(kepler_equation,
*mean_anomaly - e * Radian,
*mean_anomaly + e * Radian);
true_anomaly = 2 * ArcTan(Sqrt(1 + e) * Sin(eccentric_anomaly / 2),
Sqrt(1 - e) * Cos(eccentric_anomaly / 2));
hyperbolic_mean_anomaly = NaN<Angle>;
Expand All @@ -657,9 +657,9 @@ void KeplerOrbit<Frame>::CompleteAnomalies(KeplerianElements<Frame>& elements) {
hyperbolic_eccentric_anomaly);
};
Angle const hyperbolic_eccentric_anomaly =
Bisect(hyperbolic_kepler_equation,
0 * Radian,
*hyperbolic_mean_anomaly / (e - 1));
Brent(hyperbolic_kepler_equation,
0 * Radian,
*hyperbolic_mean_anomaly / (e - 1));
true_anomaly =
2 * ArcTan(Sqrt(e + 1) * Sinh(hyperbolic_eccentric_anomaly / 2),
Sqrt(e - 1) * Cosh(hyperbolic_eccentric_anomaly / 2));
Expand Down

0 comments on commit 08171b5

Please sign in to comment.