astrodynamics-gnss 0.8.0

//! RINEX 3.x and 4.xx navigation-message parsing (GPS LNAV, Galileo I/NAV and
//! F/NAV, BeiDou D1/D2).
//!
//! Version 4 wraps each record in a `> EPH|STO|EOP|ION SVNN MSG` frame marker but
//! keeps the same fixed-column broadcast-orbit layout, so the two versions share
//! the block parser; only the record grouping differs. CNAV-family messages
//! (CNAV/CNV1/CNV2/CNV3) reorder the orbit roster and are recognized but not
//! parsed.
//!
//! Reads broadcast ephemeris records out of a RINEX navigation file into the
//! typed [`BroadcastRecord`]s the [`crate::broadcast`] evaluator consumes. This
//! is deterministic byte-to-record parsing of a fixed-column text format, not a
//! float recipe: there is no 0-ULP claim here, and a small in-house parser is
//! used in preference to a heavyweight RINEX dependency (the published `rinex`
//! crate pulls ~90 transitive crates, including computational-geometry stacks,
//! for what is a fixed-width text read).
//!
//! Scope: the GPS, Galileo, and BeiDou Keplerian record layouts (eight lines:
//! the SV/epoch/clock line plus seven broadcast-orbit lines), plus the GLONASS
//! four-line state-vector layout (parsed by [`parse_glonass`] and evaluated by
//! the [`crate::glonass`] RK4 propagator, not the Keplerian path). Other
//! constellations' records (SBAS, QZSS) are recognized as record boundaries and
//! skipped, so a mixed file parses without error but yields only the supported
//! systems.

use crate::broadcast::{
    satellite_state, ClockPolynomial, ConstellationConstants, KeplerianElements, SECONDS_PER_WEEK,
};
use crate::glonass;
use crate::id::{GnssSatelliteId, GnssSystem};
use crate::spp::EphemerisSource;

/// Seconds from the GPS epoch (1980-01-06 00:00) to J2000 (2000-01-01 12:00).
/// Used to map an SPP J2000 second onto the GPS/Galileo continuous time the
/// broadcast elements are referenced to.
const GPS_EPOCH_TO_J2000_S: f64 = 630_763_200.0;

/// Fallback half-window (seconds, either side of `toe`) for a record that does
/// not broadcast a fit interval (Galileo, BeiDou). A coarse validity guard — a
/// stale or wrong-week product is off by at least a week, so this rejects it as
/// "no ephemeris" rather than silently extrapolating. GPS records carry an
/// explicit curve-fit interval (see [`BroadcastRecord::fit_interval_s`]) and use
/// half of that instead.
const MAX_EPHEMERIS_AGE_S: f64 = 4.0 * 3600.0;

/// GLONASS broadcast records are valid +/-15 minutes around their reference
/// epoch (the nominal half-hour upload cadence), so a query farther than this
/// reports no ephemeris rather than extrapolating the RK4 integration.
const GLONASS_MAX_AGE_S: f64 = 15.0 * 60.0;

/// GPS Time minus BeiDou Time (seconds): a constant +14 s (BDT = GPST - 14 s;
/// BDT was aligned to UTC at the BeiDou epoch, 2006-01-01, when GPS-UTC was 14 s).
/// So `bdt = gpst - GPST_MINUS_BDT_S`.
const GPST_MINUS_BDT_S: f64 = 14.0;

/// Seconds from the GPS epoch (1980-01-06) to the BeiDou epoch (2006-01-01),
/// which are exactly 1356 weeks apart, so the week boundaries align.
const BDS_EPOCH_MINUS_GPS_EPOCH_S: f64 = 1356.0 * SECONDS_PER_WEEK;

/// Which broadcast navigation message a record carries.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum NavMessage {
    /// GPS legacy navigation message.
    GpsLnav,
    /// Galileo integrity navigation message (E1/E5b dual, E1 single-frequency).
    GalileoInav,
    /// Galileo F/NAV message (E5a).
    GalileoFnav,
    /// BeiDou D1 message (MEO/IGSO satellites).
    BeidouD1,
    /// BeiDou D2 message (geostationary satellites).
    BeidouD2,
}

/// Whether a BeiDou PRN is a geostationary satellite (BDS-2 C01-C05, BDS-3
/// C59-C61), which take the geostationary orbit-evaluation branch.
pub fn is_beidou_geo(sat: GnssSatelliteId) -> bool {
    sat.system == GnssSystem::BeiDou && (sat.prn <= 5 || (59..=61).contains(&sat.prn))
}

/// A Klobuchar-8 broadcast ionosphere coefficient set (the eight alpha/beta
/// values transmitted by GPS and BeiDou; the same model serves both, evaluated
/// per carrier — see [`crate::ionex::klobuchar_native`]).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct KlobucharAlphaBeta {
    /// Cosine-amplitude polynomial coefficients (a0..a3).
    pub alpha: [f64; 4],
    /// Period polynomial coefficients (b0..b3).
    pub beta: [f64; 4],
}

