spicekit 0.1.0

Pure-Rust reader for NASA/NAIF SPICE kernel formats (DAF, SPK, PCK, text kernels). Independent of the CSpice toolkit.
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
//! Binary PCK (Planetary Constants Kernel) reader.
//!
//! Adam-core only needs time-varying planetary body orientations for the
//! Earth-associated `ITRF93` frame. Binary PCK shares the DAF container
//! with SPK, but the summary layout differs (NI=5 instead of 6, with
//! `center` absent: segments describe a single body-fixed frame wrt a
//! reference inertial frame). The payload is the standard Type 2
//! Chebyshev record, but with 3 Euler-angle channels rather than 3
//! position channels. The angles follow the NAIF 3-1-3 convention:
//! `(t1, t2, t3)` such that `R_ref→body = Rz(t3) · Rx(t2) · Rz(t1)`.
//! The `ref` frame is the segment's reference inertial frame — the
//! standard Earth PCKs use `ECLIPJ2000` (NAIF frame ID 17), so
//! consumers must compose with the J2000 ↔ ECLIPJ2000 static rotation
//! when producing `sxform("J2000","ITRF93",et)`.
//!
//! The evaluation output is `[t1, t2, t3, dt1, dt2, dt3]` at a
//! requested TDB ET. Frame consumers assemble a 3x3 rotation (and its
//! time derivative) from these angles via
//! [`crate::frame::pck_euler_rotation_and_derivative`].
//!
//! Kernel precedence: three PCKs (`predict`, `historical`, `high_prec`)
//! may cover overlapping epoch ranges. Last-loaded-wins matches
//! CSPICE `furnsh` semantics: the caller is expected to load kernels in
//! precedence order (least- to most-precise) and the file they opened
//! most recently "wins" for overlapping coverage. This crate offers a
//! per-file reader; multi-file composition is handled at the Python
//! layer with a simple list-of-readers dispatcher.

use std::path::Path;

use thiserror::Error;

use crate::daf::{DafError, DafFile};
use crate::spk::cheby_val_and_deriv;

