use core::marker::PhantomData;
use crate::state::{TranslationalState, TranslationalStateTyped};
use astrodyn_math::quaternion::NORM_LIMIT;
use astrodyn_math::JeodQuat;
use astrodyn_quantities::aliases::AngularVelocity;
use astrodyn_quantities::body_attitude::BodyAttitude;
use astrodyn_quantities::frame::{BodyFrame, Frame, RootInertial, Vehicle};
use glam::{DMat3, DVec3};
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RotationalState {
pub quaternion: JeodQuat,
pub ang_vel_body: DVec3,
}
impl Default for RotationalState {
fn default() -> Self {
Self {
quaternion: JeodQuat::identity(),
ang_vel_body: DVec3::ZERO,
}
}
}
#[derive(Debug, Clone, Copy, PartialEq, Default)]
pub struct SixDofState {
pub trans: TranslationalState,
pub rot: RotationalState,
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SixDofStateTyped<V: Vehicle, F: Frame = RootInertial> {
pub trans: TranslationalStateTyped<F>,
pub rot: RotationalStateTyped<V>,
}
impl<V: Vehicle, F: Frame> Default for SixDofStateTyped<V, F> {
#[inline]
fn default() -> Self {
Self {
trans: TranslationalStateTyped::<F>::default(),
rot: RotationalStateTyped::<V>::default(),
}
}
}
impl<V: Vehicle, F: Frame> SixDofStateTyped<V, F> {
#[inline]
pub fn relabel_to<F2: Frame>(self) -> SixDofStateTyped<V, F2> {
SixDofStateTyped {
trans: self.trans.relabel_to::<F2>(),
rot: self.rot,
}
}
#[inline]
pub fn to_untyped(&self) -> SixDofState {
SixDofState {
trans: self.trans.to_untyped(),
rot: self.rot.to_untyped(),
}
}
#[inline]
pub fn from_untyped_unchecked(s: &SixDofState) -> Self {
Self {
trans: TranslationalStateTyped::<F>::from_untyped_unchecked(&s.trans),
rot: RotationalStateTyped::<V>::from_untyped_unchecked(&s.rot),
}
}
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RotationalStateTyped<V: Vehicle> {
pub q_inertial_body: BodyAttitude<V>,
pub ang_vel_body: AngularVelocity<BodyFrame<V>>,
_v: PhantomData<V>,
}
impl<V: Vehicle> Default for RotationalStateTyped<V> {
#[inline]
fn default() -> Self {
Self {
q_inertial_body: BodyAttitude::identity(),
ang_vel_body: AngularVelocity::<BodyFrame<V>>::zero(),
_v: PhantomData,
}
}
}
impl<V: Vehicle> RotationalStateTyped<V> {
#[inline]
pub fn new(
q_inertial_body: BodyAttitude<V>,
ang_vel_body: AngularVelocity<BodyFrame<V>>,
) -> Self {
Self {
q_inertial_body,
ang_vel_body,
_v: PhantomData,
}
}
#[inline]
pub fn to_untyped(&self) -> RotationalState {
RotationalState {
quaternion: self.q_inertial_body.to_jeod_quat(),
ang_vel_body: self.ang_vel_body.raw_si(),
}
}
#[inline]
pub fn from_untyped_unchecked(s: &RotationalState) -> Self {
Self {
q_inertial_body: BodyAttitude::<V>::from_jeod_quat(s.quaternion),
ang_vel_body: AngularVelocity::<BodyFrame<V>>::from_raw_si(s.ang_vel_body),
_v: PhantomData,
}
}
}
pub fn compute_rotational_acceleration(
inertia: &DMat3,
inverse_inertia: &DMat3,
ang_vel: DVec3,
extern_torq_body: DVec3,
) -> DVec3 {
let ang_mom = *inertia * ang_vel;
let inertial_torq = ang_vel.cross(ang_mom);
let torque_body = extern_torq_body - inertial_torq;
let rot_accel = *inverse_inertia * torque_body;
zero_small(rot_accel)
}
fn zero_small(v: DVec3) -> DVec3 {
const THRESHOLD: f64 = 1e-20;
DVec3::new(
if v.x.abs() < THRESHOLD { 0.0 } else { v.x },
if v.y.abs() < THRESHOLD { 0.0 } else { v.y },
if v.z.abs() < THRESHOLD { 0.0 } else { v.z },
)
}
pub fn compute_left_quat_deriv(q: &JeodQuat, ang_vel: DVec3) -> [f64; 4] {
let mhang_vel = -0.5 * ang_vel;
let qv = q.vector();
let qs = q.scalar();
let qdot_s = -qv.dot(mhang_vel);
let qdot_v = qs * mhang_vel + mhang_vel.cross(qv);
[qdot_s, qdot_v.x, qdot_v.y, qdot_v.z]
}
pub fn normalize_integ(q: &mut JeodQuat) {
let qmagsq = q.norm_sq();
assert!(
qmagsq > 0.0,
"normalize_integ called with zero-magnitude quaternion (norm_sq == 0.0)"
);
let diff1 = 1.0 - qmagsq;
let fact = if diff1 > -NORM_LIMIT && diff1 < NORM_LIMIT {
2.0 / (1.0 + qmagsq)
} else {
1.0 / qmagsq.sqrt()
};
q.data[0] *= fact;
q.data[1] *= fact;
q.data[2] *= fact;
q.data[3] *= fact;
}
#[cfg(test)]
mod tests {
use super::*;
use std::f64::consts::PI;
const TOL: f64 = 1e-14;
#[test]
fn euler_zero_torque_zero_omega() {
let inertia = DMat3::from_diagonal(DVec3::new(10.0, 20.0, 30.0));
let inv_inertia = DMat3::from_diagonal(DVec3::new(0.1, 0.05, 1.0 / 30.0));
let alpha =
compute_rotational_acceleration(&inertia, &inv_inertia, DVec3::ZERO, DVec3::ZERO);
assert!(
alpha.length() < TOL,
"Expected zero rot_accel, got {:?}",
alpha,
);
}
#[test]
fn euler_zero_torque_principal_axis() {
let inertia = DMat3::from_diagonal(DVec3::new(10.0, 20.0, 30.0));
let inv_inertia = DMat3::from_diagonal(DVec3::new(0.1, 0.05, 1.0 / 30.0));
let omega = DVec3::new(0.0, 0.0, 5.0);
let alpha = compute_rotational_acceleration(&inertia, &inv_inertia, omega, DVec3::ZERO);
assert!(
alpha.length() < TOL,
"Expected zero rot_accel for principal-axis spin, got {:?}",
alpha,
);
}
#[test]
fn euler_zero_torque_off_axis() {
let inertia = DMat3::from_diagonal(DVec3::new(10.0, 20.0, 30.0));
let inv_inertia = DMat3::from_diagonal(DVec3::new(0.1, 0.05, 1.0 / 30.0));
let omega = DVec3::new(1.0, 2.0, 3.0);
let alpha = compute_rotational_acceleration(&inertia, &inv_inertia, omega, DVec3::ZERO);
let expected = DVec3::new(-6.0, 3.0, -20.0 / 30.0);
let diff = (alpha - expected).length();
assert!(
diff < 1e-12,
"Off-axis gyroscopic coupling: expected {:?}, got {:?}, diff={}",
expected,
alpha,
diff,
);
}
#[test]
fn quat_deriv_zero_omega() {
let q = JeodQuat::left_quat_from_eigen_rotation(0.5, DVec3::Z);
let qdot = compute_left_quat_deriv(&q, DVec3::ZERO);
for (i, &val) in qdot.iter().enumerate() {
assert!(
val.abs() < TOL,
"qdot[{}] = {} should be zero for zero angular velocity",
i,
val,
);
}
}
#[test]
fn quat_deriv_z_rotation_identity() {
let q = JeodQuat::identity();
let omega_z = 0.1;
let qdot = compute_left_quat_deriv(&q, DVec3::new(0.0, 0.0, omega_z));
assert!(
qdot[0].abs() < TOL,
"qdot scalar should be ~0, got {}",
qdot[0],
);
assert!(qdot[1].abs() < TOL, "qdot vx should be ~0, got {}", qdot[1],);
assert!(qdot[2].abs() < TOL, "qdot vy should be ~0, got {}", qdot[2],);
let expected_vz = -omega_z * 0.5;
assert!(
(qdot[3] - expected_vz).abs() < TOL,
"qdot vz expected {}, got {}",
expected_vz,
qdot[3],
);
}
#[test]
fn quat_deriv_z_rotation_nontrivial() {
let angle = PI / 6.0; let q = JeodQuat::left_quat_from_eigen_rotation(angle, DVec3::Z);
let omega_z = 2.0;
let qdot = compute_left_quat_deriv(&q, DVec3::new(0.0, 0.0, omega_z));
let half = angle * 0.5;
let expected_s = -half.sin();
let expected_vz = -half.cos();
assert!(
(qdot[0] - expected_s).abs() < 1e-12,
"qdot scalar expected {}, got {}",
expected_s,
qdot[0],
);
assert!(
qdot[1].abs() < 1e-12,
"qdot vx should be ~0, got {}",
qdot[1],
);
assert!(
qdot[2].abs() < 1e-12,
"qdot vy should be ~0, got {}",
qdot[2],
);
assert!(
(qdot[3] - expected_vz).abs() < 1e-12,
"qdot vz expected {}, got {}",
expected_vz,
qdot[3],
);
}
#[test]
fn normalize_integ_near_unity() {
let mut q = JeodQuat::new(1.0 + 1e-10, 0.0, 0.0, 0.0);
normalize_integ(&mut q);
let norm_err = (q.norm_sq() - 1.0).abs();
assert!(
norm_err < 1e-14,
"After normalize_integ, |q|^2 - 1 = {}",
norm_err,
);
assert!(
(q.scalar() - 1.0).abs() < 1e-9,
"scalar should be ~1.0, got {}",
q.scalar(),
);
}
#[test]
fn normalize_integ_far_from_unity() {
let mut q = JeodQuat::new(2.0, 0.0, 0.0, 0.0);
normalize_integ(&mut q);
let norm_err = (q.norm_sq() - 1.0).abs();
assert!(
norm_err < 1e-14,
"After normalize_integ, |q|^2 - 1 = {}",
norm_err,
);
assert!(
(q.scalar() - 1.0).abs() < 1e-14,
"scalar should be 1.0, got {}",
q.scalar(),
);
}
#[test]
fn normalize_integ_preserves_sign() {
let mut q = JeodQuat::new(-1.0, 0.0, 0.0, 0.0);
normalize_integ(&mut q);
assert!(
q.scalar() < 0.0,
"normalize_integ should NOT flip sign, scalar = {}",
q.scalar(),
);
let norm_err = (q.norm_sq() - 1.0).abs();
assert!(
norm_err < 1e-14,
"After normalize_integ, |q|^2 - 1 = {}",
norm_err,
);
}
#[test]
fn normalize_integ_nontrivial() {
let mut q = JeodQuat::new(3.0, 4.0, 0.0, 0.0);
normalize_integ(&mut q);
let norm_err = (q.norm_sq() - 1.0).abs();
assert!(
norm_err < 1e-14,
"After normalize_integ, |q|^2 - 1 = {}",
norm_err,
);
assert!(
(q.scalar() - 0.6).abs() < 1e-14,
"scalar should be 0.6, got {}",
q.scalar(),
);
assert!(
(q.data[1] - 0.8).abs() < 1e-14,
"vx should be 0.8, got {}",
q.data[1],
);
}
#[test]
fn typed_rotational_state_round_trips() {
use astrodyn_quantities::aliases::AngularVelocity;
use astrodyn_quantities::frame::TestVehicle;
let q = JeodQuat::left_quat_from_eigen_rotation(0.7, DVec3::Z);
let ang_vel = DVec3::new(0.01, 0.02, 0.03);
let typed = RotationalStateTyped::<TestVehicle>::new(
BodyAttitude::from_jeod_quat(q),
AngularVelocity::<BodyFrame<TestVehicle>>::from_raw_si(ang_vel),
);
assert_eq!(typed.q_inertial_body.to_jeod_quat(), q);
assert_eq!(typed.ang_vel_body.raw_si(), ang_vel);
}
#[test]
fn typed_rotational_default_is_identity() {
use astrodyn_quantities::frame::TestVehicle;
let s = RotationalStateTyped::<TestVehicle>::default();
assert_eq!(s.q_inertial_body.to_jeod_quat(), JeodQuat::identity());
assert_eq!(s.ang_vel_body.raw_si(), DVec3::ZERO);
}
use crate::state::TranslationalState;
use astrodyn_quantities::frame::{RootInertial, TestVehicle};
use astrodyn_quantities::quat::{LeftTransform, NormalizedQuat, ScalarFirst};
use proptest::prelude::*;
fn arb_finite_bounded() -> impl Strategy<Value = f64> {
prop_oneof![
(1.0e-9_f64..1.0e9_f64),
(1.0e-9_f64..1.0e9_f64).prop_map(|x| -x),
]
}
fn arb_dvec3_finite() -> impl Strategy<Value = DVec3> {
(
arb_finite_bounded(),
arb_finite_bounded(),
arb_finite_bounded(),
)
.prop_map(|(x, y, z)| DVec3::new(x, y, z))
}
fn arb_unit_jeod_quat() -> impl Strategy<Value = JeodQuat> {
(
arb_finite_bounded(),
arb_finite_bounded(),
arb_finite_bounded(),
arb_finite_bounded(),
)
.prop_filter("non-zero norm", |(a, b, c, d)| {
let n2 = a * a + b * b + c * c + d * d;
n2.is_finite() && n2 > 1.0e-18
})
.prop_filter_map("renormalize succeeds", |(a, b, c, d)| {
let q = JeodQuat::from_array([a, b, c, d]);
NormalizedQuat::<ScalarFirst, LeftTransform>::renormalize(q).map(|n| n.inner())
})
}
fn arb_rotational_state() -> impl Strategy<Value = RotationalState> {
(arb_unit_jeod_quat(), arb_dvec3_finite()).prop_map(|(quaternion, ang_vel_body)| {
RotationalState {
quaternion,
ang_vel_body,
}
})
}
fn arb_translational_state_for_six_dof() -> impl Strategy<Value = TranslationalState> {
(arb_dvec3_finite(), arb_dvec3_finite())
.prop_map(|(position, velocity)| TranslationalState { position, velocity })
}
fn arb_six_dof_state() -> impl Strategy<Value = SixDofState> {
(
arb_translational_state_for_six_dof(),
arb_rotational_state(),
)
.prop_map(|(trans, rot)| SixDofState { trans, rot })
}
proptest! {
#[test]
fn round_trip_rotational_untyped_typed_untyped(orig in arb_rotational_state()) {
let typed = RotationalStateTyped::<TestVehicle>::from_untyped_unchecked(&orig);
prop_assert_eq!(typed.to_untyped(), orig);
}
#[test]
fn round_trip_rotational_typed_untyped_typed(orig in arb_rotational_state()) {
let typed = RotationalStateTyped::<TestVehicle>::from_untyped_unchecked(&orig);
let lifted = RotationalStateTyped::<TestVehicle>::from_untyped_unchecked(&typed.to_untyped());
prop_assert_eq!(lifted.to_untyped(), typed.to_untyped());
}
#[test]
fn round_trip_six_dof_untyped_typed_untyped(orig in arb_six_dof_state()) {
let typed = SixDofStateTyped::<TestVehicle, RootInertial>::from_untyped_unchecked(&orig);
prop_assert_eq!(typed.to_untyped(), orig);
}
#[test]
fn round_trip_six_dof_typed_untyped_typed(orig in arb_six_dof_state()) {
let typed = SixDofStateTyped::<TestVehicle, RootInertial>::from_untyped_unchecked(&orig);
let lifted = SixDofStateTyped::<TestVehicle, RootInertial>::from_untyped_unchecked(&typed.to_untyped());
prop_assert_eq!(lifted.to_untyped(), typed.to_untyped());
}
}
}