Skip to main content

pleiades_data/
backend.rs

1use std::sync::Arc;
2
3#[cfg(feature = "packaged-artifact-path")]
4use std::path::Path;
5
6use pleiades_apparent::{
7    precess_ecliptic_date_to_j2000, precess_ecliptic_vector_j2000_to_date, ApparentPlaceError,
8};
9use pleiades_apsides::{
10    apsides, elements_from_state, mean_lunar_elements_of_date, points_from_elements,
11    MOON_MEAN_SEMI_MAJOR_AU, MU_EARTH_MOON_AU3_PER_DAY2,
12};
13use pleiades_backend::{
14    validate_observer_policy, validate_request_policy, validate_zodiac_policy, AccuracyClass,
15    BackendCapabilities, BackendFamily, BackendId, BackendMetadata, BackendProvenance,
16    CelestialBody, CoordinateFrame, EclipticCoordinates, EphemerisBackend, EphemerisError,
17    EphemerisErrorKind, EphemerisRequest, EphemerisResult, Instant, JulianDay, Latitude, Longitude,
18    Motion, QualityAnnotation, TimeScale, ZodiacMode,
19};
20use pleiades_compression::{
21    spherical_state_to_cartesian, CartesianState, CompressedArtifact, SphericalState,
22};
23
24use crate::coverage::packaged_body_coverage_summary_details;
25#[cfg(feature = "packaged-artifact-path")]
26use crate::data::packaged_artifact_from_path;
27use crate::data::{packaged_artifact, packaged_artifact_from_bytes, PackagedArtifactLoadError};
28use crate::lookup::{
29    packaged_artifact_access_summary, packaged_artifact_storage_summary,
30    packaged_frame_treatment_summary, packaged_request_policy_summary,
31};
32use crate::regenerate::{artifact_time_range, map_artifact_error, normalize_lookup_instant};
33use crate::PACKAGE_NAME;
34
35/// A packaged compressed-data backend.
36#[derive(Debug, Clone)]
37pub struct PackagedDataBackend {
38    artifact: Arc<CompressedArtifact>,
39}
40
41impl Default for PackagedDataBackend {
42    fn default() -> Self {
43        Self::new()
44    }
45}
46
47impl PackagedDataBackend {
48    /// Creates a new packaged-data backend backed by the checked-in fixture.
49    pub fn new() -> Self {
50        Self::from_artifact(packaged_artifact().clone())
51    }
52
53    /// Creates a packaged-data backend from an explicit artifact.
54    pub fn from_artifact(artifact: CompressedArtifact) -> Self {
55        Self {
56            artifact: Arc::new(artifact),
57        }
58    }
59
60    /// Creates a packaged-data backend from decoded artifact bytes.
61    ///
62    /// See [`crate::packaged_backend_from_bytes`] for an end-to-end example that
63    /// encodes the checked-in artifact and reloads it through this constructor.
64    pub fn from_bytes(bytes: &[u8]) -> Result<Self, PackagedArtifactLoadError> {
65        Ok(Self::from_artifact(
66            packaged_artifact_from_bytes(bytes).map_err(PackagedArtifactLoadError::Decode)?,
67        ))
68    }
69
70    #[cfg(feature = "packaged-artifact-path")]
71    /// Creates a packaged-data backend from an artifact file.
72    ///
73    /// See [`crate::packaged_backend_from_path`] for an end-to-end example that writes
74    /// the checked-in artifact to a temporary file and reloads it through this
75    /// constructor.
76    pub fn from_path(path: impl AsRef<Path>) -> Result<Self, PackagedArtifactLoadError> {
77        Ok(Self::from_artifact(packaged_artifact_from_path(path)?))
78    }
79
80    fn artifact(&self) -> &CompressedArtifact {
81        &self.artifact
82    }
83
84    /// Assembles a derived lunar point (osculating or mean apsis or node) into a
85    /// backend result: J2000 ecliptic from `eval`, mean-obliquity equatorial,
86    /// extrapolated central-difference motion, `Interpolated` quality.
87    fn derived_point_position(
88        &self,
89        req: &EphemerisRequest,
90        eval: &dyn Fn(Instant) -> Result<EclipticCoordinates, EphemerisError>,
91    ) -> Result<EphemerisResult, EphemerisError> {
92        let ecliptic = eval(req.instant)?;
93        let equatorial = ecliptic.to_equatorial(req.instant.mean_obliquity());
94        let motion = self.derived_point_motion(req.instant, eval)?;
95
96        let mut result = EphemerisResult::new(
97            BackendId::new(PACKAGE_NAME),
98            req.body.clone(),
99            req.instant,
100            req.frame,
101            req.zodiac_mode.clone(),
102            req.apparent,
103        );
104        result.ecliptic = Some(ecliptic);
105        result.equatorial = Some(equatorial);
106        result.motion = Some(motion);
107        result.quality = QualityAnnotation::Interpolated;
108        Ok(result)
109    }
110
111    /// The packaged Moon's geocentric J2000 Cartesian state at `instant`: the
112    /// shared input of every derived lunar point.
113    fn moon_state_j2000(&self, instant: Instant) -> Result<CartesianState, EphemerisError> {
114        let li = normalize_lookup_instant(instant);
115        let ecl = self
116            .artifact
117            .lookup_ecliptic(&CelestialBody::Moon, li)
118            .map_err(map_artifact_error)?;
119        let mot = self
120            .artifact
121            .lookup_motion(&CelestialBody::Moon, li)
122            .map_err(map_artifact_error)?;
123        let dist = ecl.distance_au.ok_or_else(|| {
124            EphemerisError::new(
125                EphemerisErrorKind::InvalidRequest,
126                "packaged Moon lacks distance for a derived lunar point",
127            )
128        })?;
129        Ok(spherical_state_to_cartesian(SphericalState {
130            lon_rad: ecl.longitude.degrees().to_radians(),
131            lat_rad: ecl.latitude.degrees().to_radians(),
132            dist_au: dist,
133            lon_rate_rad_per_day: mot.longitude_deg_per_day.unwrap_or(0.0).to_radians(),
134            lat_rate_rad_per_day: mot.latitude_deg_per_day.unwrap_or(0.0).to_radians(),
135            dist_rate_au_per_day: mot.distance_au_per_day.unwrap_or(0.0),
136        }))
137    }
138
139    fn osculating_apsis_ecliptic(
140        &self,
141        body: &CelestialBody,
142        instant: Instant,
143    ) -> Result<EclipticCoordinates, EphemerisError> {
144        let cart = self.moon_state_j2000(instant)?;
145        let aps = apsides(cart.pos_au, cart.vel_au_per_day, MU_EARTH_MOON_AU3_PER_DAY2).map_err(
146            |_| {
147                EphemerisError::new(
148                    EphemerisErrorKind::InvalidRequest,
149                    "osculating apsis undefined for the lunar state at this instant",
150                )
151            },
152        )?;
153        let point = match body {
154            CelestialBody::TrueApogee => aps.apogee,
155            CelestialBody::TruePerigee => aps.perigee,
156            _ => {
157                return Err(EphemerisError::new(
158                    EphemerisErrorKind::InvalidRequest,
159                    "not an osculating-apsis body",
160                ))
161            }
162        };
163        Ok(EclipticCoordinates::new(
164            Longitude::from_degrees(point.longitude_deg),
165            Latitude::from_degrees(point.latitude_deg),
166            Some(point.distance_au),
167        ))
168    }
169
170    /// Osculating ascending node of the geocentric lunar orbit, J2000 boundary
171    /// frame.
172    ///
173    /// The node is the intersection of the orbit plane with the *reference*
174    /// plane, so it must be formed in the plane consumers will read it in: the
175    /// mean ecliptic of date (spec §1; forming it in J2000 and rotating the
176    /// point would misplace it by ≈ tilt/sin(i) ≈ 0.04°). Position and velocity
177    /// are rotated J2000 → mean-of-date, the ellipse is formed there, and the
178    /// node point is precessed back to J2000 like the ELP point channels
179    /// (issue #57). The chart layer's forward precession + Δψ then reproduces
180    /// Swiss Ephemeris `SE_TRUE_NODE`.
181    fn osculating_node_ecliptic(
182        &self,
183        instant: Instant,
184    ) -> Result<EclipticCoordinates, EphemerisError> {
185        let jd_tt = instant.julian_day.days();
186        let cart = self.moon_state_j2000(instant)?;
187        let pos = rotate_j2000_to_mean_of_date(cart.pos_au, jd_tt)?;
188        let vel = rotate_j2000_to_mean_of_date(cart.vel_au_per_day, jd_tt)?;
189        let undefined = |_| {
190            EphemerisError::new(
191                EphemerisErrorKind::InvalidRequest,
192                "osculating node undefined for the lunar state at this instant",
193            )
194        };
195        let elements =
196            elements_from_state(pos, vel, MU_EARTH_MOON_AU3_PER_DAY2).map_err(undefined)?;
197        let node = points_from_elements(&elements, false)
198            .map_err(undefined)?
199            .ascending;
200        let j2000 = precess_ecliptic_date_to_j2000(node.longitude_deg, node.latitude_deg, jd_tt)
201            .map_err(map_precession_error)?;
202        Ok(EclipticCoordinates::new(
203            Longitude::from_degrees(j2000.longitude_deg),
204            Latitude::from_degrees(j2000.latitude_deg),
205            Some(node.distance_au),
206        ))
207    }
208
209    /// A mean lunar point (mean node, mean apogee or mean perigee), J2000
210    /// boundary frame.
211    ///
212    /// The point is placed on the Moon's mean orbit in the mean ecliptic of
213    /// date — the apsides through the inclined orbit, as Swiss Ephemeris'
214    /// `SE_MEAN_APOG` is, not as the raw longitude-of-perigee element — and
215    /// precessed back to J2000 like the osculating node. The elements are
216    /// analytic, but the point is served only inside the packaged window so
217    /// the backend's advertised range holds for every body it answers.
218    fn mean_lunar_point_ecliptic(
219        &self,
220        body: &CelestialBody,
221        instant: Instant,
222    ) -> Result<EclipticCoordinates, EphemerisError> {
223        // Window probe: the Moon series spans the packaged window.
224        self.artifact
225            .lookup_ecliptic(&CelestialBody::Moon, normalize_lookup_instant(instant))
226            .map_err(map_artifact_error)?;
227
228        let jd_tt = instant.julian_day.days();
229        let points =
230            points_from_elements(&mean_lunar_elements_of_date(jd_tt), false).map_err(|_| {
231                EphemerisError::new(
232                    EphemerisErrorKind::InvalidRequest,
233                    "mean lunar point undefined at this instant",
234                )
235            })?;
236        let (point, distance_au) = match body {
237            // Swiss Ephemeris reports the mean lunar distance for SE_MEAN_NODE,
238            // not the orbit radius at the node.
239            CelestialBody::MeanNode => (points.ascending, MOON_MEAN_SEMI_MAJOR_AU),
240            CelestialBody::MeanApogee => (points.aphelion, points.aphelion.distance_au),
241            CelestialBody::MeanPerigee => (points.perihelion, points.perihelion.distance_au),
242            _ => {
243                return Err(EphemerisError::new(
244                    EphemerisErrorKind::InvalidRequest,
245                    "not a mean lunar point",
246                ))
247            }
248        };
249        let j2000 = precess_ecliptic_date_to_j2000(point.longitude_deg, point.latitude_deg, jd_tt)
250            .map_err(map_precession_error)?;
251        Ok(EclipticCoordinates::new(
252            Longitude::from_degrees(j2000.longitude_deg),
253            Latitude::from_degrees(j2000.latitude_deg),
254            Some(distance_au),
255        ))
256    }
257
258    /// Motion of a derived lunar point: central differences over ±0.5 and
259    /// ±0.25 day, Richardson-extrapolated to cancel their leading truncation
260    /// term, `(4·D(0.25) − D(0.5)) / 3`.
261    ///
262    /// The osculating points oscillate within a fortnight, so a plain ±0.5-day
263    /// difference reads their speed low wherever it peaks: the true node by
264    /// about 3″/day in its short direct spells, enough to hide them (#108).
265    /// A shorter plain difference is no better, because it amplifies the small
266    /// steps in the packaged Moon's velocity between fitted segments. Measured
267    /// against the derivative of the node and apogee formed from the exact
268    /// DE440 Moon, 2017–2027 every 0.05 day: true node max 0.71″/day (rms
269    /// 0.07″) against 3.9″/day for the ±0.5-day difference; true apogee max
270    /// 13″/day (rms 2.3″) against 161″/day.
271    ///
272    /// A probe that falls outside the packaged window degrades to `None`
273    /// channels rather than failing the position. The outer probes sit at
274    /// ±0.5 day as before, so the same instants degrade.
275    fn derived_point_motion(
276        &self,
277        instant: Instant,
278        eval: &dyn Fn(Instant) -> Result<EclipticCoordinates, EphemerisError>,
279    ) -> Result<Motion, EphemerisError> {
280        const OUTER_HALF_SPAN_DAYS: f64 = 0.5;
281        const INNER_HALF_SPAN_DAYS: f64 = 0.25;
282        let probe = |days: f64| -> Result<Option<EclipticCoordinates>, EphemerisError> {
283            let shifted = Instant::new(
284                JulianDay::from_days(instant.julian_day.days() + days),
285                instant.scale,
286            );
287            match eval(shifted) {
288                Ok(e) => Ok(Some(e)),
289                Err(ref e) if e.kind == EphemerisErrorKind::OutOfRangeInstant => Ok(None),
290                Err(e) => Err(e),
291            }
292        };
293        let (Some(outer_before), Some(inner_before), Some(inner_after), Some(outer_after)) = (
294            probe(-OUTER_HALF_SPAN_DAYS)?,
295            probe(-INNER_HALF_SPAN_DAYS)?,
296            probe(INNER_HALF_SPAN_DAYS)?,
297            probe(OUTER_HALF_SPAN_DAYS)?,
298        ) else {
299            return Ok(Motion::new(None, None, None));
300        };
301
302        // One rate per span, then the extrapolation.
303        let extrapolate = |outer: f64, inner: f64| (4.0 * inner - outer) / 3.0;
304        let rate = |before: f64, after: f64, half_span: f64| (after - before) / (2.0 * half_span);
305        let lon_rate = |before: &EclipticCoordinates, after: &EclipticCoordinates, half_span| {
306            let mut dlon = after.longitude.degrees() - before.longitude.degrees();
307            while dlon > 180.0 {
308                dlon -= 360.0;
309            }
310            while dlon < -180.0 {
311                dlon += 360.0;
312            }
313            dlon / (2.0 * half_span)
314        };
315        let lat_rate = |before: &EclipticCoordinates, after: &EclipticCoordinates, half_span| {
316            rate(
317                before.latitude.degrees(),
318                after.latitude.degrees(),
319                half_span,
320            )
321        };
322        let dist_rate =
323            |before: &EclipticCoordinates, after: &EclipticCoordinates, half_span| match (
324                before.distance_au,
325                after.distance_au,
326            ) {
327                (Some(b), Some(a)) => Some(rate(b, a, half_span)),
328                _ => None,
329            };
330
331        let dlon_per_day = extrapolate(
332            lon_rate(&outer_before, &outer_after, OUTER_HALF_SPAN_DAYS),
333            lon_rate(&inner_before, &inner_after, INNER_HALF_SPAN_DAYS),
334        );
335        let dlat_per_day = extrapolate(
336            lat_rate(&outer_before, &outer_after, OUTER_HALF_SPAN_DAYS),
337            lat_rate(&inner_before, &inner_after, INNER_HALF_SPAN_DAYS),
338        );
339        let ddist_per_day = match (
340            dist_rate(&outer_before, &outer_after, OUTER_HALF_SPAN_DAYS),
341            dist_rate(&inner_before, &inner_after, INNER_HALF_SPAN_DAYS),
342        ) {
343            (Some(outer), Some(inner)) => Some(extrapolate(outer, inner)),
344            _ => None,
345        };
346        Ok(Motion::new(
347            Some(dlon_per_day),
348            Some(dlat_per_day),
349            ddist_per_day,
350        ))
351    }
352}
353
354/// Rotates a J2000 mean-ecliptic vector into the mean ecliptic of date at
355/// `jd_tt`, preserving its magnitude. The same map applies to position and
356/// velocity vectors alike; the ELP backend's osculating node shares it.
357fn rotate_j2000_to_mean_of_date(v: [f64; 3], jd_tt: f64) -> Result<[f64; 3], EphemerisError> {
358    precess_ecliptic_vector_j2000_to_date(v, jd_tt).map_err(map_precession_error)
359}
360
361fn map_precession_error(e: ApparentPlaceError) -> EphemerisError {
362    EphemerisError::new(
363        EphemerisErrorKind::InvalidRequest,
364        format!("precession failed for derived lunar point: {e}"),
365    )
366}
367
368impl EphemerisBackend for PackagedDataBackend {
369    fn metadata(&self) -> BackendMetadata {
370        let artifact = self.artifact();
371        let bodies = artifact
372            .bodies
373            .iter()
374            .map(|series| series.body.clone())
375            .filter(|body| !crate::is_carried_but_unserved(body))
376            .collect::<Vec<_>>();
377        let range = artifact_time_range(artifact);
378
379        BackendMetadata {
380            id: BackendId::new(PACKAGE_NAME),
381            version: format!(
382                "{} checksum:{:016x}",
383                artifact.header.version, artifact.checksum
384            ),
385            family: BackendFamily::CompressedData,
386            provenance: BackendProvenance {
387                summary: artifact.header.source.clone(),
388                data_sources: vec![
389                    packaged_body_coverage_summary_details()
390                        .validated_summary_line()
391                        .unwrap_or_else(|error| {
392                            format!("Packaged body set: unavailable ({error})")
393                        }),
394                    packaged_request_policy_summary().to_string(),
395                    packaged_frame_treatment_summary().to_string(),
396                    packaged_artifact_storage_summary().to_string(),
397                    packaged_artifact_access_summary().to_string(),
398                ],
399            },
400            nominal_range: range,
401            supported_time_scales: vec![TimeScale::Tt, TimeScale::Tdb],
402            body_claims: {
403                let declared = crate::packaged_body_claims();
404                let mut claims: Vec<_> = bodies
405                    .iter()
406                    .map(|body| {
407                        declared
408                            .iter()
409                            .find(|c| &c.body == body)
410                            .cloned()
411                            .unwrap_or_else(|| pleiades_backend::BodyClaim::from(body.clone()))
412                    })
413                    .collect();
414                claims.extend(crate::apsis_body_claims());
415                claims.extend(crate::true_node_body_claims());
416                claims.extend(crate::mean_lunar_point_body_claims());
417                claims
418            },
419            supported_frames: vec![CoordinateFrame::Ecliptic, CoordinateFrame::Equatorial],
420            capabilities: BackendCapabilities {
421                geocentric: true,
422                topocentric: false,
423                apparent: false,
424                mean: true,
425                batch: true,
426                native_sidereal: false,
427            },
428            accuracy: AccuracyClass::Approximate,
429            deterministic: true,
430            offline: true,
431        }
432    }
433
434    fn supports_body(&self, body: CelestialBody) -> bool {
435        matches!(
436            body,
437            CelestialBody::TrueApogee
438                | CelestialBody::TruePerigee
439                | CelestialBody::TrueNode
440                | CelestialBody::MeanNode
441                | CelestialBody::MeanApogee
442                | CelestialBody::MeanPerigee
443        ) || (!crate::is_carried_but_unserved(&body)
444            && self
445                .artifact
446                .bodies
447                .iter()
448                .any(|series| series.body == body))
449    }
450
451    fn position(&self, req: &EphemerisRequest) -> Result<EphemerisResult, EphemerisError> {
452        if !matches!(req.instant.scale, TimeScale::Tt | TimeScale::Tdb) {
453            return Err(EphemerisError::new(
454                EphemerisErrorKind::UnsupportedTimeScale,
455                "packaged data only supports TT or TDB requests",
456            ));
457        }
458
459        validate_request_policy(
460            req,
461            "packaged data",
462            &[TimeScale::Tt, TimeScale::Tdb],
463            &[CoordinateFrame::Ecliptic, CoordinateFrame::Equatorial],
464            true,
465            false,
466        )?;
467
468        validate_zodiac_policy(req, "packaged data", &[ZodiacMode::Tropical])?;
469
470        validate_observer_policy(req, "packaged data", false)?;
471
472        if matches!(
473            req.body,
474            CelestialBody::TrueApogee | CelestialBody::TruePerigee
475        ) {
476            let body = req.body.clone();
477            return self.derived_point_position(req, &|i| self.osculating_apsis_ecliptic(&body, i));
478        }
479        if req.body == CelestialBody::TrueNode {
480            return self.derived_point_position(req, &|i| self.osculating_node_ecliptic(i));
481        }
482        if matches!(
483            req.body,
484            CelestialBody::MeanNode | CelestialBody::MeanApogee | CelestialBody::MeanPerigee
485        ) {
486            let body = req.body.clone();
487            return self.derived_point_position(req, &|i| self.mean_lunar_point_ecliptic(&body, i));
488        }
489
490        if crate::is_carried_but_unserved(&req.body) {
491            return Err(EphemerisError::new(
492                EphemerisErrorKind::UnsupportedBody,
493                format!(
494                    "packaged data carries {} but does not serve it: its segments are fitted \
495                     to rows too sparse to interpolate. Serve it from pleiades_jpl::SpkBackend \
496                     with a JPL kernel",
497                    req.body
498                ),
499            ));
500        }
501
502        let lookup_instant = normalize_lookup_instant(req.instant);
503        let ecliptic = self
504            .artifact
505            .lookup_ecliptic(&req.body, lookup_instant)
506            .map_err(map_artifact_error)?;
507        let equatorial = ecliptic.to_equatorial(req.instant.mean_obliquity());
508        let motion = self
509            .artifact
510            .lookup_motion(&req.body, lookup_instant)
511            .map_err(map_artifact_error)?;
512
513        let mut result = EphemerisResult::new(
514            BackendId::new(PACKAGE_NAME),
515            req.body.clone(),
516            req.instant,
517            req.frame,
518            req.zodiac_mode.clone(),
519            req.apparent,
520        );
521        result.ecliptic = Some(ecliptic);
522        result.equatorial = Some(equatorial);
523        result.motion = Some(motion);
524        result.quality = QualityAnnotation::Interpolated;
525        Ok(result)
526    }
527}
528
529#[cfg(test)]
530mod coupling_fixture_tests {
531    use super::*;
532
533    /// Golden fixture pinning `PackagedDataBackend::metadata()`'s
534    /// `provenance.data_sources` to its pre-Slice-C rendered values.
535    ///
536    /// This test exists to prove byte-identity when the five report-prose
537    /// renderers backing `data_sources` are replaced with an inline rebuild
538    /// from retained structured accessors (Slice C, Task 2). It must keep
539    /// passing after that rebuild — if it fails, the rebuild drifted; fix
540    /// the rebuild, never this fixture.
541    #[test]
542    fn backend_metadata_data_sources_is_stable() {
543        use pleiades_backend::EphemerisBackend;
544        let metadata = PackagedDataBackend::default().metadata();
545        let expected: &[&str] = &[
546            "Packaged body set: 11 bundled bodies (Sun, Moon, Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, asteroid:433-Eros)",
547            "Packaged request policy: geocentric-only; frames=Ecliptic, Equatorial; time scales=TT, TDB; zodiac modes=Tropical; apparentness=Mean; topocentric observer=false; lookup epoch policy=TT-grid retag without relativistic correction; TDB lookup epochs are re-tagged onto the TT grid without applying a relativistic correction",
548            "checked-in compressed artifact stores J2000 ecliptic coordinates directly; equatorial coordinates are reconstructed from the stored channels and mean-obliquity transform",
549            "Quantized linear segments stored in pleiades-compression artifact format; body-indexed segment tables support random access by body and lookup time across the advertised range; ecliptic and equatorial coordinates are reconstructed at runtime from stored channels; apparent, topocentric, and sidereal outputs remain unsupported; motion/speed is derived from fitted segment derivatives",
550            "packaged artifact access: checked-in fixture only; explicit artifact-path loading disabled",
551        ];
552        assert_eq!(
553            metadata.provenance.data_sources, expected,
554            "backend metadata data_sources drifted:\n{:#?}",
555            metadata.provenance.data_sources
556        );
557    }
558}