#[derive(Debug, Error)]
pub enum PckError {
    #[error(transparent)]
    Daf(#[from] DafError),
    #[error("no segment covers body {body} in frame {ref_frame} at et {et}")]
    NoCoverage { body: i32, ref_frame: i32, et: f64 },
    #[error("unsupported PCK data type {0}")]
    UnsupportedType(i32),
    #[error("malformed PCK Type 2 segment: {0}")]
    BadType2(&'static str),
}

#[derive(Clone)]
pub struct PckSegment {
    /// Body-fixed frame ID described by this segment (e.g. 3000 for ITRF93).
    pub body_frame: i32,
    /// Inertial reference frame ID (always 1 = J2000 in practice).
    pub ref_frame: i32,
    pub data_type: i32,
    pub start_et: f64,
    pub end_et: f64,
    pub name: String,
    payload: PckPayload,
}

#[derive(Clone)]
enum PckPayload {
    Type2(PckType2),
    Unsupported,
}

/// PCK Type 2: 3-channel Chebyshev (RA, DEC, W) with analytic derivative.
///
/// Record layout mirrors SPK Type 2 but the three angle channels
/// replace position (x, y, z). Trailer `[INIT, INTLEN, RSIZE, N]` is
/// identical.
#[derive(Clone)]
struct PckType2 {
    file: DafFile,
    init: f64,
    intlen: f64,
    rsize: usize,
    n_records: usize,
    n_coef: usize,
    start_addr: u32,
}

impl PckType2 {
    fn from_segment(file: &DafFile, start_addr: u32, end_addr: u32) -> Result<Self, PckError> {
        let trailer = file.read_doubles(end_addr - 3, end_addr)?;
        let init = trailer[0];
        let intlen = trailer[1];
        let rsize = trailer[2] as usize;
        let n_records = trailer[3] as usize;
        if rsize < 2 || (rsize - 2) % 3 != 0 {
            return Err(PckError::BadType2("RSIZE not 2 + 3N"));
        }
        let n_coef = (rsize - 2) / 3;
        if n_coef == 0 || intlen <= 0.0 {
            return Err(PckError::BadType2("degree<0 or INTLEN<=0"));
        }
        Ok(PckType2 {
            file: file.clone(),
            init,
            intlen,
            rsize,
            n_records,
            n_coef,
            start_addr,
        })
    }

    fn evaluate(&self, et: f64) -> Result<[f64; 6], PckError> {
        let raw_idx = ((et - self.init) / self.intlen).floor() as isize;
        let idx = raw_idx.clamp(0, self.n_records as isize - 1) as usize;
        let rec_start = self.start_addr + (idx * self.rsize) as u32;
        let rec_end = rec_start + self.rsize as u32 - 1;
        let rec = self.file.read_doubles(rec_start, rec_end)?;

        let mid = rec[0];
        let radius = rec[1];
        if radius == 0.0 {
            return Err(PckError::BadType2("RADIUS == 0"));
        }
        let s = (et - mid) / radius;

        let n = self.n_coef;
        let ra_c = &rec[2..2 + n];
        let dec_c = &rec[2 + n..2 + 2 * n];
        let w_c = &rec[2 + 2 * n..2 + 3 * n];

        let (ra, dra) = cheby_val_and_deriv(ra_c, s);
        let (dec, ddec) = cheby_val_and_deriv(dec_c, s);
        let (w, dw) = cheby_val_and_deriv(w_c, s);

        let inv_r = 1.0 / radius;
        Ok([ra, dec, w, dra * inv_r, ddec * inv_r, dw * inv_r])
    }
}

pub struct PckFile {
    segments: Vec<PckSegment>,
}

impl PckFile {
    pub fn open<P: AsRef<Path>>(path: P) -> Result<Self, PckError> {
        let daf = DafFile::open(path)?;
        Self::from_daf(daf)
    }

    pub fn from_daf(daf: DafFile) -> Result<Self, PckError> {
        let mut segments = Vec::new();
        for summary in daf.summaries()? {
            if summary.doubles.len() < 2 || summary.integers.len() < 5 {
                continue;
            }
            let start_et = summary.doubles[0];
            let end_et = summary.doubles[1];
            let body_frame = summary.integers[0];
            let ref_frame = summary.integers[1];
            let data_type = summary.integers[2];
            let start_addr = summary.integers[3] as u32;
            let end_addr = summary.integers[4] as u32;

            let payload = match data_type {
                2 => PckPayload::Type2(PckType2::from_segment(&daf, start_addr, end_addr)?),
                _ => PckPayload::Unsupported,
            };

            segments.push(PckSegment {
                body_frame,
                ref_frame,
                data_type,
                start_et,
                end_et,
                name: summary.name,
                payload,
            });
        }
        Ok(PckFile { segments })
    }

    pub fn segments(&self) -> &[PckSegment] {
        &self.segments
    }

    /// Raw evaluation: returns `[RA, DEC, W, dRA/dt, dDEC/dt, dW/dt]` at
    /// ET for the requested body-fixed frame wrt its reference frame.
    /// Returns `NoCoverage` if no segment covers the epoch.
    pub fn euler_state(&self, body_frame: i32, et: f64) -> Result<[f64; 6], PckError> {
        Ok(self.euler_state_with_ref(body_frame, et)?.1)
    }

    /// Raw evaluation returning `(ref_frame_id, [RA, DEC, W, dRA, dDEC, dW])`.
    /// The reference frame is usually 1 (J2000) or 17 (ECLIPJ2000); the
    /// caller composes with the inter-inertial rotation to land in a
    /// chosen inertial frame.
    pub fn euler_state_with_ref(
        &self,
        body_frame: i32,
        et: f64,
    ) -> Result<(i32, [f64; 6]), PckError> {
        // Scan newest-first so late segments (typically higher precision)
        // win when segments overlap within a single file.
        for seg in self.segments.iter().rev() {
            if seg.body_frame != body_frame {
                continue;
            }
            if et < seg.start_et || et > seg.end_et {
                continue;
            }
            let state = match &seg.payload {
                PckPayload::Type2(t) => t.evaluate(et)?,
                PckPayload::Unsupported => return Err(PckError::UnsupportedType(seg.data_type)),
            };
            return Ok((seg.ref_frame, state));
        }
        Err(PckError::NoCoverage {
            body: body_frame,
            ref_frame: 1,
            et,
        })
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::daf::{DOUBLE_BYTES, RECORD_BYTES};
    use std::io::Write;
    use tempfile::NamedTempFile;

    /// ITRF93 NAIF frame ID (what Earth PCKs actually record).
    const ITRF93_FRAME: i32 = 3000;
    /// ECLIPJ2000 frame ID (the reference frame the standard Earth PCKs use).
    const ECLIPJ2000_FRAME: i32 = 17;

    /// Build a PCK file with one Type 2 segment. `coefs_ra`, `coefs_dec`,
    /// and `coefs_w` must all have the same length (= n_coef) and will
    /// appear in a single data record spanning `[mid-radius, mid+radius]`.
    // Test helper that lays out a PCK byte-for-byte; each argument is a
    // distinct file field, so a struct would add noise without clarity.
    #[allow(clippy::too_many_arguments)]
    fn build_single_segment_pck(
        body_frame: i32,
        ref_frame: i32,
        start_et: f64,
        end_et: f64,
        mid: f64,
        radius: f64,
        coefs_ra: &[f64],
        coefs_dec: &[f64],
        coefs_w: &[f64],
    ) -> NamedTempFile {
        assert_eq!(coefs_ra.len(), coefs_dec.len());
        assert_eq!(coefs_ra.len(), coefs_w.len());
        let n_coef = coefs_ra.len();
        let rsize = 2 + 3 * n_coef; // [mid, radius, RA..., DEC..., W...]
        let n_records = 1;

        // One data record + 4-double trailer [INIT, INTLEN, RSIZE, N].
        let mut data: Vec<f64> = Vec::with_capacity(rsize + 4);
        data.push(mid);
        data.push(radius);
        data.extend_from_slice(coefs_ra);
        data.extend_from_slice(coefs_dec);
        data.extend_from_slice(coefs_w);
        // trailer
        data.push(start_et); // INIT
        data.push(end_et - start_et); // INTLEN = full coverage for 1 record
        data.push(rsize as f64);
        data.push(n_records as f64);

        // Data starts at record 4 in our layout.
        let data_start_record = 4u32;
        let doubles_per_record = (RECORD_BYTES / DOUBLE_BYTES) as u32;
        let data_start_addr = 1 + (data_start_record - 1) * doubles_per_record;
        let data_end_addr = data_start_addr + data.len() as u32 - 1;

        // Build file: record 1 (file record), record 2 (summary), record 3 (names),
        // record 4 (data).
        let mut record1 = vec![0u8; RECORD_BYTES];
        record1[0..8].copy_from_slice(b"DAF/PCK ");
        record1[8..12].copy_from_slice(&2u32.to_le_bytes()); // nd
        record1[12..16].copy_from_slice(&5u32.to_le_bytes()); // ni
        for b in &mut record1[16..76] {
            *b = b' ';
        }
        record1[76..80].copy_from_slice(&2u32.to_le_bytes()); // fward
        record1[80..84].copy_from_slice(&2u32.to_le_bytes()); // bward
        record1[84..88].copy_from_slice(&0u32.to_le_bytes()); // free
        record1[88..96].copy_from_slice(b"LTL-IEEE");

        let mut record2 = vec![0u8; RECORD_BYTES];
        record2[0..8].copy_from_slice(&0.0f64.to_le_bytes()); // NEXT
        record2[8..16].copy_from_slice(&0.0f64.to_le_bytes()); // PREV
        record2[16..24].copy_from_slice(&1.0f64.to_le_bytes()); // NSUM

        // nd=2 leading doubles (start_et, end_et), ni=5 integers
        // (body_frame, ref_frame, dtype, start_addr, end_addr).
        let soff = 24;
        record2[soff..soff + 8].copy_from_slice(&start_et.to_le_bytes());
        record2[soff + 8..soff + 16].copy_from_slice(&end_et.to_le_bytes());
        let int_start = soff + 2 * DOUBLE_BYTES;
        record2[int_start..int_start + 4].copy_from_slice(&body_frame.to_le_bytes());
        record2[int_start + 4..int_start + 8].copy_from_slice(&ref_frame.to_le_bytes());
        record2[int_start + 8..int_start + 12].copy_from_slice(&2i32.to_le_bytes()); // dtype=2
        record2[int_start + 12..int_start + 16]
            .copy_from_slice(&(data_start_addr as i32).to_le_bytes());
        record2[int_start + 16..int_start + 20]
            .copy_from_slice(&(data_end_addr as i32).to_le_bytes());

        let mut record3 = vec![b' '; RECORD_BYTES];
        let name = b"SYNTHETIC TEST SEGMENT";
        record3[..name.len()].copy_from_slice(name);

        let mut data_bytes = Vec::new();
        for d in &data {
            data_bytes.extend_from_slice(&d.to_le_bytes());
        }
        while data_bytes.len() % RECORD_BYTES != 0 {
            data_bytes.push(0);
        }

        let mut all = Vec::new();
        all.extend_from_slice(&record1);
        all.extend_from_slice(&record2);
        all.extend_from_slice(&record3);
        all.extend_from_slice(&data_bytes);

        let mut tmp = NamedTempFile::new().unwrap();
        tmp.write_all(&all).unwrap();
        tmp.flush().unwrap();
        tmp
    }

    #[test]
    fn from_daf_recovers_segment_metadata() {
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -3e8,
            3e8,
            0.0,
            3e8,
            &[0.1, 0.0],
            &[0.2, 0.0],
            &[0.3, 0.5],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        let segs = pck.segments();
        assert_eq!(segs.len(), 1);
        let seg = &segs[0];
        assert_eq!(seg.body_frame, ITRF93_FRAME);
        assert_eq!(seg.ref_frame, ECLIPJ2000_FRAME);
        assert_eq!(seg.data_type, 2);
        assert_eq!(seg.start_et, -3e8);
        assert_eq!(seg.end_et, 3e8);
        assert_eq!(seg.name, "SYNTHETIC TEST SEGMENT");
    }

    #[test]
    fn euler_state_at_record_midpoint_returns_degree_zero_coefs() {
        // At et = mid, s = 0, so Chebyshev T0(0)=1, T1(0)=0, T2(0)=-1, ...
        // With coefs [c0, c1] => value = c0 - 0*c1 = c0.
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -50.0,
            50.0,
            0.0,
            50.0,
            &[1.25, 0.0],
            &[2.5, 0.0],
            &[3.75, 0.5],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        let state = pck.euler_state(ITRF93_FRAME, 0.0).expect("eval at mid");
        assert!((state[0] - 1.25).abs() < 1e-15);
        assert!((state[1] - 2.5).abs() < 1e-15);
        assert!((state[2] - 3.75).abs() < 1e-15);
        // dW = c1 * dT1/ds * 1/radius = 0.5 * 1 * 1/50 = 0.01
        assert!((state[3] - 0.0).abs() < 1e-15);
        assert!((state[4] - 0.0).abs() < 1e-15);
        assert!((state[5] - 0.01).abs() < 1e-15);
    }

    #[test]
    fn euler_state_with_ref_returns_reference_frame_id() {
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -10.0,
            10.0,
            0.0,
            10.0,
            &[0.0, 0.0],
            &[0.0, 0.0],
            &[0.0, 0.0],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        let (ref_id, _state) = pck.euler_state_with_ref(ITRF93_FRAME, 0.0).expect("eval");
        assert_eq!(ref_id, ECLIPJ2000_FRAME);
    }

    #[test]
    fn derivative_matches_central_finite_difference() {
        // Use non-trivial coefficients; verify dW/dt against symmetric FD.
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -1000.0,
            1000.0,
            0.0,
            1000.0,
            &[0.1, -0.05, 0.02, -0.003],
            &[0.2, 0.04, -0.01, 0.005],
            &[0.3, 0.08, 0.015, -0.004],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        for &et in &[-500.0_f64, -100.0, 0.0, 250.0, 750.0] {
            let h = 1e-3_f64;
            let [_, _, _, dra, ddec, dw] = pck.euler_state(ITRF93_FRAME, et).unwrap();
            let [rap, decp, wp, ..] = pck.euler_state(ITRF93_FRAME, et + h).unwrap();
            let [ram, decm, wm, ..] = pck.euler_state(ITRF93_FRAME, et - h).unwrap();
            let dra_fd = (rap - ram) / (2.0 * h);
            let ddec_fd = (decp - decm) / (2.0 * h);
            let dw_fd = (wp - wm) / (2.0 * h);
            assert!(
                (dra - dra_fd).abs() < 1e-9,
                "dra et={et}: {dra} vs {dra_fd}"
            );
            assert!(
                (ddec - ddec_fd).abs() < 1e-9,
                "ddec et={et}: {ddec} vs {ddec_fd}"
            );
            assert!((dw - dw_fd).abs() < 1e-9, "dw et={et}: {dw} vs {dw_fd}");
        }
    }

    #[test]
    fn no_coverage_when_body_frame_not_present() {
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -10.0,
            10.0,
            0.0,
            10.0,
            &[0.0, 0.0],
            &[0.0, 0.0],
            &[0.0, 0.0],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        let err = pck.euler_state(9999, 0.0).unwrap_err();
        match err {
            PckError::NoCoverage { body, et, .. } => {
                assert_eq!(body, 9999);
                assert_eq!(et, 0.0);
            }
            other => panic!("expected NoCoverage, got {other:?}"),
        }
    }

    #[test]
    fn no_coverage_when_epoch_outside_segment_range() {
        // Segment covers [-10, 10]; query at et=50 must miss.
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -10.0,
            10.0,
            0.0,
            10.0,
            &[0.0, 0.0],
            &[0.0, 0.0],
            &[0.0, 0.0],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        let err = pck.euler_state(ITRF93_FRAME, 50.0).unwrap_err();
        assert!(matches!(err, PckError::NoCoverage { .. }));
        let err_neg = pck.euler_state(ITRF93_FRAME, -100.0).unwrap_err();
        assert!(matches!(err_neg, PckError::NoCoverage { .. }));
    }

    #[test]
    fn segment_value_matches_chebyshev_by_hand() {
        // RA coefs [1, 2, 3] → f(s) = 1 + 2s + 3*(2s² - 1) = 6s² + 2s - 2.
        // At et=25, mid=0, radius=50 → s = 0.5 → f(0.5) = 6*0.25 + 1 - 2 = 0.5.
        // df/ds = 12s + 2 → 8 at s=0.5; df/dt = df/ds / radius = 0.16.
        let tmp = build_single_segment_pck(
            ITRF93_FRAME,
            ECLIPJ2000_FRAME,
            -50.0,
            50.0,
            0.0,
            50.0,
            &[1.0, 2.0, 3.0],
            &[0.0, 0.0, 0.0],
            &[0.0, 0.0, 0.0],
        );
        let pck = PckFile::open(tmp.path()).expect("open");
        let state = pck.euler_state(ITRF93_FRAME, 25.0).expect("eval");
        assert!((state[0] - 0.5).abs() < 1e-14, "RA={}", state[0]);
        assert!((state[3] - 0.16).abs() < 1e-14, "dRA={}", state[3]);
    }

    /// Descriptor for one PCK segment in a multi-segment file.
    struct SegSpec {
        body_frame: i32,
        ref_frame: i32,
        start_et: f64,
        end_et: f64,
        mid: f64,
        radius: f64,
        ra: Vec<f64>,
        dec: Vec<f64>,
        w: Vec<f64>,
    }

    /// Build a DAF/PCK file containing multiple Type 2 segments. Each
    /// segment owns one data record packed back-to-back after the name
    /// record; summaries appear in a single summary record in the
    /// order supplied.
    fn build_multi_segment_pck(segs: &[SegSpec]) -> NamedTempFile {
        let ss_doubles = 2 + 5_usize.div_ceil(2); // nd=2, ni=5 → 5 doubles
        let summary_bytes = ss_doubles * DOUBLE_BYTES; // 40

        let doubles_per_record = (RECORD_BYTES / DOUBLE_BYTES) as u32;
        let mut cursor = 1 + 3 * doubles_per_record;
        let mut segment_payloads: Vec<(u32, u32, Vec<f64>)> = Vec::new();
        for s in segs {
            assert_eq!(s.ra.len(), s.dec.len());
            assert_eq!(s.ra.len(), s.w.len());
            let n_coef = s.ra.len();
            let rsize = 2 + 3 * n_coef;
            let mut data = Vec::with_capacity(rsize + 4);
            data.push(s.mid);
            data.push(s.radius);
            data.extend_from_slice(&s.ra);
            data.extend_from_slice(&s.dec);
            data.extend_from_slice(&s.w);
            data.push(s.start_et);
            data.push(s.end_et - s.start_et);
            data.push(rsize as f64);
            data.push(1.0);
            let start_addr = cursor;
            let end_addr = start_addr + data.len() as u32 - 1;
            cursor = end_addr + 1;
            segment_payloads.push((start_addr, end_addr, data));
        }

        let mut record1 = vec![0u8; RECORD_BYTES];
        record1[0..8].copy_from_slice(b"DAF/PCK ");
        record1[8..12].copy_from_slice(&2u32.to_le_bytes());
        record1[12..16].copy_from_slice(&5u32.to_le_bytes());
        for b in &mut record1[16..76] {
            *b = b' ';
        }
        record1[76..80].copy_from_slice(&2u32.to_le_bytes());
        record1[80..84].copy_from_slice(&2u32.to_le_bytes());
        record1[88..96].copy_from_slice(b"LTL-IEEE");

        let mut record2 = vec![0u8; RECORD_BYTES];
        record2[16..24].copy_from_slice(&(segs.len() as f64).to_le_bytes());
        for (i, s) in segs.iter().enumerate() {
            let (start_addr, end_addr, _) = &segment_payloads[i];
            let soff = 24 + i * summary_bytes;
            record2[soff..soff + 8].copy_from_slice(&s.start_et.to_le_bytes());
            record2[soff + 8..soff + 16].copy_from_slice(&s.end_et.to_le_bytes());
            let int_start = soff + 2 * DOUBLE_BYTES;
            record2[int_start..int_start + 4].copy_from_slice(&s.body_frame.to_le_bytes());
            record2[int_start + 4..int_start + 8].copy_from_slice(&s.ref_frame.to_le_bytes());
            record2[int_start + 8..int_start + 12].copy_from_slice(&2i32.to_le_bytes());
            record2[int_start + 12..int_start + 16]
                .copy_from_slice(&(*start_addr as i32).to_le_bytes());
            record2[int_start + 16..int_start + 20]
                .copy_from_slice(&(*end_addr as i32).to_le_bytes());
        }

        let record3 = vec![b' '; RECORD_BYTES];

        let mut data_bytes = Vec::new();
        for (_, _, data) in &segment_payloads {
            for d in data {
                data_bytes.extend_from_slice(&d.to_le_bytes());
            }
        }
        while data_bytes.len() % RECORD_BYTES != 0 {
            data_bytes.push(0);
        }

        let mut all = Vec::new();
        all.extend_from_slice(&record1);
        all.extend_from_slice(&record2);
        all.extend_from_slice(&record3);
        all.extend_from_slice(&data_bytes);

        let mut tmp = NamedTempFile::new().unwrap();
        tmp.write_all(&all).unwrap();
        tmp.flush().unwrap();
        tmp
    }

    #[test]
    fn later_segment_wins_when_coverage_overlaps() {
        // Two segments for the same body frame with overlapping
        // coverage. The one added last (stored later in the file)
        // reports a different constant RA; the newest-first scan
        // means lookups in the overlap region return the later value.
        let earlier = SegSpec {
            body_frame: ITRF93_FRAME,
            ref_frame: ECLIPJ2000_FRAME,
            start_et: -100.0,
            end_et: 100.0,
            mid: 0.0,
            radius: 100.0,
            ra: vec![1.0, 0.0],
            dec: vec![0.0, 0.0],
            w: vec![0.0, 0.0],
        };
        let later = SegSpec {
            body_frame: ITRF93_FRAME,
            ref_frame: ECLIPJ2000_FRAME,
            start_et: -50.0,
            end_et: 50.0,
            mid: 0.0,
            radius: 50.0,
            ra: vec![2.0, 0.0],
            dec: vec![0.0, 0.0],
            w: vec![0.0, 0.0],
        };
        let tmp = build_multi_segment_pck(&[earlier, later]);
        let pck = PckFile::open(tmp.path()).expect("open");
        // In the overlap: later wins (RA=2.0).
        let overlap = pck.euler_state(ITRF93_FRAME, 0.0).expect("overlap");
        assert!(
            (overlap[0] - 2.0).abs() < 1e-14,
            "RA in overlap={}",
            overlap[0]
        );
        // Outside the later segment but inside the earlier: earlier wins (RA=1.0).
        let outside = pck.euler_state(ITRF93_FRAME, 75.0).expect("outside");
        assert!(
            (outside[0] - 1.0).abs() < 1e-14,
            "RA outside={}",
            outside[0]
        );
    }
}