Five conversions panicked on a condition the caller may have no control over: elements from a non-inertial state, geodetic coordinates from a non-body-fixed one, a look direction across frames, and equinoctial elements for a retrograde orbit. Inclination and frame are data, so a catalog with a retrograde object or a state read from a file would abort the process. ConversionError covers the three cases, and reads as a failure of the value rather than of the code. FrameError stays for a rotation that could not be produced, which is a different thing. The four should_panic tests became error assertions, which say more than they did: the retrograde one now also checks that Keplerian elements still work, since that is the fallback the error leaves open. The debug_asserts in state_derivative are untouched. ordinary-diffeq's right-hand side returns SVector and not Result, so a derivative cannot report anything; moving those checks to assembly time is filed separately. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01M9aqCoSmjgoaTbewDHAZhh
74 lines
2.1 KiB
Rust
74 lines
2.1 KiB
Rust
//! Two-body propagation over one orbital period.
|
|
//!
|
|
//! Run with `cargo run --example two_body`.
|
|
|
|
use astrokit::prelude::*;
|
|
use nalgebra::Vector3;
|
|
|
|
fn main() -> Result<(), Box<dyn std::error::Error>> {
|
|
let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
|
|
let frame = Frame::EARTH_J2000;
|
|
|
|
let initial = Cartesian::new(
|
|
Vector3::new(
|
|
4.263868426884883e6,
|
|
5.146189057155391e6,
|
|
1.1310208421331816e6,
|
|
),
|
|
Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558),
|
|
frame,
|
|
);
|
|
|
|
let elements = initial.keplerian()?;
|
|
let period = Duration::from_seconds(initial.period()?);
|
|
let forces = vec![TwoBody { body: Body::EARTH }];
|
|
|
|
let mut propagation = Propagation::from_cartesian(epoch, initial, &forces)
|
|
.duration(period)
|
|
.accuracy(Accuracy::High);
|
|
let solution = propagation.run();
|
|
|
|
let end = Orbit
|
|
.state(*solution.states.last().unwrap())
|
|
.cartesian(frame);
|
|
|
|
println!("epoch: {epoch}");
|
|
println!("body: {}", initial.frame.center.name());
|
|
println!("period: {period}");
|
|
println!(
|
|
"elements: a = {:.3} km, e = {:.6}, i = {:.3} deg ({:?})",
|
|
elements.semi_major_axis / 1e3,
|
|
elements.eccentricity,
|
|
elements.inclination.to_degrees(),
|
|
elements.degeneracy
|
|
);
|
|
println!(
|
|
"altitude: {:.3} km",
|
|
initial.spherical_altitude() / 1e3
|
|
);
|
|
println!("steps: {}", solution.states.len() - 1);
|
|
println!(
|
|
"position error: {:.3e} m",
|
|
(end.position - initial.position).norm()
|
|
);
|
|
println!(
|
|
"velocity error: {:.3e} m/s",
|
|
(end.velocity - initial.velocity).norm()
|
|
);
|
|
println!(
|
|
"energy drift: {:.3e} J/kg",
|
|
end.specific_energy() - initial.specific_energy()
|
|
);
|
|
|
|
// Dense output: the state at any time in the arc, not only at a step.
|
|
let half = Orbit
|
|
.state(solution.interpolate(period.to_seconds() / 2.0))
|
|
.cartesian(frame);
|
|
println!(
|
|
"half a period: {:.3} km altitude",
|
|
half.spherical_altitude() / 1e3
|
|
);
|
|
|
|
Ok(())
|
|
}
|