/// Broadcast ionosphere-correction coefficients from a RINEX header's
/// `IONOSPHERIC CORR` lines.
///
/// Only the Klobuchar-8 sets the toolkit's ionosphere model uses are captured:
/// GPS (`GPSA`/`GPSB`) and BeiDou (`BDSA`/`BDSB`). Both feed the same Klobuchar
/// algorithm (frequency-scaled per carrier), so a BeiDou-only solve can use
/// `beidou`'s native coefficients rather than the GPS set. The Galileo NeQuick
/// (`GAL`) coefficients are a different model and are not captured here; QZSS and
/// NavIC Klobuchar sets are likewise not retained.
#[derive(Debug, Clone, Copy, PartialEq, Default)]
pub struct IonoCorrections {
    /// GPS broadcast Klobuchar coefficients (`GPSA`/`GPSB`), if present.
    pub gps: Option<KlobucharAlphaBeta>,
    /// BeiDou broadcast Klobuchar coefficients (`BDSA`/`BDSB`), if present.
    pub beidou: Option<KlobucharAlphaBeta>,
}

/// One parsed GLONASS broadcast record: a PZ-90.11 ECEF state vector and the
/// clock terms, evaluated by the [`crate::glonass`] RK4 propagator (GLONASS is
/// not Keplerian, so it does not use [`BroadcastRecord`]).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GlonassRecord {
    /// The transmitting satellite.
    pub satellite_id: GnssSatelliteId,
    /// Reference epoch as seconds past J2000 in **UTC** (leap-second-independent;
    /// the store adds the GPS−UTC offset to compare with the GPST-aligned query).
    pub toe_utc_j2000_s: f64,
    /// PZ-90.11 ECEF position at the reference epoch (meters).
    pub pos_m: [f64; 3],
    /// PZ-90.11 ECEF velocity at the reference epoch (meters/second).
    pub vel_m_s: [f64; 3],
    /// Lunisolar acceleration at the reference epoch (meters/second^2).
    pub acc_m_s2: [f64; 3],
    /// Clock bias broadcast field (−TauN, seconds).
    pub clk_bias: f64,
    /// Relative frequency offset (+GammaN, dimensionless).
    pub gamma_n: f64,
    /// Satellite health (0 is healthy).
    pub sv_health: f64,
    /// FDMA frequency-channel number.
    pub freq_channel: i32,
}

/// One parsed broadcast navigation record.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct BroadcastRecord {
    /// The transmitting satellite.
    pub satellite_id: GnssSatelliteId,
    /// The navigation message the record carries.
    pub message: NavMessage,
    /// GPS/Galileo continuous week number (from the broadcast record).
    pub week: u32,
    /// Keplerian orbital elements (`toe_sow` is seconds of week).
    pub elements: KeplerianElements,
    /// Clock polynomial (`toc_sow` is the record's own epoch, seconds of week).
    pub clock: ClockPolynomial,
    /// Single-frequency broadcast group delay (GPS L1 TGD / Galileo E1-E5a BGD), seconds.
    pub group_delay_s: f64,
    /// Satellite health word (0 is healthy for the GPS/Galileo nominal case).
    pub sv_health: f64,
    /// Signal-in-space accuracy: GPS URA (m) / Galileo SISA (m).
    pub sv_accuracy_m: f64,
    /// GPS curve-fit interval in seconds, centered on `toe` (IS-GPS-200): the
    /// record is valid for `toe ± fit_interval_s / 2`. `None` for Galileo and
    /// BeiDou, which do not broadcast a fit interval in the RINEX record; those
    /// fall back to [`MAX_EPHEMERIS_AGE_S`].
    pub fit_interval_s: Option<f64>,
}

impl BroadcastRecord {
    /// The per-constellation constants this record evaluates with.
    pub const fn constants(&self) -> ConstellationConstants {
        match self.satellite_id.system {
            GnssSystem::Galileo => ConstellationConstants::GALILEO,
            GnssSystem::BeiDou => ConstellationConstants::BEIDOU,
            // GPS (and any other Keplerian system) use the GPS constants.
            _ => ConstellationConstants::GPS,
        }
    }
}

/// Why a RINEX NAV file could not be parsed.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum NavParseError {
    /// The header did not declare a supported RINEX 3 navigation file.
    UnsupportedHeader(String),
    /// No `END OF HEADER` line was found.
    MissingHeaderEnd,
    /// A record was shorter than its message layout requires.
    TruncatedRecord(String),
    /// A required numeric field was missing or unparseable.
    BadField {
        /// The satellite whose record holds the bad field.
        satellite: String,
        /// Which field failed.
        field: &'static str,
    },
}

impl core::fmt::Display for NavParseError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            NavParseError::UnsupportedHeader(s) => write!(f, "unsupported RINEX NAV header: {s}"),
            NavParseError::MissingHeaderEnd => write!(f, "no END OF HEADER line"),
            NavParseError::TruncatedRecord(s) => write!(f, "truncated navigation record for {s}"),
            NavParseError::BadField { satellite, field } => {
                write!(f, "bad/missing {field} field in record for {satellite}")
            }
        }
    }
}

impl std::error::Error for NavParseError {}

