sidereon_core/astro/math/
mat3.rs

1//! 3x3 matrix operations with exact floating-point rounding control.
2//!
3//! These routines are ported from the C++ inline_rxr, inline_tr, and
4//! inline_mxmxm functions which were compiled under `fp-contract=off`
5//! to match Python/Skyfield arithmetic.  Rust's default arithmetic
6//! operators do not fuse multiply-adds, so no special annotation is
7//! needed -- just avoid `f64::mul_add`.
8
9/// A row-major 3x3 matrix.
10pub type Mat3 = [[f64; 3]; 3];
11
12/// Matrix-vector product: `result = m * v`.
13pub fn mul_vec3(m: &Mat3, v: [f64; 3]) -> [f64; 3] {
14    [
15        m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2],
16        m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2],
17        m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2],
18    ]
19}
20
21/// Standard matrix multiply: `result = a * b`.
22pub fn inline_rxr(a: &Mat3, b: &Mat3) -> Mat3 {
23    let mut w = [[0.0_f64; 3]; 3];
24    for i in 0..3 {
25        for j in 0..3 {
26            let mut s = 0.0_f64;
27            for k in 0..3 {
28                s += a[i][k] * b[k][j];
29            }
30            w[i][j] = s;
31        }
32    }
33    w
34}
35
36/// Matrix transpose: `result = r^T`.
37pub fn inline_tr(r: &Mat3) -> Mat3 {
38    let mut rt = [[0.0_f64; 3]; 3];
39    for i in 0..3 {
40        for j in 0..3 {
41            rt[i][j] = r[j][i];
42        }
43    }
44    rt
45}
46
47/// Triple matrix product with Kahan-compensated summation.
48///
49/// Computes `result = A * B * C` without materialising the intermediate
50/// A*B matrix.  Each output element accumulates all 9 terms
51/// (matching numpy `einsum('ij,jk,kl->il')`) using Kahan summation to
52/// reduce floating-point drift.
53pub fn inline_mxmxm(a: &Mat3, b: &Mat3, c: &Mat3) -> Mat3 {
54    let mut w = [[0.0_f64; 3]; 3];
55    for i in 0..3 {
56        for l in 0..3 {
57            let mut s = 0.0_f64;
58            let mut comp = 0.0_f64; // Kahan compensation
59            for j in 0..3 {
60                for k in 0..3 {
61                    let term = a[i][j] * b[j][k] * c[k][l];
62                    let y = term - comp;
63                    let t = s + y;
64                    comp = (t - s) - y;
65                    s = t;
66                }
67            }
68            w[i][l] = s;
69        }
70    }
71    w
72}
sidereon_core/astro/math/mat3.rs

sidereon_core/astro/math/
mat3.rs
