solar_line_core/
vec3.rs

1/// Generic 3D vector for use with unit newtypes.
2use std::fmt;
3use std::ops::{Add, Neg, Sub};
4
5#[derive(Debug, Clone, Copy, PartialEq)]
6pub struct Vec3<T> {
7    pub x: T,
8    pub y: T,
9    pub z: T,
10}
11
12impl<T: Copy> Vec3<T> {
13    pub fn new(x: T, y: T, z: T) -> Self {
14        Self { x, y, z }
15    }
16}
17
18/// Vec3 of f64-newtypes that support the operations we need.
19/// We implement common operations generically where the inner type supports them.
20impl<T> Add for Vec3<T>
21where
22    T: Add<Output = T>,
23{
24    type Output = Self;
25    fn add(self, rhs: Self) -> Self {
26        Self {
27            x: self.x + rhs.x,
28            y: self.y + rhs.y,
29            z: self.z + rhs.z,
30        }
31    }
32}
33
34impl<T> Sub for Vec3<T>
35where
36    T: Sub<Output = T>,
37{
38    type Output = Self;
39    fn sub(self, rhs: Self) -> Self {
40        Self {
41            x: self.x - rhs.x,
42            y: self.y - rhs.y,
43            z: self.z - rhs.z,
44        }
45    }
46}
47
48impl<T> Neg for Vec3<T>
49where
50    T: Neg<Output = T>,
51{
52    type Output = Self;
53    fn neg(self) -> Self {
54        Self {
55            x: -self.x,
56            y: -self.y,
57            z: -self.z,
58        }
59    }
60}
61
62impl<T> Vec3<T>
63where
64    T: Copy + std::ops::Mul<f64, Output = T>,
65{
66    /// Scalar multiplication
67    pub fn scale(self, s: f64) -> Self {
68        Self {
69            x: self.x * s,
70            y: self.y * s,
71            z: self.z * s,
72        }
73    }
74}
75
76/// Trait for types that wrap an f64 value.
77/// All implementors are `Copy`, so consuming `self` is intentional.
78#[allow(clippy::wrong_self_convention)]
79pub trait AsF64 {
80    fn as_f64(self) -> f64;
81}
82
83impl AsF64 for f64 {
84    fn as_f64(self) -> f64 {
85        self
86    }
87}
88
89impl<T: Copy + AsF64> Vec3<T> {
90    /// Dot product returning f64 (in squared units of T).
91    pub fn dot_raw(self, rhs: Self) -> f64 {
92        self.x.as_f64() * rhs.x.as_f64()
93            + self.y.as_f64() * rhs.y.as_f64()
94            + self.z.as_f64() * rhs.z.as_f64()
95    }
96
97    /// Euclidean norm (magnitude) as f64, in units of T.
98    pub fn norm_raw(self) -> f64 {
99        self.dot_raw(self).sqrt()
100    }
101
102    /// Cross product with same type, returning Vec3<f64> (units are T²).
103    pub fn cross_raw(self, rhs: Self) -> Vec3<f64> {
104        Vec3 {
105            x: self.y.as_f64() * rhs.z.as_f64() - self.z.as_f64() * rhs.y.as_f64(),
106            y: self.z.as_f64() * rhs.x.as_f64() - self.x.as_f64() * rhs.z.as_f64(),
107            z: self.x.as_f64() * rhs.y.as_f64() - self.y.as_f64() * rhs.x.as_f64(),
108        }
109    }
110
111    /// Cross product with a different type, returning Vec3<f64> (units are T·U).
112    pub fn cross_raw_with<U: Copy + AsF64>(self, rhs: Vec3<U>) -> Vec3<f64> {
113        Vec3 {
114            x: self.y.as_f64() * rhs.z.as_f64() - self.z.as_f64() * rhs.y.as_f64(),
115            y: self.z.as_f64() * rhs.x.as_f64() - self.x.as_f64() * rhs.z.as_f64(),
116            z: self.x.as_f64() * rhs.y.as_f64() - self.y.as_f64() * rhs.x.as_f64(),
117        }
118    }
119}
120
121impl Vec3<f64> {
122    /// Normalize to unit vector. Returns zero vector if norm is zero.
123    pub fn normalize(self) -> Self {
124        let n = self.norm_raw();
125        if n < 1e-15 {
126            Self::new(0.0, 0.0, 0.0)
127        } else {
128            self.scale(1.0 / n)
129        }
130    }
131}
132
133impl<T: fmt::Display> fmt::Display for Vec3<T> {
134    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
135        write!(f, "({}, {}, {})", self.x, self.y, self.z)
136    }
137}
138
139// Implement AsF64 for our unit types
140use crate::units::{Km, KmPerSec};
141
142impl AsF64 for Km {
143    fn as_f64(self) -> f64 {
144        self.value()
145    }
146}
147
148impl AsF64 for KmPerSec {
149    fn as_f64(self) -> f64 {
150        self.value()
151    }
152}
153
154#[cfg(test)]
155mod tests {
156    use super::*;
157    use crate::units::{Km, KmPerSec};
158
159    #[test]
160    fn test_vec3_f64_operations() {
161        let a = Vec3::new(1.0_f64, 2.0, 3.0);
162        let b = Vec3::new(4.0_f64, 5.0, 6.0);
163        let sum = a + b;
164        assert_eq!(sum, Vec3::new(5.0, 7.0, 9.0));
165        let diff = b - a;
166        assert_eq!(diff, Vec3::new(3.0, 3.0, 3.0));
167    }
168
169    #[test]
170    fn test_vec3_km_operations() {
171        let a = Vec3::new(Km(1.0), Km(0.0), Km(0.0));
172        let b = Vec3::new(Km(0.0), Km(1.0), Km(0.0));
173        let sum = a + b;
174        assert_eq!(sum, Vec3::new(Km(1.0), Km(1.0), Km(0.0)));
175    }
176
177    #[test]
178    fn test_vec3_dot_and_norm() {
179        let v = Vec3::new(Km(3.0), Km(4.0), Km(0.0));
180        assert!((v.norm_raw() - 5.0).abs() < 1e-15);
181        assert!((v.dot_raw(v) - 25.0).abs() < 1e-15);
182    }
183
184    #[test]
185    fn test_vec3_scale() {
186        let v = Vec3::new(Km(1.0), Km(2.0), Km(3.0));
187        let scaled = v.scale(2.0);
188        assert_eq!(scaled, Vec3::new(Km(2.0), Km(4.0), Km(6.0)));
189    }
190
191    #[test]
192    fn test_vec3_speed_norm() {
193        // ISS orbital velocity roughly ~7.66 km/s
194        let v = Vec3::new(KmPerSec(7.0), KmPerSec(3.0), KmPerSec(1.0));
195        let speed = v.norm_raw();
196        let expected = (49.0 + 9.0 + 1.0_f64).sqrt();
197        assert!((speed - expected).abs() < 1e-15);
198    }
199
200    #[test]
201    fn test_vec3_neg() {
202        let v = Vec3::new(Km(1.0), Km(-2.0), Km(3.0));
203        let neg_v = -v;
204        assert_eq!(neg_v, Vec3::new(Km(-1.0), Km(2.0), Km(-3.0)));
205    }
206
207    #[test]
208    fn test_vec3_cross_product() {
209        // x × y = z
210        let x = Vec3::new(Km(1.0), Km(0.0), Km(0.0));
211        let y = Vec3::new(Km(0.0), Km(1.0), Km(0.0));
212        let z = x.cross_raw(y);
213        assert!((z.x - 0.0).abs() < 1e-15);
214        assert!((z.y - 0.0).abs() < 1e-15);
215        assert!((z.z - 1.0).abs() < 1e-15);
216    }
217
218    #[test]
219    fn test_vec3_cross_product_anticommutative() {
220        let a = Vec3::new(Km(1.0), Km(2.0), Km(3.0));
221        let b = Vec3::new(Km(4.0), Km(5.0), Km(6.0));
222        let ab = a.cross_raw(b);
223        let ba = b.cross_raw(a);
224        assert!((ab.x + ba.x).abs() < 1e-10);
225        assert!((ab.y + ba.y).abs() < 1e-10);
226        assert!((ab.z + ba.z).abs() < 1e-10);
227    }
228
229    #[test]
230    fn test_vec3_normalize() {
231        let v = Vec3::new(3.0_f64, 4.0, 0.0);
232        let n = v.normalize();
233        assert!((n.norm_raw() - 1.0).abs() < 1e-15);
234        assert!((n.x - 0.6).abs() < 1e-15);
235        assert!((n.y - 0.8).abs() < 1e-15);
236    }
237
238    #[test]
239    fn test_vec3_normalize_zero() {
240        let v = Vec3::new(0.0_f64, 0.0, 0.0);
241        let n = v.normalize();
242        assert!((n.norm_raw()).abs() < 1e-15);
243    }
244
245    #[test]
246    fn test_vec3_cross_raw_with_heterogeneous() {
247        // cross_raw_with: Km × KmPerSec (angular momentum direction)
248        let r = Vec3::new(Km(1.0), Km(0.0), Km(0.0));
249        let v = Vec3::new(KmPerSec(0.0), KmPerSec(1.0), KmPerSec(0.0));
250        let h = r.cross_raw_with(v);
251        // (1,0,0) × (0,1,0) = (0,0,1)
252        assert!((h.x - 0.0).abs() < 1e-15);
253        assert!((h.y - 0.0).abs() < 1e-15);
254        assert!((h.z - 1.0).abs() < 1e-15);
255    }
256
257    #[test]
258    fn test_vec3_dot_raw_orthogonal() {
259        // Orthogonal vectors should have zero dot product
260        let a = Vec3::new(Km(1.0), Km(0.0), Km(0.0));
261        let b = Vec3::new(Km(0.0), Km(1.0), Km(0.0));
262        assert!(
263            a.dot_raw(b).abs() < 1e-15,
264            "orthogonal dot product should be 0"
265        );
266    }
267
268    #[test]
269    fn test_vec3_scale_negative() {
270        let v = Vec3::new(Km(1.0), Km(2.0), Km(3.0));
271        let neg = v.scale(-1.0);
272        assert_eq!(neg, Vec3::new(Km(-1.0), Km(-2.0), Km(-3.0)));
273        let zero = v.scale(0.0);
274        assert_eq!(zero, Vec3::new(Km(0.0), Km(0.0), Km(0.0)));
275    }
276}