/// Parse a RINEX 3.x or 4.xx navigation file into the supported (GPS, Galileo,
/// BeiDou) Keplerian records.
///
/// Records of other constellations (GLONASS state-vector, SBAS) are skipped, as
/// are version-4 CNAV-family messages (CNAV/CNV1/CNV2/CNV3): their broadcast-orbit
/// roster reorders the fixed columns (`t_op` for `toe`, `wn_op` for `week`, extra
/// `adot`/`deltaN0Dot` terms), so they are recognized but not parsed rather than
/// fed wrong values. The records are returned in file order; selection (by epoch,
/// health, message type) is the caller's job.
pub fn parse_nav(text: &str) -> Result<Vec<BroadcastRecord>, NavParseError> {
    let mut lines = text.lines();
    let major = verify_and_skip_header(&mut lines)?;
    if major >= 4 {
        parse_nav_v4(lines)
    } else {
        parse_nav_v3(lines)
    }
}

/// Version-3 body: a record starts at a line whose first three columns are a
/// system letter followed by two digits; continuation lines are column-indented.
fn parse_nav_v3<'a, I>(lines: I) -> Result<Vec<BroadcastRecord>, NavParseError>
where
    I: Iterator<Item = &'a str>,
{
    let mut blocks: Vec<Vec<&str>> = Vec::new();
    for line in lines {
        if is_record_start(line) {
            blocks.push(vec![line]);
        } else if let Some(last) = blocks.last_mut() {
            last.push(line);
        }
    }

    let mut records = Vec::new();
    for block in &blocks {
        let letter = block[0].as_bytes()[0] as char;
        match GnssSystem::from_letter(letter) {
            Some(GnssSystem::Gps) | Some(GnssSystem::Galileo) | Some(GnssSystem::BeiDou) => {
                records.push(parse_keplerian_block(block)?);
            }
            // Recognized boundary, unsupported model (GLONASS state-vector, SBAS): skip.
            _ => {}
        }
    }
    Ok(records)
}

/// Version-4 body: each record is introduced by a `> EPH|STO|EOP|ION SVNN MSG`
/// frame marker. Only `EPH` frames carrying a supported Keplerian message are
/// parsed; the broadcast-orbit lines that follow the marker have the same
/// fixed-column layout as version 3, so they reuse [`parse_keplerian_block`].
/// The message type is taken from the marker token (so I/NAV vs F/NAV and D1 vs
/// D2 are explicit, not inferred). STO/EOP/ION frames, other constellations, and
/// CNAV-family messages are skipped.
fn parse_nav_v4<'a, I>(lines: I) -> Result<Vec<BroadcastRecord>, NavParseError>
where
    I: Iterator<Item = &'a str>,
{
    // Group by marker line: each frame is its marker plus the body lines up to
    // the next marker.
    let mut frames: Vec<(&str, Vec<&str>)> = Vec::new();
    for line in lines {
        if is_v4_frame_marker(line) {
            frames.push((line, Vec::new()));
        } else if let Some((_, body)) = frames.last_mut() {
            body.push(line);
        }
    }

    let mut records = Vec::new();
    for (marker, body) in &frames {
        let Some((frame_type, sv, msg_token)) = parse_v4_marker(marker) else {
            continue;
        };
        if frame_type != "EPH" {
            continue; // STO/EOP/ION carry no ephemeris.
        }
        let letter = sv.as_bytes().first().map(|b| *b as char).unwrap_or(' ');
        let supported = matches!(
            GnssSystem::from_letter(letter),
            Some(GnssSystem::Gps) | Some(GnssSystem::Galileo) | Some(GnssSystem::BeiDou)
        );
        if !supported {
            continue; // GLONASS/SBAS/QZSS/NavIC: not a supported Keplerian system.
        }
        // Only messages whose orbit roster matches the version-3 column layout
        // are parsed; CNAV-family (and anything unrecognized) is skipped.
        if let Some(message) = nav_message_from_v4_token(msg_token) {
            let mut record = parse_keplerian_block(body)?;
            record.message = message;
            records.push(record);
        }
    }
    Ok(records)
}

/// Whether a version-4 line is a frame marker (`> ...`).
fn is_v4_frame_marker(line: &str) -> bool {
    line.starts_with("> ")
}

/// Split a version-4 frame marker `> EPH G01 LNAV` into (frame type, SV, message
/// token), or `None` if it is malformed. Mirrors the RINEX-4 marker layout:
/// `>` then the 4-column frame class, the SV, and the message-type token.
fn parse_v4_marker(line: &str) -> Option<(&str, &str, &str)> {
    let rest = line.strip_prefix('>')?;
    let mut fields = rest.split_whitespace();
    let frame_type = fields.next()?;
    let sv = fields.next()?;
    let msg_token = fields.next()?;
    Some((frame_type, sv, msg_token))
}

/// Map a version-4 EPH message token to the [`NavMessage`] for the supported
/// Keplerian messages, or `None` for a message whose orbit layout does not match
/// the version-3 columns (CNAV-family) or is otherwise unsupported here.
fn nav_message_from_v4_token(token: &str) -> Option<NavMessage> {
    match token {
        "LNAV" => Some(NavMessage::GpsLnav),
        "INAV" => Some(NavMessage::GalileoInav),
        "FNAV" => Some(NavMessage::GalileoFnav),
        "D1" => Some(NavMessage::BeidouD1),
        "D2" => Some(NavMessage::BeidouD2),
        _ => None,
    }
}

