solar_line_core/
units.rs

1/// Type-safe unit newtypes for orbital mechanics.
2///
3/// All units use km-based system consistent with standard astrodynamics:
4/// - Distance: km
5/// - Speed: km/s
6/// - Time: seconds
7/// - Angles: radians
8/// - Gravitational parameter: km³/s²
9use std::fmt;
10use std::ops::{Add, Div, Mul, Neg, Sub};
11
12macro_rules! unit_newtype {
13    ($name:ident, $unit_str:expr) => {
14        #[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
15        pub struct $name(pub f64);
16
17        impl $name {
18            pub fn value(self) -> f64 {
19                self.0
20            }
21        }
22
23        impl fmt::Display for $name {
24            fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
25                write!(f, "{} {}", self.0, $unit_str)
26            }
27        }
28
29        impl Add for $name {
30            type Output = Self;
31            fn add(self, rhs: Self) -> Self {
32                Self(self.0 + rhs.0)
33            }
34        }
35
36        impl Sub for $name {
37            type Output = Self;
38            fn sub(self, rhs: Self) -> Self {
39                Self(self.0 - rhs.0)
40            }
41        }
42
43        impl Neg for $name {
44            type Output = Self;
45            fn neg(self) -> Self {
46                Self(-self.0)
47            }
48        }
49
50        impl Mul<f64> for $name {
51            type Output = Self;
52            fn mul(self, rhs: f64) -> Self {
53                Self(self.0 * rhs)
54            }
55        }
56
57        impl Mul<$name> for f64 {
58            type Output = $name;
59            fn mul(self, rhs: $name) -> $name {
60                $name(self * rhs.0)
61            }
62        }
63
64        impl Div<f64> for $name {
65            type Output = Self;
66            fn div(self, rhs: f64) -> Self {
67                Self(self.0 / rhs)
68            }
69        }
70
71        impl Div<$name> for $name {
72            type Output = f64;
73            fn div(self, rhs: $name) -> f64 {
74                self.0 / rhs.0
75            }
76        }
77    };
78}
79
80unit_newtype!(Km, "km");
81unit_newtype!(KmPerSec, "km/s");
82unit_newtype!(Seconds, "s");
83unit_newtype!(Radians, "rad");
84unit_newtype!(Mu, "km³/s²");
85
86impl Km {
87    pub fn abs(self) -> Self {
88        Self(self.0.abs())
89    }
90}
91
92impl KmPerSec {
93    pub fn abs(self) -> Self {
94        Self(self.0.abs())
95    }
96}
97
98impl Radians {
99    /// Normalize angle to [0, 2π)
100    pub fn normalize(self) -> Self {
101        let two_pi = std::f64::consts::TAU;
102        let mut v = self.0 % two_pi;
103        if v < 0.0 {
104            v += two_pi;
105        }
106        Self(v)
107    }
108
109    /// Normalize angle to (-π, π]
110    pub fn normalize_signed(self) -> Self {
111        let two_pi = std::f64::consts::TAU;
112        let pi = std::f64::consts::PI;
113        let mut v = self.0 % two_pi;
114        if v > pi {
115            v -= two_pi;
116        } else if v <= -pi {
117            v += two_pi;
118        }
119        Self(v)
120    }
121
122    pub fn sin(self) -> f64 {
123        self.0.sin()
124    }
125
126    pub fn cos(self) -> f64 {
127        self.0.cos()
128    }
129
130    pub fn tan(self) -> f64 {
131        self.0.tan()
132    }
133}
134
135/// Eccentricity with validation.
136/// For elliptical orbits: 0 <= e < 1
137/// For parabolic: e == 1
138/// For hyperbolic: e > 1
139#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
140pub struct Eccentricity(f64);
141
142impl Eccentricity {
143    /// Create a new eccentricity value. Returns None if negative.
144    pub fn new(e: f64) -> Option<Self> {
145        if e >= 0.0 {
146            Some(Self(e))
147        } else {
148            None
149        }
150    }
151
152    /// Create eccentricity for an elliptical orbit (0 <= e < 1).
153    /// Returns None if out of range.
154    pub fn elliptical(e: f64) -> Option<Self> {
155        if (0.0..1.0).contains(&e) {
156            Some(Self(e))
157        } else {
158            None
159        }
160    }
161
162    pub fn value(self) -> f64 {
163        self.0
164    }
165
166    pub fn is_circular(&self) -> bool {
167        self.0 == 0.0
168    }
169
170    pub fn is_elliptical(&self) -> bool {
171        self.0 < 1.0
172    }
173
174    pub fn is_parabolic(&self) -> bool {
175        (self.0 - 1.0).abs() < f64::EPSILON
176    }
177
178    pub fn is_hyperbolic(&self) -> bool {
179        self.0 > 1.0
180    }
181}
182
183impl fmt::Display for Eccentricity {
184    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
185        write!(f, "e={}", self.0)
186    }
187}
188
189#[cfg(test)]
190mod tests {
191    use super::*;
192    use std::f64::consts::{FRAC_PI_2, PI, TAU};
193
194    #[test]
195    fn test_km_arithmetic() {
196        let a = Km(100.0);
197        let b = Km(50.0);
198        assert_eq!((a + b).value(), 150.0);
199        assert_eq!((a - b).value(), 50.0);
200        assert_eq!((a * 2.0).value(), 200.0);
201        assert_eq!((3.0 * b).value(), 150.0);
202        assert_eq!((a / 2.0).value(), 50.0);
203        assert_eq!(a / b, 2.0); // dimensionless ratio
204    }
205
206    #[test]
207    fn test_radians_normalize() {
208        // Already normalized
209        assert!((Radians(1.0).normalize().value() - 1.0).abs() < 1e-15);
210
211        // Negative angle
212        let r = Radians(-FRAC_PI_2).normalize();
213        assert!((r.value() - (TAU - FRAC_PI_2)).abs() < 1e-15);
214
215        // Greater than 2π
216        let r = Radians(TAU + 1.0).normalize();
217        assert!((r.value() - 1.0).abs() < 1e-14);
218    }
219
220    #[test]
221    fn test_radians_normalize_signed() {
222        // Positive within range
223        assert!((Radians(1.0).normalize_signed().value() - 1.0).abs() < 1e-15);
224
225        // Should wrap to negative
226        let r = Radians(PI + 0.5).normalize_signed();
227        assert!((r.value() - (-PI + 0.5)).abs() < 1e-14);
228
229        // Already negative, within range
230        let r = Radians(-1.0).normalize_signed();
231        assert!((r.value() - (-1.0)).abs() < 1e-15);
232    }
233
234    #[test]
235    fn test_radians_trig() {
236        let r = Radians(FRAC_PI_2);
237        assert!((r.sin() - 1.0).abs() < 1e-15);
238        assert!(r.cos().abs() < 1e-15);
239    }
240
241    #[test]
242    fn test_eccentricity_validation() {
243        assert!(Eccentricity::new(-0.1).is_none());
244        assert!(Eccentricity::new(0.0).is_some());
245        assert!(Eccentricity::new(0.5).is_some());
246        assert!(Eccentricity::new(1.5).is_some()); // hyperbolic is valid
247
248        assert!(Eccentricity::elliptical(-0.1).is_none());
249        assert!(Eccentricity::elliptical(0.5).is_some());
250        assert!(Eccentricity::elliptical(1.0).is_none()); // not elliptical
251        assert!(Eccentricity::elliptical(1.5).is_none());
252    }
253
254    #[test]
255    fn test_eccentricity_classification() {
256        let circular = Eccentricity::new(0.0).unwrap();
257        assert!(circular.is_circular());
258        assert!(circular.is_elliptical());
259
260        let elliptical = Eccentricity::new(0.5).unwrap();
261        assert!(!elliptical.is_circular());
262        assert!(elliptical.is_elliptical());
263
264        let parabolic = Eccentricity::new(1.0).unwrap();
265        assert!(parabolic.is_parabolic());
266
267        let hyperbolic = Eccentricity::new(1.5).unwrap();
268        assert!(hyperbolic.is_hyperbolic());
269    }
270
271    #[test]
272    fn test_display() {
273        assert_eq!(format!("{}", Km(42164.0)), "42164 km");
274        assert_eq!(format!("{}", KmPerSec(7.784)), "7.784 km/s");
275        assert_eq!(format!("{}", Seconds(3600.0)), "3600 s");
276    }
277
278    #[test]
279    fn test_negation() {
280        assert_eq!((-KmPerSec(3.0)).value(), -3.0);
281        assert_eq!((-Km(100.0)).value(), -100.0);
282    }
283
284    // --- Seconds type tests ---
285
286    #[test]
287    fn seconds_arithmetic() {
288        let a = Seconds(3600.0);
289        let b = Seconds(1800.0);
290        assert_eq!((a + b).value(), 5400.0);
291        assert_eq!((a - b).value(), 1800.0);
292        assert_eq!((a * 2.0).value(), 7200.0);
293        assert_eq!((2.0 * b).value(), 3600.0);
294        assert_eq!((a / 2.0).value(), 1800.0);
295        assert_eq!(a / b, 2.0);
296    }
297
298    #[test]
299    fn seconds_negation() {
300        assert_eq!((-Seconds(60.0)).value(), -60.0);
301        assert_eq!((-Seconds(-30.0)).value(), 30.0);
302    }
303
304    #[test]
305    fn seconds_display() {
306        assert_eq!(format!("{}", Seconds(86400.0)), "86400 s");
307        assert_eq!(format!("{}", Seconds(0.001)), "0.001 s");
308    }
309
310    // --- Mu type tests ---
311
312    #[test]
313    fn mu_arithmetic() {
314        let a = Mu(3.986e5);
315        let b = Mu(1.327e11);
316        assert_eq!((a + b).value(), 3.986e5 + 1.327e11);
317        assert_eq!((b - a).value(), 1.327e11 - 3.986e5);
318        assert_eq!((a * 2.0).value(), 2.0 * 3.986e5);
319        assert_eq!((0.5 * a).value(), 0.5 * 3.986e5);
320        assert_eq!((a / 2.0).value(), 3.986e5 / 2.0);
321        assert!((b / a - 1.327e11 / 3.986e5).abs() < 1e-6);
322    }
323
324    #[test]
325    fn mu_negation() {
326        assert_eq!((-Mu(100.0)).value(), -100.0);
327    }
328
329    #[test]
330    fn mu_display() {
331        assert_eq!(format!("{}", Mu(398600.4418)), "398600.4418 km³/s²");
332    }
333
334    // --- KmPerSec arithmetic tests ---
335
336    #[test]
337    fn km_per_sec_arithmetic() {
338        let a = KmPerSec(10.0);
339        let b = KmPerSec(3.0);
340        assert_eq!((a + b).value(), 13.0);
341        assert_eq!((a - b).value(), 7.0);
342        assert_eq!((a * 0.5).value(), 5.0);
343        assert_eq!((4.0 * b).value(), 12.0);
344        assert_eq!((a / 5.0).value(), 2.0);
345        assert!((a / b - 10.0 / 3.0).abs() < 1e-15);
346    }
347
348    // --- Radians arithmetic tests ---
349
350    #[test]
351    fn radians_arithmetic() {
352        let a = Radians(PI);
353        let b = Radians(FRAC_PI_2);
354        assert!((a + b).value() - (PI + FRAC_PI_2) < 1e-15);
355        assert!((a - b).value() - FRAC_PI_2 < 1e-15);
356        assert_eq!((-a).value(), -PI);
357        assert!((a * 2.0).value() - TAU < 1e-15);
358        assert!((2.0 * b).value() - PI < 1e-15);
359        assert!((a / 2.0).value() - FRAC_PI_2 < 1e-15);
360        assert!((a / b - 2.0).abs() < 1e-15);
361    }
362
363    #[test]
364    fn radians_display() {
365        assert_eq!(format!("{}", Radians(3.14159)), "3.14159 rad");
366    }
367
368    // --- Km::abs / KmPerSec::abs edge cases ---
369
370    #[test]
371    fn km_abs() {
372        assert_eq!(Km(-42.0).abs().value(), 42.0);
373        assert_eq!(Km(42.0).abs().value(), 42.0);
374        assert_eq!(Km(0.0).abs().value(), 0.0);
375    }
376
377    #[test]
378    fn km_per_sec_abs() {
379        assert_eq!(KmPerSec(-7.8).abs().value(), 7.8);
380        assert_eq!(KmPerSec(7.8).abs().value(), 7.8);
381        assert_eq!(KmPerSec(0.0).abs().value(), 0.0);
382    }
383
384    // --- Radians::tan() ---
385
386    #[test]
387    fn radians_tan() {
388        // tan(π/4) = 1
389        assert!((Radians(std::f64::consts::FRAC_PI_4).tan() - 1.0).abs() < 1e-15);
390        // tan(0) = 0
391        assert!(Radians(0.0).tan().abs() < 1e-15);
392    }
393
394    // --- normalize_signed boundary at exactly -π ---
395
396    #[test]
397    fn radians_normalize_signed_at_neg_pi() {
398        // At exactly -π, should wrap to positive π (since range is (-π, π])
399        let r = Radians(-PI).normalize_signed();
400        assert!((r.value() - PI).abs() < 1e-14);
401    }
402
403    #[test]
404    fn radians_normalize_signed_at_pi() {
405        // At exactly π, should stay at π
406        let r = Radians(PI).normalize_signed();
407        assert!((r.value() - PI).abs() < 1e-14);
408    }
409
410    #[test]
411    fn radians_normalize_large_negative() {
412        // -3π should normalize to -π + 2π = π ... but modulo: -3π % 2π = -π → +2π = π
413        let r = Radians(-3.0 * PI).normalize_signed();
414        assert!((r.value() - PI).abs() < 1e-13);
415    }
416
417    #[test]
418    fn radians_normalize_zero() {
419        assert_eq!(Radians(0.0).normalize().value(), 0.0);
420        assert_eq!(Radians(0.0).normalize_signed().value(), 0.0);
421    }
422
423    // --- Eccentricity edge cases ---
424
425    #[test]
426    fn eccentricity_value_and_display() {
427        let e = Eccentricity::new(0.0167).unwrap();
428        assert!((e.value() - 0.0167).abs() < 1e-15);
429        assert_eq!(format!("{}", e), "e=0.0167");
430    }
431
432    #[test]
433    fn eccentricity_elliptical_boundary() {
434        // Exactly 0 is valid elliptical (circular)
435        assert!(Eccentricity::elliptical(0.0).is_some());
436        // Just below 1 is valid elliptical
437        assert!(Eccentricity::elliptical(0.999999).is_some());
438        // Exactly 1 is NOT elliptical
439        assert!(Eccentricity::elliptical(1.0).is_none());
440    }
441
442    #[test]
443    fn eccentricity_partial_ord() {
444        let low = Eccentricity::new(0.1).unwrap();
445        let high = Eccentricity::new(0.9).unwrap();
446        assert!(low < high);
447    }
448
449    // --- Zero and special float tests ---
450
451    #[test]
452    fn unit_zero_operations() {
453        let z = Km(0.0);
454        assert_eq!((z + z).value(), 0.0);
455        assert_eq!((z - z).value(), 0.0);
456        assert_eq!((z * 100.0).value(), 0.0);
457        assert_eq!((-z).value(), 0.0);
458    }
459
460    #[test]
461    fn unit_division_by_itself() {
462        assert_eq!(Km(42.0) / Km(42.0), 1.0);
463        assert_eq!(KmPerSec(7.784) / KmPerSec(7.784), 1.0);
464        assert_eq!(Seconds(3600.0) / Seconds(3600.0), 1.0);
465        assert_eq!(Mu(3.986e5) / Mu(3.986e5), 1.0);
466        assert_eq!(Radians(PI) / Radians(PI), 1.0);
467    }
468}
solar_line_core/units.rs

solar_line_core/
units.rs
