113 lines
4.2 KiB
Rust
113 lines
4.2 KiB
Rust
//! Orbit and attitude propagated together, with the area that follows from them.
|
|
//!
|
|
//! Run with `cargo run --example attitude`.
|
|
|
|
use astrokit::attitude;
|
|
use astrokit::prelude::*;
|
|
use nalgebra::Vector3;
|
|
|
|
fn main() {
|
|
let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
|
|
let frame = Frame::EARTH_J2000;
|
|
let forces = vec![TwoBody { body: Body::EARTH }];
|
|
|
|
// A 400 km circular orbit at the inclination of the space station. Circular, so
|
|
// the argument of periapsis is undefined and the last angle is the argument of
|
|
// latitude; `Keplerian::new` works that out rather than being told.
|
|
let orbit = Keplerian::new(6.7781e6, 0.0, 51.6_f64.to_radians(), 0.4, 0.0, 0.0, frame);
|
|
let spacecraft = Spacecraft::rectangular_box(500.0, Vector3::new(2.0, 1.0, 3.0), 2.2, 1.3);
|
|
// A tumble about no principal axis, so all three rates change.
|
|
let rate = Vector3::new(0.0012, 0.0028, -0.0007);
|
|
|
|
let initial = Full
|
|
.with_keplerian(orbit)
|
|
.with_mass(spacecraft.mass)
|
|
.with_angular_rate(rate);
|
|
|
|
let period = orbit.period();
|
|
let mut propagation = Propagation::new(epoch, initial, &forces)
|
|
.frame(frame)
|
|
.spacecraft(spacecraft.clone())
|
|
.duration(Duration::from_seconds(period))
|
|
.accuracy(Accuracy::High);
|
|
let solution = propagation.run();
|
|
|
|
let state = Full.state(*solution.states.last().unwrap());
|
|
// Switches of the parameters to their shadow set, which the propagation does for
|
|
// itself whenever the layout carries an attitude.
|
|
let switches = jump_times(&solution.times).len();
|
|
|
|
// The same arc without the attitude, to confirm the orbit is untouched by it.
|
|
let mut translational = Propagation::new(epoch, Orbit.with_keplerian(orbit), &forces)
|
|
.frame(frame)
|
|
.duration(Duration::from_seconds(period))
|
|
.accuracy(Accuracy::High);
|
|
let reference = Orbit
|
|
.state(*translational.run().states.last().unwrap())
|
|
.cartesian(frame);
|
|
|
|
let inertia = &spacecraft.inertia;
|
|
let momentum = |w: Vector3<f64>| (inertia * w).norm();
|
|
let energy = |w: Vector3<f64>| (w.transpose() * inertia * w)[0] * 0.5;
|
|
|
|
// Drag sees the face turned into the flow, which is the velocity direction
|
|
// taken into body components.
|
|
let presented_area = |s: &State<f64, 14, Full>| {
|
|
let incident = attitude::to_dcm(s.attitude()) * -s.velocity().normalize();
|
|
spacecraft.surface.drag_coefficient_area(incident)
|
|
};
|
|
|
|
println!("epoch: {epoch}");
|
|
println!(
|
|
"orbit: a = {:.1} km, i = {:.1} deg ({:?})",
|
|
orbit.semi_major_axis / 1e3,
|
|
orbit.inclination.to_degrees(),
|
|
orbit.degeneracy
|
|
);
|
|
println!("period: {}", Duration::from_seconds(period));
|
|
println!(
|
|
"inertia: diag({:.1}, {:.1}, {:.1}) kg m^2",
|
|
inertia[(0, 0)],
|
|
inertia[(1, 1)],
|
|
inertia[(2, 2)]
|
|
);
|
|
println!(
|
|
"steps: {}, shadow-set switches: {switches}",
|
|
solution.times.len() - 1 - switches
|
|
);
|
|
let final_rate = state.angular_rate();
|
|
println!(
|
|
"rate: {:.4} deg/s about ({:.2}, {:.2}, {:.2})",
|
|
rate.norm().to_degrees(),
|
|
rate[0] / rate.norm(),
|
|
rate[1] / rate.norm(),
|
|
rate[2] / rate.norm()
|
|
);
|
|
println!(
|
|
" -> {:.4} deg/s about ({:.2}, {:.2}, {:.2})",
|
|
final_rate.norm().to_degrees(),
|
|
final_rate[0] / final_rate.norm(),
|
|
final_rate[1] / final_rate.norm(),
|
|
final_rate[2] / final_rate.norm()
|
|
);
|
|
// An asymmetric body under no torque changes both its axis and its rate. What
|
|
// it keeps is the angular momentum and the rotational energy.
|
|
println!(
|
|
"momentum drift: {:.3e} relative",
|
|
(momentum(final_rate) - momentum(rate)).abs() / momentum(rate)
|
|
);
|
|
println!(
|
|
"energy drift: {:.3e} relative",
|
|
(energy(final_rate) - energy(rate)).abs() / energy(rate)
|
|
);
|
|
println!(
|
|
"presented area: {:.3} -> {:.3} m^2 of C_d A",
|
|
presented_area(&initial),
|
|
presented_area(&state)
|
|
);
|
|
println!(
|
|
"orbit, against a translational propagation of the same arc: {:.3e} m",
|
|
(state.position() - reference.position).norm()
|
|
);
|
|
}
|