/// Parse the broadcast ionosphere coefficients from a RINEX header's
/// `IONOSPHERIC CORR` lines (GPS `GPSA`/`GPSB` and BeiDou `BDSA`/`BDSB`).
///
/// Header-only and version-agnostic: it scans up to `END OF HEADER` and ignores
/// the body. A label with all four values present yields the coefficient set; a
/// missing label yields `None` for that system. Parsing is deterministic text,
/// not a 0-ULP target.
pub fn parse_iono_corrections(text: &str) -> IonoCorrections {
    // The IONOSPHERIC CORR line is `A4,1X,4(D12.4)`: a 4-char label, a space,
    // then four coefficients in 12-wide columns.
    let row = |line: &str| -> Option<[f64; 4]> {
        Some([
            parse_f64(line, 5, 17)?,
            parse_f64(line, 17, 29)?,
            parse_f64(line, 29, 41)?,
            parse_f64(line, 41, 53)?,
        ])
    };
    let (mut gpsa, mut gpsb, mut bdsa, mut bdsb) = (None, None, None, None);
    for line in text.lines() {
        if line.contains("END OF HEADER") {
            break;
        }
        if !line.contains("IONOSPHERIC CORR") {
            continue;
        }
        match line.get(0..4).map(str::trim) {
            Some("GPSA") => gpsa = row(line),
            Some("GPSB") => gpsb = row(line),
            Some("BDSA") => bdsa = row(line),
            Some("BDSB") => bdsb = row(line),
            _ => {}
        }
    }
    let pair = |a: Option<[f64; 4]>, b: Option<[f64; 4]>| match (a, b) {
        (Some(alpha), Some(beta)) => Some(KlobucharAlphaBeta { alpha, beta }),
        _ => None,
    };
    IonoCorrections {
        gps: pair(gpsa, gpsb),
        beidou: pair(bdsa, bdsb),
    }
}

/// The leap-second count (GPS − UTC) from the header's `LEAP SECONDS` line, used
/// to map a GLONASS (UTC) reference epoch onto the GPST-aligned query time. The
/// value is the first field; `None` if the line is absent.
pub fn parse_leap_seconds(text: &str) -> Option<f64> {
    for line in text.lines() {
        if line.contains("END OF HEADER") {
            break;
        }
        if line.contains("LEAP SECONDS") {
            return line.get(0..6).and_then(|s| s.trim().parse::<f64>().ok());
        }
    }
    None
}

/// Days from the civil epoch 1970-01-01 to a proleptic-Gregorian date
/// (Howard Hinnant's algorithm), used to place a GLONASS UTC epoch on an
/// absolute timeline.
fn days_from_civil(y: i64, m: i64, d: i64) -> i64 {
    let y = if m <= 2 { y - 1 } else { y };
    let era = if y >= 0 { y } else { y - 399 } / 400;
    let yoe = y - era * 400;
    let doy = (153 * (if m > 2 { m - 3 } else { m + 9 }) + 2) / 5 + d - 1;
    let doe = yoe * 365 + yoe / 4 - yoe / 100 + doy;
    era * 146097 + doe - 719468
}

/// Seconds from the J2000 epoch (2000-01-01 12:00) to a UTC calendar instant.
/// (2000-01-01 is 10957 days after 1970-01-01; J2000 is noon that day.)
fn j2000_seconds_utc(y: i64, mo: i64, d: i64, h: i64, mi: i64, s: i64) -> f64 {
    let day = days_from_civil(y, mo, d) - 10957;
    (day * 86_400 + h * 3600 + mi * 60 + s - 43_200) as f64
}

/// Parse the GLONASS epoch line (`Rnn YYYY MM DD HH MM SS`) to a UTC second past
/// J2000.
fn parse_glonass_epoch(l0: &str) -> Option<f64> {
    let y: i64 = field(l0, 4, 8)?.parse().ok()?;
    let mo: i64 = field(l0, 9, 11)?.parse().ok()?;
    let d: i64 = field(l0, 12, 14)?.parse().ok()?;
    let h: i64 = field(l0, 15, 17)?.parse().ok()?;
    let mi: i64 = field(l0, 18, 20)?.parse().ok()?;
    let s: i64 = field(l0, 21, 23)?.trim().parse().ok()?;
    Some(j2000_seconds_utc(y, mo, d, h, mi, s))
}

