astroceleste_engine/ephemeris/
kernels.rs

1//! A set of JPL kernels in preference order, and barycentric states of bodies.
2
3use crate::constants::{AU_KM, DAY_S, T0};
4use crate::frames::Vec3;
5use crate::time::Time;
6
7use super::spk::{Segment, Spk, SpkError};
8
9/// NAIF code of the Solar System Barycenter.
10pub const SSB: i32 = 0;
11
12/// One loaded kernel and the Julian-date span that *all* its segments cover (a chart
13/// needs every body, so the usable span is their intersection).
14pub struct Kernel {
15    /// Label used in error messages and diagnostics, e.g. "de440s.bsp".
16    pub name: String,
17    spk: Spk,
18    /// First Julian date (TDB) every segment covers.
19    pub start_jd: f64,
20    /// Last Julian date (TDB) every segment covers.
21    pub end_jd: f64,
22}
23
24impl Kernel {
25    /// Wrap a loaded SPK file; fails if it has no segments.
26    pub fn new(name: impl Into<String>, spk: Spk) -> Result<Self, SpkError> {
27        let segments = spk.segments();
28        if segments.is_empty() {
29            return Err(SpkError::Format("kernel has no segments".into()));
30        }
31        let start_jd = segments
32            .iter()
33            .map(|s| s.start_et / DAY_S + T0)
34            .fold(f64::NEG_INFINITY, f64::max);
35        let end_jd = segments
36            .iter()
37            .map(|s| s.end_et / DAY_S + T0)
38            .fold(f64::INFINITY, f64::min);
39        Ok(Kernel {
40            name: name.into(),
41            spk,
42            start_jd,
43            end_jd,
44        })
45    }
46
47    /// Whether every segment covers Julian date `jd`.
48    pub fn covers(&self, jd: f64) -> bool {
49        self.start_jd <= jd && jd <= self.end_jd
50    }
51
52    /// Whether the kernel has a segment for this body (Skyfield `code in ephemeris`).
53    pub fn contains(&self, code: i32) -> bool {
54        self.spk
55            .segments()
56            .iter()
57            .any(|s| s.target == code || s.center == code)
58    }
59
60    fn segment_for(&self, target: i32) -> Option<&Segment> {
61        self.spk.segments().iter().find(|s| s.target == target)
62    }
63
64    /// Position (au) and velocity (au/day) of `code` relative to the Solar System
65    /// Barycenter, summing the chain of segments outward from the SSB like Skyfield's
66    /// `VectorSum`.
67    #[doc(hidden)] // takes the internal `Time`
68    pub fn barycentric(&self, code: i32, t: &Time) -> Result<(Vec3, Vec3), SpkError> {
69        let mut chain = Vec::new();
70        let mut current = code;
71        while current != SSB {
72            let segment = self.segment_for(current).ok_or(SpkError::NoSegment {
73                target: current,
74                center: SSB,
75            })?;
76            chain.push(segment);
77            current = segment.center;
78            if chain.len() > 8 {
79                return Err(SpkError::Format(format!("segment chain for {code} loops")));
80            }
81        }
82        let mut position = [0.0; 3];
83        let mut velocity = [0.0; 3];
84        for segment in chain.iter().rev() {
85            let (p, v) = self
86                .spk
87                .segment_state_split(segment, t.whole, t.tdb_fraction)?;
88            for k in 0..3 {
89                position[k] += p[k] / AU_KM;
90                velocity[k] += v[k] * DAY_S / AU_KM;
91            }
92        }
93        Ok((position, velocity))
94    }
95}
96
97/// Kernels in preference order: a date is computed with the first kernel covering it.
98#[derive(Default)]
99pub struct KernelSet {
100    kernels: Vec<Kernel>,
101}
102
103impl KernelSet {
104    /// An empty set.
105    pub fn new() -> Self {
106        Self::default()
107    }
108
109    /// Add a kernel, after (so less preferred than) the ones already loaded.
110    pub fn push(&mut self, kernel: Kernel) {
111        self.kernels.push(kernel);
112    }
113
114    /// The kernels, in preference order.
115    pub fn kernels(&self) -> &[Kernel] {
116        &self.kernels
117    }
118
119    /// Whether no kernel is loaded.
120    pub fn is_empty(&self) -> bool {
121        self.kernels.is_empty()
122    }
123
124    /// The preferred kernel covering this Julian date.
125    pub fn for_jd(&self, jd: f64) -> Option<&Kernel> {
126        self.kernels.iter().find(|k| k.covers(jd))
127    }
128
129    /// Combined span of every loaded kernel, as Julian dates.
130    pub fn coverage(&self) -> Option<(f64, f64)> {
131        if self.kernels.is_empty() {
132            return None;
133        }
134        let start = self
135            .kernels
136            .iter()
137            .map(|k| k.start_jd)
138            .fold(f64::INFINITY, f64::min);
139        let end = self
140            .kernels
141            .iter()
142            .map(|k| k.end_jd)
143            .fold(f64::NEG_INFINITY, f64::max);
144        Some((start, end))
145    }
146}
astroceleste_engine/ephemeris/kernels.rs

astroceleste_engine/ephemeris/
kernels.rs
