astrodynamics-gnss 0.9.4

//! SP3 serialization — the inverse of the parser ([`super::Sp3::parse`]).
//!
//! Pure and deterministic: the same [`Sp3`] always produces byte-identical text.
//! No I/O. A read -> (merge) -> write pipeline round-trips: re-parsing the output
//! yields the same epochs, satellites, positions, and clocks to SP3 format
//! precision (mm / sub-ns). Header fields are derived from the product, never
//! hardcoded; the quality-metadata fields the reader does not retain (per-satellite
//! accuracy exponents, the `%f`/`%i` base descriptors) are emitted as standard
//! defaults.
//!
//! A satellite absent at an epoch is written as the SP3 missing-orbit sentinel
//! (`0.0 0.0 0.0`, bad clock), never a fabricated position — so a quarantined
//! `(sat, epoch)` cell from [`super::merge`] re-reads as missing, not zero.

use core::fmt::Write as _;

use astrodynamics::time::model::TimeScale;

use super::{
    Sp3, Sp3DataType, Sp3Version, BAD_CLOCK_US, CLOCK_RATE_TO_S_PER_S, DM_S_TO_M_S, KM_TO_M,
    US_TO_S,
};

/// Maximum SP3 satellite-id slots per `+` / `++` header line.
const SATS_PER_LINE: usize = 17;
/// SP3-c fixes five `+`/`++` lines (85 slots); SP3-d may use more.
const MIN_PLUS_LINES: usize = 5;

impl Sp3 {
    /// Serialize this product to standard SP3 text (the format named by its
    /// header version, `c` or `d`).
    ///
    /// Pure and deterministic. See the [module docs](self) for the round-trip
    /// and missing-satellite guarantees.
    pub fn to_sp3_string(&self) -> String {
        let mut out =
            String::with_capacity(self.epochs.len() * (self.header.satellites.len() + 4) * 61);
        self.write_header(&mut out);
        self.write_records(&mut out);
        out.push_str("EOF\n");
        out
    }

    fn write_header(&self, out: &mut String) {
        let h = &self.header;
        let version = match h.version {
            Sp3Version::A => 'a',
            Sp3Version::B => 'b',
            Sp3Version::C => 'c',
            Sp3Version::D => 'd',
        };
        let dtype = match h.data_type {
            Sp3DataType::Position => 'P',
            Sp3DataType::Velocity => 'V',
        };

        // Line 1: version/type, first-epoch calendar (cosmetic — the parser reads
        // epochs from the `*` lines), epoch count, data descriptor, coordinate
        // system, orbit type, agency. Columns match the parser's field offsets.
        let (y, mo, d, hh, mi, ss) = self
            .epochs
            .first()
            .map(julian_to_civil)
            .unwrap_or((2000, 1, 1, 0, 0, 0.0));
        let dt = format_calendar(y, mo, d, hh, mi, ss);
        let _ = writeln!(
            out,
            "#{version}{dtype}{dt} {n:>7} {data:<5}{coord:>6}{orbit:>4} {agency}",
            n = self.epochs.len(),
            data = "ORBIT",
            coord = h.coordinate_system,
            orbit = h.orbit_type,
            agency = h.agency,
        );

        // Line 2 (`##`): GPS week, seconds-of-week, epoch interval, MJD, MJD frac.
        let _ = writeln!(
            out,
            "## {wk:>4} {sow:15.8} {interval:14.8} {mjd:>5} {frac:.13}",
            wk = h.gnss_week,
            sow = h.seconds_of_week,
            interval = h.epoch_interval_s,
            mjd = h.mjd,
            frac = h.mjd_fraction,
        );

        // `+` satellite-id lines and `++` accuracy-exponent lines. The reader does
        // not retain per-satellite accuracy, so exponents are emitted as 0.
        let sats = &h.satellites;
        let n_lines = MIN_PLUS_LINES.max(sats.len().div_ceil(SATS_PER_LINE));
        for line in 0..n_lines {
            // `+` line: first carries the count in columns 3-5; all start ids at 9.
            if line == 0 {
                let _ = write!(out, "+  {:>3}   ", sats.len());
            } else {
                out.push_str("+        ");
            }
            for slot in 0..SATS_PER_LINE {
                match sats.get(line * SATS_PER_LINE + slot) {
                    Some(sat) => {
                        let _ = write!(out, "{sat}");
                    }
                    None => out.push_str("  0"),
                }
            }
            out.push('\n');
        }
        for _ in 0..n_lines {
            out.push_str("++       ");
            for _ in 0..SATS_PER_LINE {
                out.push_str("  0");
            }
            out.push('\n');
        }

        // `%c` descriptors — the first carries the time system at columns 9-11
        // (the only `%c` content the parser reads). `%f`/`%i` are standard
        // base/accuracy descriptors the parser skips.
        let tsys = sp3_time_label(h.time_scale);
        let _ = writeln!(
            out,
            "%c M  cc {tsys} ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc"
        );
        out.push_str("%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n");
        out.push_str("%f  1.2500000  1.025000000  0.00000000000  0.000000000000000\n");
        out.push_str("%f  0.0000000  0.000000000  0.00000000000  0.000000000000000\n");
        out.push_str("%i    0    0    0    0      0      0      0      0         0\n");
        out.push_str("%i    0    0    0    0      0      0      0      0         0\n");

        // Provenance comments (e.g. merge derivation). SP3 conventionally carries
        // at least one `/*` line.
        if self.comments.is_empty() {
            out.push_str("/* Generated by orbis\n");
        } else {
            for comment in &self.comments {
                let _ = writeln!(out, "/* {comment}");
            }
        }
    }