/// Parse a 4-line RINEX 3 GLONASS record block into a [`GlonassRecord`]
/// (km/(km/s)/(km/s^2) state converted to SI). A missing or unparseable field is
/// a [`NavParseError`], not a silently dropped record.
fn parse_glonass_block(block: &[&str]) -> Result<GlonassRecord, NavParseError> {
    let l0 = block[0];
    let sat = l0.get(0..3).unwrap_or("").trim().to_string();
    if block.len() < 4 {
        return Err(NavParseError::TruncatedRecord(sat));
    }
    let bad = |what: &'static str| NavParseError::BadField {
        satellite: sat.clone(),
        field: what,
    };
    let prn: u8 = l0
        .get(1..3)
        .and_then(|s| s.trim().parse().ok())
        .ok_or_else(|| bad("prn"))?;
    let satellite_id = GnssSatelliteId::new(GnssSystem::Glonass, prn);
    let toe_utc_j2000_s = parse_glonass_epoch(l0).ok_or_else(|| bad("epoch"))?;
    let clk_bias = parse_f64(l0, 23, 42).ok_or_else(|| bad("clock bias"))?;
    let gamma_n = parse_f64(l0, 42, 61).ok_or_else(|| bad("gamma_n"))?;
    let o1 = orbit_row(block[1]);
    let o2 = orbit_row(block[2]);
    let o3 = orbit_row(block[3]);
    let km = |v: Option<f64>, what: &'static str| v.map(|x| x * 1000.0).ok_or_else(|| bad(what));
    Ok(GlonassRecord {
        satellite_id,
        toe_utc_j2000_s,
        pos_m: [km(o1[0], "x")?, km(o2[0], "y")?, km(o3[0], "z")?],
        vel_m_s: [km(o1[1], "vx")?, km(o2[1], "vy")?, km(o3[1], "vz")?],
        acc_m_s2: [km(o1[2], "ax")?, km(o2[2], "ay")?, km(o3[2], "az")?],
        clk_bias,
        gamma_n,
        sv_health: o1[3].unwrap_or(0.0),
        freq_channel: o2[3].unwrap_or(0.0) as i32,
    })
}

/// Parse all GLONASS (`R`) records from a RINEX 3.x navigation file, in file
/// order; selection is the caller's job. A malformed supported record is a
/// [`NavParseError`] rather than a silently dropped one. (Version-4 GLONASS
/// frames are not yet parsed.)
pub fn parse_glonass(text: &str) -> Result<Vec<GlonassRecord>, NavParseError> {
    let mut lines = text.lines();
    if verify_and_skip_header(&mut lines).is_err() {
        return Ok(Vec::new());
    }
    let mut blocks: Vec<Vec<&str>> = Vec::new();
    for line in lines {
        if is_record_start(line) {
            blocks.push(vec![line]);
        } else if let Some(last) = blocks.last_mut() {
            last.push(line);
        }
    }
    blocks
        .iter()
        .filter(|b| b[0].starts_with('R'))
        .map(|b| parse_glonass_block(b))
        .collect()
}

/// Skip the header, returning the RINEX major version (3 or 4). Both share the
/// fixed-column orbit layout; version 4 wraps each record in a frame marker line
/// (see [`parse_v4_marker`]), which is why `parse_nav` dispatches on it.
fn verify_and_skip_header<'a, I>(lines: &mut I) -> Result<u8, NavParseError>
where
    I: Iterator<Item = &'a str>,
{
    let mut major: Option<u8> = None;
    for line in lines.by_ref() {
        if line.contains("RINEX VERSION / TYPE") {
            // Column 0-8 holds the version; column 20 the file type ('N' = NAV).
            let version = line.get(0..9).unwrap_or("").trim();
            let detected = match version.chars().next() {
                Some('3') => Some(3u8),
                Some('4') => Some(4u8),
                _ => None,
            };
            let is_nav = line.get(20..21) == Some("N");
            match (detected, is_nav) {
                (Some(v), true) => major = Some(v),
                _ => {
                    return Err(NavParseError::UnsupportedHeader(
                        line.trim_end().to_string(),
                    ))
                }
            }
        }
        if line.contains("END OF HEADER") {
            return major.ok_or_else(|| {
                NavParseError::UnsupportedHeader("no RINEX VERSION / TYPE".to_string())
            });
        }
    }
    Err(NavParseError::MissingHeaderEnd)
}

fn is_record_start(line: &str) -> bool {
    let b = line.as_bytes();
    b.len() >= 3 && b[0].is_ascii_alphabetic() && b[1].is_ascii_digit() && b[2].is_ascii_digit()
}

/// A fixed-column numeric field, trimmed; `None` if out of range or blank.
fn field(line: &str, start: usize, end: usize) -> Option<&str> {
    let s = line.get(start..end.min(line.len()))?.trim();
    if s.is_empty() {
        None
    } else {
        Some(s)
    }
}

/// Parse a RINEX numeric field, tolerating the Fortran `D`/`d` exponent marker.
fn parse_f64(line: &str, start: usize, end: usize) -> Option<f64> {
    let s = field(line, start, end)?;
    s.replace(['D', 'd'], "e").parse::<f64>().ok()
}

/// The four broadcast-orbit values of a continuation line (columns 4/23/42/61).
fn orbit_row(line: &str) -> [Option<f64>; 4] {
    [
        parse_f64(line, 4, 23),
        parse_f64(line, 23, 42),
        parse_f64(line, 42, 61),
        parse_f64(line, 61, 80),
    ]
}

/// Seconds-of-week of a calendar epoch in its own system time (Sunday 00:00
/// origin). Sakamoto's day-of-week gives 0 = Sunday = the GPS/GST day-of-week.
fn epoch_seconds_of_week(
    year: i64,
    month: i64,
    day: i64,
    hour: i64,
    minute: i64,
    second: i64,
) -> f64 {
    const T: [i64; 12] = [0, 3, 2, 5, 0, 3, 5, 1, 4, 6, 2, 4];
    let y = if month < 3 { year - 1 } else { year };
    let dow = (y + y / 4 - y / 100 + y / 400 + T[(month - 1) as usize] + day).rem_euclid(7);
    (dow * 86_400 + hour * 3600 + minute * 60 + second) as f64
}

