1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
// SPDX-License-Identifier: AGPL-3.0-only
//! IAU 2009 Mars body-fixed orientation — the rotation from Mars-centred inertial (≈J2000/ICRF)
//! to the Mars body-fixed frame.
//!
//! Mars has a constant pole (the IAU 2009/WGCCRE Mars model carries no periodic nutation terms,
//! unlike the Moon's libration series), so the orientation is fully described by the fixed pole
//! right ascension `α₀` and declination `δ₀` plus the linearly-advancing prime meridian
//!
//! ```text
//! α₀ = 317.681° (constant)
//! δ₀ = 52.886° (constant)
//! W = 176.630° + 350.89198226°/day · (jd_tdb − JD_J2000)
//! ```
//!
//! The body-fixed rotation is the standard IAU 3-1-3 sequence
//! `R = R_z(W)·R_x(90°−δ₀)·R_z(90°+α₀)`, so `r_bodyfixed = R · r_inertial` — the same
//! inertial→body-fixed convention [`crate::lunar_frame::icrf_to_iau_moon`] uses for the Moon and
//! [`crate::cio::gcrs_to_itrs_matrix`] uses for the Earth, so the geopotential evaluation path
//! ([`crate::gravity_sh`]) is reused unchanged: rotate the inertial position in, evaluate the
//! Mars-fixed field, rotate the acceleration back.
//!
//! The pole and prime-meridian constants are read from the [`Body`] passed in
//! ([`Body::mars`](crate::body::Body::mars) carries the IAU 2009 Mars values), so this works for
//! any rotating body whose `pole_ra0`/`pole_dec0`/`prime_w0`/`prime_w_dot` are populated.
//!
//! ## Scope (honest)
//!
//! This realizes the *mean* IAU 2009 Mars orientation (the published `BODY499_*` linear model,
//! NAIF generic PCK). Mars carries no IAU periodic-pole terms, so unlike the lunar frame there is
//! no libration/nutation series to approximate; the residual is the difference between the linear
//! IAU model and a numerically-integrated Mars orientation, well below the gravity-field accuracy
//! at this degree. TDB is approximated by TT (sub-2-ms agreement, a sub-µas orientation error).
use crate::body::Body;
use crate::precession::{mat_vec, transpose, Mat3};
type Vec3 = [f64; 3];
/// Julian Date of the J2000.0 epoch — the reference instant the prime-meridian advance `Ẇ` is
/// measured from.
const JD_J2000: f64 = 2_451_545.0;
/// Frame rotation about +z by `theta` (rotates the coordinate axes; the
/// [`crate::lunar_frame`]/`crate::lunar::rot3` convention).
fn rz(theta: f64) -> Mat3 {
let (s, c) = theta.sin_cos();
[[c, s, 0.0], [-s, c, 0.0], [0.0, 0.0, 1.0]]
}
/// Frame rotation about +x by `theta`.
fn rx(theta: f64) -> Mat3 {
let (s, c) = theta.sin_cos();
[[1.0, 0.0, 0.0], [0.0, c, s], [0.0, -s, c]]
}
/// 3×3 matrix product `a·b`.
fn matmul(a: &Mat3, b: &Mat3) -> Mat3 {
let mut m = [[0.0; 3]; 3];
for (i, mrow) in m.iter_mut().enumerate() {
for (j, e) in mrow.iter_mut().enumerate() {
*e = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
m
}
/// The Mars prime-meridian angle `W` (radians) at `jd_tdb`, from the body's epoch value
/// `prime_w0` and rate `prime_w_dot` (rad/day): `W = W₀ + Ẇ·(jd_tdb − J2000)`. Not reduced to
/// `[0, 2π)` (the rotation builder takes sin/cos).
pub fn mars_prime_meridian(body: &Body, jd_tdb: f64) -> f64 {
body.prime_w0 + body.prime_w_dot * (jd_tdb - JD_J2000)
}
/// The rotation from Mars-centred inertial (≈J2000/ICRF) to the Mars body-fixed frame at
/// `jd_tdb`, as the IAU 3-1-3 sequence `R_z(W)·R_x(90°−δ₀)·R_z(90°+α₀)` built from the body's
/// fixed pole `(pole_ra0, pole_dec0)` and the prime meridian `W` of [`mars_prime_meridian`].
/// Apply with [`crate::precession::mat_vec`]: `r_bodyfixed = R · r_inertial`; the inverse
/// (body-fixed → inertial) is its transpose. The rows of `R` are the Mars-fixed axes expressed
/// in ICRF, so row 2 is the Mars pole `(cos δ₀ cos α₀, cos δ₀ sin α₀, sin δ₀)`.
pub fn iau_mars_rotation(body: &Body, jd_tdb: f64) -> Mat3 {
let ra = body.pole_ra0;
let dec = body.pole_dec0;
let w = mars_prime_meridian(body, jd_tdb);
let half_pi = std::f64::consts::FRAC_PI_2;
// R = Rz(W) · Rx(90°−δ) · Rz(90°+α).
matmul(&matmul(&rz(w), &rx(half_pi - dec)), &rz(half_pi + ra))
}
/// Rotate an inertial position (≈J2000/ICRF, Mars-centred) into the Mars body-fixed frame at
/// `jd_tdb`: `r_bodyfixed = R · r_inertial`.
pub fn inertial_to_bodyfixed(r_inertial: Vec3, body: &Body, jd_tdb: f64) -> Vec3 {
mat_vec(&iau_mars_rotation(body, jd_tdb), r_inertial)
}
/// Rotate a Mars body-fixed vector back into the inertial frame at `jd_tdb`:
/// `r_inertial = Rᵀ · r_bodyfixed`. Used to carry a body-fixed gravity acceleration back to the
/// inertial integration frame.
pub fn bodyfixed_to_inertial(r_bodyfixed: Vec3, body: &Body, jd_tdb: f64) -> Vec3 {
mat_vec(&transpose(&iau_mars_rotation(body, jd_tdb)), r_bodyfixed)
}
#[cfg(test)]
mod tests {
use super::*;
fn norm(v: Vec3) -> f64 {
(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt()
}
/// A body-fixed orientation must be a proper rotation: Rᵀ·R = I (to 1e-12) and det = +1.
#[test]
fn rotation_is_orthonormal_and_proper() {
let body = Body::mars();
let jd = 2_459_580.5; // 2022-01-01 00:00
let m = iau_mars_rotation(&body, jd);
let mt = transpose(&m);
let prod = matmul(&mt, &m);
for (i, row) in prod.iter().enumerate() {
for (j, &e) in row.iter().enumerate() {
let want = if i == j { 1.0 } else { 0.0 };
assert!(
(e - want).abs() < 1e-12,
"RᵀR[{i}][{j}] = {e} (want {want})"
);
}
}
let det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
- m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
+ m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
assert!(
(det - 1.0).abs() < 1e-12,
"det = {det} (want +1, proper rotation)"
);
}
/// Applying the rotation then its transpose returns the original vector (round-trip), and the
/// helpers `inertial_to_bodyfixed`/`bodyfixed_to_inertial` are exact inverses.
#[test]
fn round_trip_recovers_the_original_vector() {
let body = Body::mars();
let jd = 2_459_580.5;
let r = [1.50e6, -0.70e6, 0.95e6];
let bf = inertial_to_bodyfixed(r, &body, jd);
let back = bodyfixed_to_inertial(bf, &body, jd);
for k in 0..3 {
assert!(
(back[k] - r[k]).abs() < 1e-6,
"axis {k}: round-trip {} vs original {}",
back[k],
r[k]
);
}
// The rotation preserves length (it is an isometry).
assert!(
(norm(bf) - norm(r)).abs() < 1e-6,
"rotation must preserve length"
);
}
/// The matrix is inertial→body-fixed, so its row 2 is the Mars pole expressed in ICRF:
/// `(cos δ₀ cos α₀, cos δ₀ sin α₀, sin δ₀)`. This ties the 3-1-3 assembly to `(α₀, δ₀)`, and
/// the pole is constant in time (Mars has no IAU periodic-pole terms).
#[test]
fn pole_row_matches_ra_dec_and_is_constant() {
let body = Body::mars();
let (ra, dec) = (body.pole_ra0, body.pole_dec0);
let want = [dec.cos() * ra.cos(), dec.cos() * ra.sin(), dec.sin()];
for &jd in &[JD_J2000, 2_459_580.5, JD_J2000 + 1000.0] {
let m = iau_mars_rotation(&body, jd);
let pole = [m[2][0], m[2][1], m[2][2]];
for k in 0..3 {
assert!(
(pole[k] - want[k]).abs() < 1e-12,
"jd {jd} pole[{k}] {} vs {}",
pole[k],
want[k]
);
}
}
// Physical sanity: the IAU 2009 Mars pole is RA ≈ 317.68°, Dec ≈ 52.89°.
assert!((ra.to_degrees() - 317.681).abs() < 1e-3);
assert!((dec.to_degrees() - 52.886).abs() < 1e-3);
}
/// `W` advances by exactly `prime_w_dot` per day, so the rotation at `jd` and `jd+1` differ by
/// the expected Mars spin angle (≈ 350.89198°/day, just over one sidereal Mars day).
#[test]
fn prime_meridian_advances_by_w_dot_per_day() {
let body = Body::mars();
let jd = 2_459_580.5;
let dw =
(mars_prime_meridian(&body, jd + 1.0) - mars_prime_meridian(&body, jd)).to_degrees();
assert!(
(dw - 350.891_982_26).abs() < 1e-9,
"W rate {dw}°/day (want 350.89198226)"
);
// And the rotation itself reorients by that angle: a vector fixed in inertial space, seen
// in the body-fixed frame, rotates about the pole by Ẇ between jd and jd+1. Equivalently,
// applying R(jd+1)·R(jd)ᵀ to a body-fixed vector is a rotation by Ẇ about +z (the pole).
let spin = matmul(
&iau_mars_rotation(&body, jd + 1.0),
&transpose(&iau_mars_rotation(&body, jd)),
);
// The trace of a rotation by angle θ about any axis is 1 + 2cos θ.
let trace = spin[0][0] + spin[1][1] + spin[2][2];
let theta = ((trace - 1.0) / 2.0).acos().to_degrees();
// 350.892° is equivalent to a −9.108° rotation; acos gives the unsigned 9.108°.
let expected = (360.0 - 350.891_982_26_f64).abs();
assert!(
(theta - expected).abs() < 1e-6,
"spin angle {theta}° between jd and jd+1 (want {expected})"
);
}
}