Check failed: std::gcd(Dᴛₒ, Cᴛₒ) == 1

[lpgagnon]

Seen twice in a row, but wasn't able to reproduce with journaling on.

    @        0x69884d844  principia::astronomy::_orbit_recurrence::internal::OrbitRecurrence::OrbitRecurrence()
    @        0x698875de1  principia::ksp_plugin::_orbit_analyser::internal::OrbitAnalyser::AnalyseOrbit()
    @        0x698879958  _ZNSt3__1L14__thread_proxyINS_5tupleIJNS_10unique_ptrINS_15__thread_structENS_14default_deleteIS3_EEEEZN9principia4base8_jthread8internalL19MakeStoppableThreadIZNS7_10ksp_plugin15_orbit_analyser8internal13OrbitAnalyser15RequestAnalysisERKNSF_10ParametersEE3$_0JEEENSA_7jthreadEOT_DpOT0_EUlRKNSA_10stop_tokenEE_SQ_EEEEEPvSV_
    @     0x7ff8142ff1d3  _pthread_start
F1231 20:12:02.344707 3476037632 orbit_recurrence_body.hpp:33] Check failed: std::gcd(Dᴛₒ, Cᴛₒ) == 1 (16 vs. 1) [-2147483648; -2147483648; -818929872]

Context: vessel trajectory is Uranus flyby, roughly in Triton's plane, with Pe near Triton altitude. Creating a flight plan for Triton capture: one maneuver plotted at Pe to achieve highly elliptical Uranus capture, and finetuning it to get a Triton encounter on the next pass.

Both times, crash occurred while fiddling to get Triton Pe just above its atmosphere.

check_failed_ 3831.zip

[eggrobin]

Unfortunately we do not have a stacktrace decoder for mac, but I think this is it.

νₒ = Dᴛₒ = -2147483648 might be std::nearbyint on a NaN or infinity. Cᴛₒ = -818929872 is strange. Maybe uninitialized?

Principia/ksp_plugin/orbit_analyser.cpp

Lines 101 to 155 in 4eb331d

	absl::Status OrbitAnalyser::AnalyseOrbit(Parameters const& parameters) {
	Analysis analysis{parameters.first_time};
	RotatingBody<Barycentric> const* primary = nullptr;
	auto smallest_osculating_period = Infinity<Time>;
	auto const primary_status =
	FindBodyWithSmallestOsculatingPeriod(parameters,
	primary,
	smallest_osculating_period);
	RETURN_IF_ERROR(primary_status);
	if (primary != nullptr) {
	DiscreteTrajectory<Barycentric> trajectory;
	trajectory.segments().front().SetDownsampling(
	OrbitAnalyserDownsamplingParameters());
	Time const analysis_duration = std::min(
	parameters.extended_mission_duration.value_or(
	parameters.mission_duration),
	std::max(2 * smallest_osculating_period, parameters.mission_duration));
	RETURN_IF_ERROR(
	FlowWithProgressBar(parameters, analysis_duration, trajectory));
	analysis.mission_duration_ = trajectory.back().time - parameters.first_time;
	// TODO(egg): |next_analysis_percentage_| only reflects the progress of
	// the integration, but the analysis itself can take a while; this results
	// in the progress bar being stuck at 100% while the elements and nodes
	// are being computed.
	BodyCentredNonRotatingReferenceFrame<Barycentric, PrimaryCentred> const
	primary_centred(ephemeris_, primary);
	auto const status_or_primary_centred_trajectory =
	ToPrimaryCentred(primary_centred, trajectory);
	RETURN_IF_ERROR(status_or_primary_centred_trajectory);
	auto const& primary_centred_trajectory =
	status_or_primary_centred_trajectory.value();
	analysis.primary_ = primary;
	analysis.radial_distance_interval_ =
	RadialDistanceInterval(primary_centred_trajectory);
	auto elements = OrbitalElements::ForTrajectory(
	primary_centred_trajectory, *primary, MasslessBody{});
	// We do not RETURN_IF_ERROR as ForTrajectory can return non-CANCELLED
	// statuses.
	RETURN_IF_STOPPED;
	if (elements.ok()) {
	analysis.elements_ = std::move(elements).value();
	// TODO(egg): max_abs_Cᴛₒ should probably depend on the number of
	// revolutions.
	analysis.closest_recurrence_ = OrbitRecurrence::ClosestRecurrence(
	analysis.elements_->nodal_period(),
	analysis.elements_->nodal_precession(),
	*primary,
	/*max_abs_Cᴛₒ=*/100);

Principia/astronomy/orbit_recurrence_body.hpp

Lines 51 to 79 in 4eb331d

	OrbitRecurrence OrbitRecurrence::ClosestRecurrence(
	Time const& nodal_period,
	AngularFrequency const& nodal_precession,
	RotatingBody<Frame> const& primary,
	int max_abs_Cᴛₒ) {
	AngularFrequency const& Ωʹ = nodal_precession;
	AngularFrequency const& Ωʹᴛ = primary.angular_frequency();
	// Nodal mean motion.
	AngularFrequency const& nd = 2 * π * Radian / nodal_period;
	// Daily recurrence frequency, see (7.41).
	double const κ = nd / (Ωʹᴛ - Ωʹ);
	// Look for the closest rational approximation Nᴛₒ / Cᴛₒ to κ whose
	// denominator is at most max_Cᴛₒ.
	// The notation follows section 11.7.2.
	int Cᴛₒ;
	double min_frac_abs_κ_J = std::numeric_limits<double>::infinity();
	for (int J = 1; J <= max_abs_Cᴛₒ; ++J) {
	double const abs_κ_J = Abs(κ * J);
	double const frac_abs_κ_J = Abs(abs_κ_J - std::nearbyint(abs_κ_J));
	if (frac_abs_κ_J < min_frac_abs_κ_J) {
	min_frac_abs_κ_J = frac_abs_κ_J;
	Cᴛₒ = Sign(κ) * J;
	}
	}
	int const νₒ = std::nearbyint(κ);
	int const Dᴛₒ = std::nearbyint((κ - νₒ) * Cᴛₒ);
	return OrbitRecurrence(νₒ, Dᴛₒ, Cᴛₒ);

Principia/astronomy/orbit_recurrence_body.hpp

Lines 23 to 33 in 4eb331d

	inline OrbitRecurrence::OrbitRecurrence(int const νₒ,
	int const Dᴛₒ,
	int const Cᴛₒ)
	: νₒ_(νₒ), Dᴛₒ_(Dᴛₒ), Cᴛₒ_(Cᴛₒ) {
	CHECK_NE(Cᴛₒ, 0) << *this;
	Sign const sign_Cᴛₒ = Sign::OfNonZero(Cᴛₒ);
	if (νₒ != 0) {
	CHECK_EQ(Sign::OfNonZero(νₒ), sign_Cᴛₒ) << *this;
	}
	CHECK_LE(Abs(2 * Dᴛₒ), Abs(Cᴛₒ)) << *this;
	CHECK_EQ(std::gcd(Dᴛₒ, Cᴛₒ), 1) << *this;