fn parse_keplerian_block(block: &[&str]) -> Result<BroadcastRecord, NavParseError> {
    let l0 = block[0];
    let sat = l0.get(0..3).unwrap_or("").trim().to_string();
    if block.len() < 8 {
        return Err(NavParseError::TruncatedRecord(sat));
    }
    let bad = |what: &'static str| NavParseError::BadField {
        satellite: sat.clone(),
        field: what,
    };

    let letter = l0.as_bytes()[0] as char;
    let system = GnssSystem::from_letter(letter).ok_or_else(|| bad("system"))?;
    let prn: u8 = l0
        .get(1..3)
        .and_then(|s| s.trim().parse().ok())
        .ok_or_else(|| bad("prn"))?;
    let satellite_id = GnssSatelliteId::new(system, prn);

    // Clock line: epoch (-> toc) and the af0/af1/af2 polynomial.
    let toc_sow = parse_toc(l0).ok_or_else(|| bad("toc epoch"))?;
    let af0 = parse_f64(l0, 23, 42).ok_or_else(|| bad("af0"))?;
    let af1 = parse_f64(l0, 42, 61).ok_or_else(|| bad("af1"))?;
    let af2 = parse_f64(l0, 61, 80).ok_or_else(|| bad("af2"))?;

    let o1 = orbit_row(block[1]);
    let o2 = orbit_row(block[2]);
    let o3 = orbit_row(block[3]);
    let o4 = orbit_row(block[4]);
    let o5 = orbit_row(block[5]);
    let o6 = orbit_row(block[6]);

    let g = |v: Option<f64>, what: &'static str| v.ok_or_else(|| bad(what));

    let elements = KeplerianElements {
        crs: g(o1[1], "crs")?,
        delta_n: g(o1[2], "deltaN")?,
        m0: g(o1[3], "m0")?,
        cuc: g(o2[0], "cuc")?,
        e: g(o2[1], "e")?,
        cus: g(o2[2], "cus")?,
        sqrt_a: g(o2[3], "sqrtA")?,
        toe_sow: g(o3[0], "toe")?,
        cic: g(o3[1], "cic")?,
        omega0: g(o3[2], "omega0")?,
        cis: g(o3[3], "cis")?,
        i0: g(o4[0], "i0")?,
        crc: g(o4[1], "crc")?,
        omega: g(o4[2], "omega")?,
        omega_dot: g(o4[3], "omegaDot")?,
        idot: g(o5[0], "idot")?,
    };
    let clock = ClockPolynomial {
        af0,
        af1,
        af2,
        toc_sow,
    };

    let week = g(o5[2], "week")? as u32;
    let sv_accuracy_m = g(o6[0], "accuracy")?;
    let sv_health = g(o6[1], "health")?;
    // ORBIT-6 field 3 (index 2): GPS L1 C/A TGD, or Galileo E1-E5a BGD.
    let group_delay_s = g(o6[2], "group delay")?;

    let message = match system {
        GnssSystem::Galileo => galileo_message(o5[1]),
        GnssSystem::BeiDou => {
            if is_beidou_geo(satellite_id) {
                NavMessage::BeidouD2
            } else {
                NavMessage::BeidouD1
            }
        }
        _ => NavMessage::GpsLnav,
    };

    // Only GPS LNAV broadcasts a curve-fit interval (ORBIT-7 field 2); Galileo
    // and BeiDou leave that column blank or spare, so they carry no fit interval.
    let fit_interval_s = match system {
        GnssSystem::Gps => Some(gps_fit_interval_s(block[7]).map_err(|()| bad("fit interval"))?),
        _ => None,
    };

    Ok(BroadcastRecord {
        satellite_id,
        message,
        week,
        elements,
        clock,
        group_delay_s,
        sv_health,
        sv_accuracy_m,
        fit_interval_s,
    })
}

/// The GPS curve-fit interval in seconds from the ORBIT-7 fit-interval field,
/// which RINEX records in hours. Per IS-GPS-200 the value is the total interval
/// centered on `toe`; a zero or absent field denotes the nominal four hours.
///
/// A blank/absent field is the legitimate nominal case (some products omit it);
/// a present but non-numeric field is a malformed record, reported as `Err` so
/// the caller can raise the same `BadField` error as for other numeric fields
/// rather than silently substituting four hours.
fn gps_fit_interval_s(orbit7: &str) -> Result<f64, ()> {
    let hours = match field(orbit7, 23, 42) {
        None => 0.0,
        Some(_) => parse_f64(orbit7, 23, 42).ok_or(())?,
    };
    let hours = if hours == 0.0 { 4.0 } else { hours };
    Ok(hours * 3600.0)
}