    fn write_records(&self, out: &mut String) {
        let with_velocity = matches!(self.header.data_type, Sp3DataType::Velocity);
        for (idx, epoch) in self.epochs.iter().enumerate() {
            let (y, mo, d, hh, mi, ss) = julian_to_civil(epoch);
            let _ = writeln!(out, "*  {}", format_calendar(y, mo, d, hh, mi, ss));

            let states = &self.states[idx];
            // Every header satellite gets a record at every epoch; an absent one
            // is the missing-orbit sentinel (so a quarantined cell is "missing",
            // never a fabricated zero position).
            for sat in &self.header.satellites {
                match states.get(sat) {
                    Some(state) => {
                        let p = state.position;
                        let clk = clock_field_us(state.clock_s);
                        let _ = writeln!(
                            out,
                            "P{sat}{:14.6}{:14.6}{:14.6}{clk:14.6}",
                            p.x_m / KM_TO_M,
                            p.y_m / KM_TO_M,
                            p.z_m / KM_TO_M,
                        );
                        if with_velocity {
                            if let Some(v) = state.velocity {
                                let rate = clock_rate_field(state.clock_rate_s_s);
                                let _ = writeln!(
                                    out,
                                    "V{sat}{:14.6}{:14.6}{:14.6}{rate:14.6}",
                                    v.vx_m_s / DM_S_TO_M_S,
                                    v.vy_m_s / DM_S_TO_M_S,
                                    v.vz_m_s / DM_S_TO_M_S,
                                );
                            }
                        }
                    }
                    None => {
                        let _ = writeln!(
                            out,
                            "P{sat}{:14.6}{:14.6}{:14.6}{:14.6}",
                            0.0, 0.0, 0.0, BAD_CLOCK_US
                        );
                    }
                }
            }
        }
    }
}

/// SP3 clock column (microseconds): the bad-clock sentinel when absent.
fn clock_field_us(clock_s: Option<f64>) -> f64 {
    match clock_s {
        Some(s) => s / US_TO_S,
        None => BAD_CLOCK_US,
    }
}

/// SP3 clock-rate column: the bad-clock sentinel when absent.
fn clock_rate_field(rate_s_s: Option<f64>) -> f64 {
    match rate_s_s {
        Some(r) => r / CLOCK_RATE_TO_S_PER_S,
        None => BAD_CLOCK_US,
    }
}

fn sp3_time_label(ts: TimeScale) -> &'static str {
    match ts {
        TimeScale::Gpst => "GPS",
        TimeScale::Gst => "GAL",
        TimeScale::Bdt => "BDT",
        TimeScale::Tai => "TAI",
        TimeScale::Utc => "UTC",
        // Other scales do not occur in SP3 products; default keeps output valid.
        _ => "GPS",
    }
}

/// `YYYY MM DD HH MM SS.SSSSSSSS` in the SP3 epoch-line / line-1 layout.
fn format_calendar(
    year: i64,
    month: i64,
    day: i64,
    hour: i64,
    minute: i64,
    seconds: f64,
) -> String {
    format!("{year:4} {month:>2} {day:>2} {hour:>2} {minute:>2} {seconds:11.8}")
}

/// Inverse of `super::civil_to_julian_split`: a Julian-date Instant back to the
/// civil `(year, month, day, hour, minute, seconds)` it was built from.
fn julian_to_civil(epoch: &astrodynamics::time::model::Instant) -> (i64, i64, i64, i64, i64, f64) {
    let split = match epoch.julian_date() {
        Some(s) => s,
        None => return (2000, 1, 1, 0, 0, 0.0),
    };
    // `jd_whole` is the `*.5` midnight boundary (jdn - 0.5); recover the integer
    // Julian Day Number, then the Fliegel-Van Flandern inverse for the calendar.
    let jdn = (split.jd_whole + 0.5).round() as i64;
    let l = jdn + 68_569;
    let n = 4 * l / 146_097;
    let l = l - (146_097 * n + 3) / 4;
    let i = 4_000 * (l + 1) / 1_461_001;
    let l = l - 1_461 * i / 4 + 31;
    let j = 80 * l / 2_447;
    let day = l - 2_447 * j / 80;
    let l = j / 11;
    let month = j + 2 - 12 * l;
    let year = 100 * (n - 49) + i + l;

    // Time of day from the fraction. Round to 10 ns (the SP3 seconds precision)
    // in integer arithmetic so float jitter never flips a minute boundary.
    let ticks = (split.fraction * 86_400.0 * 1.0e8).round() as i64;
    let hour = ticks / 360_000_000_000;
    let rem = ticks % 360_000_000_000;
    let minute = rem / 6_000_000_000;
    let seconds = (rem % 6_000_000_000) as f64 / 1.0e8;
    (year, month, day, hour, minute, seconds)
}