/// Classify a Galileo record from its data-source word (orbit-5 field 1): bit 9
/// (0x200) is the I/NAV message, bit 8 (0x100) is F/NAV.
fn galileo_message(data_sources: Option<f64>) -> NavMessage {
    let word = data_sources.map(|v| v as u32).unwrap_or(0);
    if word & 0x200 != 0 {
        NavMessage::GalileoInav
    } else if word & 0x100 != 0 {
        NavMessage::GalileoFnav
    } else {
        // No source bit set: default to I/NAV (the operational E1 message).
        NavMessage::GalileoInav
    }
}

/// Parse the clock reference epoch from the SV/epoch line into seconds of week.
fn parse_toc(l0: &str) -> Option<f64> {
    let year: i64 = field(l0, 4, 8)?.parse().ok()?;
    let month: i64 = field(l0, 9, 11)?.parse().ok()?;
    let day: i64 = field(l0, 12, 14)?.parse().ok()?;
    let hour: i64 = field(l0, 15, 17)?.parse().ok()?;
    let minute: i64 = field(l0, 18, 20)?.parse().ok()?;
    let second: i64 = field(l0, 21, 23)?.trim().parse().ok()?;
    Some(epoch_seconds_of_week(
        year, month, day, hour, minute, second,
    ))
}

/// A queryable set of parsed broadcast records, usable as an SPP
/// [`EphemerisSource`].
///
/// For a satellite and epoch it selects the record whose reference time `toe` is
/// nearest in **continuous** GPS time (week number times the week length plus
/// the seconds of week, not the seconds of week alone), and rejects the query as
/// having no ephemeris if it falls outside that record's validity window — half
/// the broadcast GPS curve-fit interval, or the coarse [`MAX_EPHEMERIS_AGE_S`]
/// fallback for systems that do not broadcast one — so a stale or wrong-week
/// product cannot silently produce a position.
///
/// [`from_nav`](BroadcastStore::from_nav) applies a default usability policy:
/// only healthy GPS LNAV and healthy Galileo I/NAV records are kept (F/NAV and
/// unhealthy satellites are dropped). [`new`](BroadcastStore::new) keeps records
/// verbatim for callers that want their own policy.
pub struct BroadcastStore {
    records: Vec<BroadcastRecord>,
    glonass: Vec<GlonassRecord>,
    leap_seconds: Option<f64>,
    iono: IonoCorrections,
}

impl BroadcastStore {
    /// Build a store from already-parsed Keplerian records, verbatim (no policy
    /// filter, no GLONASS records, no leap-second offset, and no ionosphere
    /// coefficients; use [`from_nav`](Self::from_nav) to capture those).
    pub fn new(records: Vec<BroadcastRecord>) -> Self {
        Self {
            records,
            glonass: Vec::new(),
            leap_seconds: None,
            iono: IonoCorrections::default(),
        }
    }

    /// Parse a RINEX 3.x/4.xx navigation file and keep the records usable for
    /// single-frequency positioning: healthy GPS LNAV, Galileo I/NAV, BeiDou
    /// D1/D2, and healthy GLONASS. The header's broadcast ionosphere coefficients
    /// (see [`iono_corrections`](Self::iono_corrections)) and leap-second offset
    /// are captured.
    pub fn from_nav(text: &str) -> Result<Self, NavParseError> {
        let records = parse_nav(text)?
            .into_iter()
            .filter(Self::is_default_usable)
            .collect();
        let glonass = parse_glonass(text)?
            .into_iter()
            .filter(|r| r.sv_health == 0.0)
            .collect();
        Ok(Self {
            records,
            glonass,
            leap_seconds: parse_leap_seconds(text),
            iono: parse_iono_corrections(text),
        })
    }

    /// The broadcast ionosphere coefficients parsed from the navigation header
    /// (GPS `GPSA`/`GPSB` and BeiDou `BDSA`/`BDSB`). Empty for a store built with
    /// [`new`](Self::new).
    pub fn iono_corrections(&self) -> IonoCorrections {
        self.iono
    }

    /// The held GLONASS records.
    pub fn glonass_records(&self) -> &[GlonassRecord] {
        &self.glonass
    }

    /// The default usability policy: healthy and a single-frequency-appropriate
    /// message — GPS LNAV, Galileo I/NAV (the E1 message), and BeiDou D1/D2;
    /// Galileo F/NAV and unhealthy satellites are excluded.
    fn is_default_usable(r: &BroadcastRecord) -> bool {
        r.sv_health == 0.0
            && matches!(
                r.message,
                NavMessage::GpsLnav
                    | NavMessage::GalileoInav
                    | NavMessage::BeidouD1
                    | NavMessage::BeidouD2
            )
    }

    /// The held records.
    pub fn records(&self) -> &[BroadcastRecord] {
        &self.records
    }

    /// Keep only the records matching a predicate (e.g. a custom message/health
    /// policy on a store built with [`new`](BroadcastStore::new)).
    pub fn retain(&mut self, keep: impl FnMut(&BroadcastRecord) -> bool) {
        self.records.retain(keep);
    }

    /// Continuous GPS time (seconds since the GPS epoch) of a record's `toe`.
    fn toe_continuous_s(rec: &BroadcastRecord) -> f64 {
        f64::from(rec.week) * SECONDS_PER_WEEK + rec.elements.toe_sow
    }

    /// The half-validity window (seconds either side of `toe`) for a record: half
    /// the broadcast GPS fit interval, or [`MAX_EPHEMERIS_AGE_S`] when no fit
    /// interval is broadcast (Galileo, BeiDou).
    fn half_window_s(rec: &BroadcastRecord) -> f64 {
        match rec.fit_interval_s {
            Some(fit) => fit / 2.0,
            None => MAX_EPHEMERIS_AGE_S,
        }
    }

    /// The record for `sat` whose `toe` is nearest `t_continuous_s` (seconds
    /// since the GPS epoch) **among those whose validity window covers the
    /// query** (see [`half_window_s`](Self::half_window_s)). Filtering by
    /// validity before choosing the nearest means a query just past one record's
    /// fit interval can still be served by a farther record whose own window is
    /// wide enough, rather than being rejected outright.
    fn select(&self, sat: GnssSatelliteId, t_continuous_s: f64) -> Option<&BroadcastRecord> {
        self.records
            .iter()
            .filter(|r| r.satellite_id == sat)
            .filter(|r| {
                (t_continuous_s - Self::toe_continuous_s(r)).abs() <= Self::half_window_s(r)
            })
            .min_by(|a, b| {
                let da = (t_continuous_s - Self::toe_continuous_s(a)).abs();
                let db = (t_continuous_s - Self::toe_continuous_s(b)).abs();
                da.partial_cmp(&db).unwrap_or(core::cmp::Ordering::Equal)
            })
    }

    /// The GLONASS record for `sat` nearest the GPST-aligned query `t_j2000_s`
    /// (within [`GLONASS_MAX_AGE_S`]), with `tk` = query − the record's reference
    /// epoch in GPS time. Returns `None` if no leap-second offset was parsed (the
    /// GLONASS UTC epoch then cannot be placed on the GPST timeline).
    fn select_glonass(
        &self,
        sat: GnssSatelliteId,
        t_j2000_s: f64,
    ) -> Option<(&GlonassRecord, f64)> {
        let leap = self.leap_seconds?;
        let toe_gpst = |r: &GlonassRecord| r.toe_utc_j2000_s + leap;
        let rec = self
            .glonass
            .iter()
            .filter(|r| r.satellite_id == sat)
            .min_by(|a, b| {
                let da = (t_j2000_s - toe_gpst(a)).abs();
                let db = (t_j2000_s - toe_gpst(b)).abs();
                da.partial_cmp(&db).unwrap_or(core::cmp::Ordering::Equal)
            })?;
        let tk = t_j2000_s - toe_gpst(rec);
        if tk.abs() <= GLONASS_MAX_AGE_S {
            Some((rec, tk))
        } else {
            None
        }
    }
}

impl EphemerisSource for BroadcastStore {
    fn position_clock_at_j2000_s(
        &self,
        sat: GnssSatelliteId,
        t_j2000_s: f64,
    ) -> Option<([f64; 3], f64)> {
        // GLONASS is not Keplerian: integrate its broadcast state vector with the
        // RK4 propagator. Its reference epoch is UTC, mapped onto the GPST-aligned
        // query via the parsed leap-second offset.
        if sat.system == GnssSystem::Glonass {
            let (rec, tk) = self.select_glonass(sat, t_j2000_s)?;
            let state0 = [
                rec.pos_m[0],
                rec.pos_m[1],
                rec.pos_m[2],
                rec.vel_m_s[0],
                rec.vel_m_s[1],
                rec.vel_m_s[2],
            ];
            let state = glonass::propagate(state0, rec.acc_m_s2, tk);
            let clock = glonass::clock_offset_s(rec.clk_bias, rec.gamma_n, tk);
            return Some(([state[0], state[1], state[2]], clock));
        }

        // Supported Keplerian systems only; a record from any other system (e.g.
        // QZSS/SBAS) reports no ephemeris rather than being evaluated with the
        // wrong model. (`from_nav` already restricts records, but `new` accepts
        // arbitrary ones.)
        if !matches!(
            sat.system,
            GnssSystem::Gps | GnssSystem::Galileo | GnssSystem::BeiDou
        ) {
            return None;
        }

        // Map the receive instant (J2000, GPST-aligned) onto the satellite
        // system's continuous time and seconds of week. BeiDou runs on BDT
        // (= GPST - 14 s) with its week epoch 1356 weeks after the GPS epoch, and
        // its geostationary satellites take the GEO orbit branch.
        let gpst_continuous = t_j2000_s + GPS_EPOCH_TO_J2000_S;
        let (t_continuous, is_geo) = if sat.system == GnssSystem::BeiDou {
            (
                gpst_continuous - GPST_MINUS_BDT_S - BDS_EPOCH_MINUS_GPS_EPOCH_S,
                is_beidou_geo(sat),
            )
        } else {
            (gpst_continuous, false)
        };

        let rec = self.select(sat, t_continuous)?;
        let sow = t_continuous.rem_euclid(SECONDS_PER_WEEK);
        let state = satellite_state(
            &rec.elements,
            &rec.clock,
            &rec.constants(),
            sow,
            rec.group_delay_s,
            is_geo,
        );
        Some((
            state.orbit.position().as_array(),
            state.clock.dt_clock_total_s,
        ))
    }
}

#[cfg(test)]
mod tests;