astrodynamics-gnss 0.9.8

//! RINEX 3.0x observation-file parser and single-frequency pseudorange
//! extraction.
//!
//! Parses a RINEX **version 3** observation file (`OBSERVATION DATA`) into a
//! typed [`RinexObs`] product: the header (including the surveyed
//! [`ObsHeader::approx_position_m`] a-priori receiver position and optional
//! [`ObsHeader::antenna_delta_hen_m`] antenna offset), the per-constellation
//! observation-code table, and the per-epoch
//! satellite→observation values. A pseudorange helper ([`pseudoranges`]) then
//! selects one single-frequency code per system and yields the
//! `(satellite, range_m)` pairs the single-point-positioning solver consumes.
//!
//! # Build vs adopt
//!
//! Like the SP3 and RINEX-NAV readers, this is a hand-rolled, fixed-column text
//! reader in the house style rather than an adoption of the MPL-2.0 `rinex`
//! crate (which would pull a parallel time stack and identifier set into the
//! GNSS layer). The grammar is small and fully specified.
//!
//! It is a **deterministic byte-to-record** parse of a fixed-column text format,
//! not a float recipe; there is no 0-ULP claim here. The pseudorange values are
//! the file's own ASCII decimals parsed to `f64` and carried through unchanged.
//!
//! # Layout (RINEX 3)
//!
//! - Header records are `cols 0..60` content + `cols 60..80` label. The
//!   load-bearing ones are `RINEX VERSION / TYPE` (must be observation, major 3),
//!   `APPROX POSITION XYZ`, `ANTENNA: DELTA H/E/N`, `SYS / # / OBS TYPES` (the
//!   per-system code list, order-preserving, with continuation lines), `TIME OF
//!   FIRST OBS` (+ time system), `INTERVAL`, and the optional
//!   `GLONASS SLOT / FRQ #`.
//! - The body is per-epoch: a `>`-prefixed epoch line carrying the civil time,
//!   an event flag, and the satellite count, then one line per satellite with
//!   each observation as a 16-column `F14.3` value + LLI + SSI field, in the
//!   order the system's `SYS / # / OBS TYPES` list declares.

use std::collections::BTreeMap;

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

use crate::id::{GnssSatelliteId, GnssSystem};
use crate::parse::raw_field as field;
use crate::{Error, Result};

/// Width of one RINEX-3 observation field (`F14.3` value + LLI + SSI).
const OBS_FIELD_WIDTH: usize = 16;
/// Width of the numeric part of one observation field (`F14.3`).
const OBS_VALUE_WIDTH: usize = 14;

/// A civil epoch as it appears on a RINEX observation epoch line, in the file's
/// own time scale (no leap-second shifting). This is the natural boundary for
/// the solver, which derives seconds-of-J2000 / second-of-day / day-of-year
/// from the civil components.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ObsEpochTime {
    /// Four-digit calendar year.
    pub year: i32,
    /// Calendar month, 1..=12.
    pub month: u8,
    /// Calendar day of month, 1..=31.
    pub day: u8,
    /// Hour of day, 0..=23.
    pub hour: u8,
    /// Minute of hour, 0..=59.
    pub minute: u8,
    /// Seconds of minute (fractional), 0.0..60.0.
    pub second: f64,
}

/// One reconstructed observation: a value (or blank) with its loss-of-lock and
/// signal-strength indicators.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ObsValue {
    /// The observed value (meters for code/`C` observables, cycles for `L`,
    /// etc.), or `None` when the field was blank.
    pub value: Option<f64>,
    /// Loss-of-lock indicator (RINEX LLI), `None` when blank.
    pub lli: Option<u8>,
    /// Signal-strength indicator (RINEX SSI), `None` when blank.
    pub ssi: Option<u8>,
}

/// One epoch record: the civil time, the event flag, and the per-satellite
/// observation values (aligned to that system's `SYS / # / OBS TYPES` order).
#[derive(Debug, Clone, PartialEq)]
pub struct ObsEpoch {
    /// Civil epoch in the header time scale.
    pub epoch: ObsEpochTime,
    /// Epoch flag: 0 = OK, 1 = power failure, >1 = an event record (skipped).
    pub flag: u8,
    /// Satellite → observation values, ascending satellite id. The value vector
    /// is index-aligned to [`ObsHeader::obs_codes`] for that satellite's system.
    pub sats: BTreeMap<GnssSatelliteId, Vec<ObsValue>>,
}

/// Parsed RINEX 3 observation header.
#[derive(Debug, Clone, PartialEq)]
pub struct ObsHeader {
    /// The full RINEX version (e.g. `3.05`); the major must be 3.
    pub version: f64,
    /// The surveyed a-priori receiver position (ECEF meters), if the file
    /// carries an `APPROX POSITION XYZ` record.
    pub approx_position_m: Option<[f64; 3]>,
    /// Antenna reference-point offset from the marker in the RINEX
    /// height/east/north convention (meters), if the file carries an
    /// `ANTENNA: DELTA H/E/N` record.
    pub antenna_delta_hen_m: Option<[f64; 3]>,
    /// Per-constellation observation-code list, in declared order.
    pub obs_codes: BTreeMap<GnssSystem, Vec<String>>,
    /// Nominal epoch spacing in seconds (`INTERVAL`), if present.
    pub interval_s: Option<f64>,
    /// First observation epoch and its time system (`TIME OF FIRST OBS`).
    pub time_of_first_obs: Option<(ObsEpochTime, TimeScale)>,
    /// GLONASS slot → frequency channel map (`GLONASS SLOT / FRQ #`), if present.
    pub glonass_slots: BTreeMap<u8, i8>,
    /// Marker (station) name, if present.
    pub marker_name: Option<String>,
}

/// A parsed RINEX 3 observation product.
///
/// Construct with [`RinexObs::parse`]. Epochs are stored in file order; access
/// the header via [`RinexObs::header`], the epochs via [`RinexObs::epochs`], and
/// per-system code lists via [`RinexObs::obs_codes`].
#[derive(Debug, Clone, PartialEq)]
pub struct RinexObs {
    /// The parsed header.
    pub header: ObsHeader,
    /// Epoch records in file order. Event records (flag > 1) are retained with
    /// an empty satellite map so epoch indices stay stable.
    pub epochs: Vec<ObsEpoch>,
}

impl RinexObs {
    /// Parse RINEX 3 observation text into a typed product.
    ///
    /// Returns [`Error::Parse`] if the file is not observation data, is not RINEX
    /// major version 3, is missing a required header record, or has a malformed
    /// epoch record.
    pub fn parse(text: &str) -> Result<Self> {
        let mut parser = Parser::new();
        let mut lines = text.lines();
        parser.parse_header(&mut lines)?;
        parser.parse_body(&mut lines)?;
        parser.finish()
    }

    /// The parsed header.
    pub fn header(&self) -> &ObsHeader {
        &self.header
    }

    /// The epoch records, in file order.
    pub fn epochs(&self) -> &[ObsEpoch] {
        &self.epochs
    }

    /// The observation-code list for a constellation, in declared order.
    pub fn obs_codes(&self, sys: GnssSystem) -> Option<&[String]> {
        self.header.obs_codes.get(&sys).map(Vec::as_slice)
    }
}

impl core::str::FromStr for RinexObs {
    type Err = Error;

    fn from_str(s: &str) -> Result<Self> {
        Self::parse(s)
    }
}

/// Per-system single-frequency code-selection policy.
///
/// For each constellation, an ordered list of observation codes to try; the
/// first one present at an epoch is used. Build the version-aware defaults with
/// [`SignalPolicy::default_for`] and adjust per system with
/// [`SignalPolicy::with_override`].
#[derive(Debug, Clone, PartialEq)]
pub struct SignalPolicy {
    /// Ordered preference list of observation codes per constellation.
    pub codes: BTreeMap<GnssSystem, Vec<String>>,
}

impl SignalPolicy {
    /// The default single-frequency pseudorange policy:
    ///
    /// - GPS `C1C` (L1 C/A),
    /// - Galileo `C1C` then `C1X` (E1),
    /// - BeiDou `C1I` for RINEX 3.02, `C2I` for 3.01 and 3.03+ (the B1I code
    ///   label changed between minor versions),
    /// - GLONASS `C1C` (G1 C/A).
    ///
    /// `version` is the file's RINEX version, which selects the BeiDou default.
    pub fn default_for(version: f64) -> Self {
        let mut codes = BTreeMap::new();
        codes.insert(GnssSystem::Gps, vec!["C1C".to_string()]);
        codes.insert(
            GnssSystem::Galileo,
            vec!["C1C".to_string(), "C1X".to_string()],
        );
        // BeiDou B1I label history: C2I in 3.01, relabelled band 1 (C1I) in
        // 3.02, then reverted to C2I in 3.03 and later. Only the narrow 3.02
        // window prefers C1I; every other version prefers C2I. Offer both, with
        // the version-appropriate one first.
        let beidou = if (3.015..3.025).contains(&version) {
            vec!["C1I".to_string(), "C2I".to_string()]
        } else {
            vec!["C2I".to_string(), "C1I".to_string()]
        };
        codes.insert(GnssSystem::BeiDou, beidou);
        codes.insert(GnssSystem::Glonass, vec!["C1C".to_string()]);
        Self { codes }
    }

    /// Replace the preference list for one constellation.
    pub fn with_override(mut self, sys: GnssSystem, codes: Vec<String>) -> Self {
        self.codes.insert(sys, codes);
        self
    }
}

/// Extract single-frequency pseudoranges for one epoch under a [`SignalPolicy`].
///
/// For each satellite in the epoch, the first code in that system's preference
/// list whose value is present at the epoch is used. Satellites whose system has
/// no policy entry, or that lack every preferred code, are skipped. The result
/// is the ascending-id `(satellite, range_m)` list the solver consumes.
pub fn pseudoranges(
    obs: &RinexObs,
    epoch: &ObsEpoch,
    policy: &SignalPolicy,
) -> Vec<(GnssSatelliteId, f64)> {
    let mut out = Vec::new();
    for (sat, values) in &epoch.sats {
        let Some(prefs) = policy.codes.get(&sat.system) else {
            continue;
        };
        let Some(code_list) = obs.header.obs_codes.get(&sat.system) else {
            continue;
        };
        for code in prefs {
            if let Some(idx) = code_list.iter().position(|c| c == code) {
                if let Some(ObsValue {
                    value: Some(range_m),
                    ..
                }) = values.get(idx)
                {
                    out.push((*sat, *range_m));
                    break;
                }
            }
        }
    }
    out
}

/// Incremental RINEX 3 observation parser state.
struct Parser {
    version: Option<f64>,
    is_observation: bool,
    approx_position_m: Option<[f64; 3]>,
    antenna_delta_hen_m: Option<[f64; 3]>,
    obs_codes: BTreeMap<GnssSystem, Vec<String>>,
    interval_s: Option<f64>,
    time_of_first_obs: Option<(ObsEpochTime, TimeScale)>,
    glonass_slots: BTreeMap<u8, i8>,
    marker_name: Option<String>,
    epochs: Vec<ObsEpoch>,
    /// The constellation whose `SYS / # / OBS TYPES` list is currently being
    /// filled (for continuation lines).
    current_obs_sys: Option<GnssSystem>,
    /// Number of codes still expected for `current_obs_sys`.
    obs_codes_remaining: usize,
}

impl Parser {
    fn new() -> Self {
        Self {
            version: None,
            is_observation: false,
            approx_position_m: None,
            antenna_delta_hen_m: None,
            obs_codes: BTreeMap::new(),
            interval_s: None,
            time_of_first_obs: None,
            glonass_slots: BTreeMap::new(),
            marker_name: None,
            epochs: Vec::new(),
            current_obs_sys: None,
            obs_codes_remaining: 0,
        }
    }

    fn parse_header<'a, I: Iterator<Item = &'a str>>(&mut self, lines: &mut I) -> Result<()> {
        let mut saw_end = false;
        for raw in lines.by_ref() {
            let line = raw.trim_end_matches(['\r', '\n']);
            let label = field(line, 60, 80).trim();
            match label {
                "RINEX VERSION / TYPE" => self.parse_version(line)?,
                "APPROX POSITION XYZ" => self.parse_approx_position(line)?,
                "ANTENNA: DELTA H/E/N" => self.parse_antenna_delta(line)?,
                "SYS / # / OBS TYPES" => self.parse_obs_types(line)?,
                "TIME OF FIRST OBS" => self.parse_time_of_first_obs(line)?,
                "INTERVAL" => {
                    self.interval_s = field(line, 0, 10).trim().parse::<f64>().ok();
                }
                "GLONASS SLOT / FRQ #" => self.parse_glonass_slots(line),
                "MARKER NAME" => {
                    let name = field(line, 0, 60).trim();
                    if !name.is_empty() {
                        self.marker_name = Some(name.to_string());
                    }
                }
                "END OF HEADER" => {
                    saw_end = true;
                    break;
                }
                // Every other header record is tolerated and skipped.
                _ => {}
            }
        }
        if !saw_end {
            return Err(Error::Parse("RINEX OBS header has no END OF HEADER".into()));
        }
        Ok(())
    }

    fn parse_version(&mut self, line: &str) -> Result<()> {
        let version = field(line, 0, 20)
            .trim()
            .parse::<f64>()
            .map_err(|_| Error::Parse(format!("RINEX OBS version unparsable in {line:?}")))?;
        // The file type letter is at column 20; observation files carry 'O'.
        let type_field = field(line, 20, 40);
        self.is_observation =
            type_field.trim_start().starts_with('O') || type_field.contains("OBSERVATION");
        if !self.is_observation {
            return Err(Error::Parse(format!(
                "RINEX file is not observation data: {type_field:?}"
            )));
        }
        if version.floor() as i64 != 3 {
            return Err(Error::Parse(format!(
                "RINEX OBS parser requires major version 3, got {version}"
            )));
        }
        self.version = Some(version);
        Ok(())
    }

    fn parse_approx_position(&mut self, line: &str) -> Result<()> {
        let body = field(line, 0, 60);
        let parts: Vec<f64> = body
            .split_whitespace()
            .filter_map(|t| t.parse::<f64>().ok())
            .collect();
        if parts.len() >= 3 {
            self.approx_position_m = Some([parts[0], parts[1], parts[2]]);
        }
        Ok(())
    }

    fn parse_antenna_delta(&mut self, line: &str) -> Result<()> {
        let body = field(line, 0, 60);
        let parts: Vec<f64> = body
            .split_whitespace()
            .filter_map(|t| t.parse::<f64>().ok())
            .collect();
        if parts.len() >= 3 {
            self.antenna_delta_hen_m = Some([parts[0], parts[1], parts[2]]);
        }
        Ok(())
    }

    fn parse_obs_types(&mut self, line: &str) -> Result<()> {
        // A new system line carries its letter at column 0 and the count at
        // columns 3..6; a continuation line has a blank system field and only
        // adds more codes to the current system.
        let sys_field = field(line, 0, 1).trim();
        if !sys_field.is_empty() {
            let letter = sys_field.chars().next().unwrap();
            let system = GnssSystem::from_letter(letter).ok_or_else(|| {
                Error::Parse(format!("RINEX OBS unknown system letter {letter:?}"))
            })?;
            let count = field(line, 3, 6).trim().parse::<usize>().map_err(|_| {
                Error::Parse(format!("RINEX OBS obs-type count unparsable in {line:?}"))
            })?;
            self.current_obs_sys = Some(system);
            self.obs_codes_remaining = count;
            self.obs_codes.entry(system).or_default();
        }
        let Some(system) = self.current_obs_sys else {
            return Ok(());
        };
        // Codes occupy 4-wide fields (" CCC") from column 7; collect up to the
        // remaining count.
        let codes_section = field(line, 7, 60);
        let list = self.obs_codes.get_mut(&system).expect("system inserted");
        for tok in codes_section.split_whitespace() {
            if self.obs_codes_remaining == 0 {
                break;
            }
            list.push(tok.to_string());
            self.obs_codes_remaining -= 1;
        }
        Ok(())
    }

    fn parse_time_of_first_obs(&mut self, line: &str) -> Result<()> {
        let body = field(line, 0, 43);
        let nums: Vec<f64> = body
            .split_whitespace()
            .filter_map(|t| t.parse::<f64>().ok())
            .collect();
        if nums.len() >= 6 {
            let epoch = ObsEpochTime {
                year: nums[0] as i32,
                month: nums[1] as u8,
                day: nums[2] as u8,
                hour: nums[3] as u8,
                minute: nums[4] as u8,
                second: nums[5],
            };
            let scale_label = field(line, 48, 51).trim();
            let scale = time_scale_from_label(scale_label);
            self.time_of_first_obs = Some((epoch, scale));
        }
        Ok(())
    }

    fn parse_glonass_slots(&mut self, line: &str) {
        // " N R01  1 R02 -4 ...": a count then 7-wide "SVNN ±k" entries.
        let body = field(line, 4, 60);
        let tokens: Vec<&str> = body.split_whitespace().collect();
        let mut i = 0;
        while i + 1 < tokens.len() {
            let sv = tokens[i];
            let ch = tokens[i + 1];
            if let Some(rest) = sv.strip_prefix('R') {
                if let (Ok(slot), Ok(channel)) = (rest.parse::<u8>(), ch.parse::<i8>()) {
                    self.glonass_slots.insert(slot, channel);
                }
            }
            i += 2;
        }
    }

    fn parse_body<'a, I: Iterator<Item = &'a str>>(&mut self, lines: &mut I) -> Result<()> {
        while let Some(raw) = lines.next() {
            let line = raw.trim_end_matches(['\r', '\n']);
            if line.is_empty() {
                continue;
            }
            if !line.starts_with('>') {
                // A stray non-epoch line outside an epoch block; tolerate.
                continue;
            }
            let (epoch_time, flag, numsat) = parse_epoch_line(line)?;

            if flag > 1 {
                // Event record: the next `numsat` lines are header/comment
                // records, not observations. Consume and skip them, keeping a
                // placeholder epoch so indices stay meaningful.
                for _ in 0..numsat {
                    let _ = lines.next();
                }
                self.epochs.push(ObsEpoch {
                    epoch: epoch_time,
                    flag,
                    sats: BTreeMap::new(),
                });
                continue;
            }

            let mut sats = BTreeMap::new();
            for _ in 0..numsat {
                let sat_line = lines.next().ok_or_else(|| {
                    Error::Parse("RINEX OBS epoch truncated: missing satellite line".into())
                })?;
                let sat_line = sat_line.trim_end_matches(['\r', '\n']);
                let (sat, values) = self.parse_sat_line(sat_line)?;
                sats.insert(sat, values);
            }
            self.epochs.push(ObsEpoch {
                epoch: epoch_time,
                flag,
                sats,
            });
        }
        Ok(())
    }

    fn parse_sat_line(&self, line: &str) -> Result<(GnssSatelliteId, Vec<ObsValue>)> {
        let token = field(line, 0, 3);
        let sat = parse_sv_token(token).ok_or_else(|| {
            Error::Parse(format!("RINEX OBS unparsable satellite token {token:?}"))
        })?;
        let n_obs = self.obs_codes.get(&sat.system).map(Vec::len).unwrap_or(0);
        let mut values = Vec::with_capacity(n_obs);
        for i in 0..n_obs {
            let start = 3 + i * OBS_FIELD_WIDTH;
            let value_str = field(line, start, start + OBS_VALUE_WIDTH).trim();
            let value = if value_str.is_empty() {
                None
            } else {
                Some(value_str.parse::<f64>().map_err(|_| {
                    Error::Parse(format!(
                        "RINEX OBS observation {value_str:?} unparsable on {line:?}"
                    ))
                })?)
            };
            let lli = digit_at(line, start + OBS_VALUE_WIDTH);
            let ssi = digit_at(line, start + OBS_VALUE_WIDTH + 1);
            values.push(ObsValue { value, lli, ssi });
        }
        Ok((sat, values))
    }

    fn finish(self) -> Result<RinexObs> {
        let version = self
            .version
            .ok_or_else(|| Error::Parse("RINEX OBS missing RINEX VERSION / TYPE".into()))?;
        if self.obs_codes.is_empty() {
            return Err(Error::Parse(
                "RINEX OBS header has no SYS / # / OBS TYPES records".into(),
            ));
        }
        let header = ObsHeader {
            version,
            approx_position_m: self.approx_position_m,
            antenna_delta_hen_m: self.antenna_delta_hen_m,
            obs_codes: self.obs_codes,
            interval_s: self.interval_s,
            time_of_first_obs: self.time_of_first_obs,
            glonass_slots: self.glonass_slots,
            marker_name: self.marker_name,
        };
        Ok(RinexObs {
            header,
            epochs: self.epochs,
        })
    }
}

/// Parse a RINEX-3 epoch line `> YYYY MM DD HH MM SS.sssssss  F NN [clock]`,
/// returning the civil time, event flag, and satellite count.
fn parse_epoch_line(line: &str) -> Result<(ObsEpochTime, u8, usize)> {
    // The date occupies a fixed width after the leading '>'; the flag is at
    // column 31 and the satellite count at columns 32..35.
    let date_body = field(line, 1, 29);
    let nums: Vec<f64> = date_body
        .split_whitespace()
        .filter_map(|t| t.parse::<f64>().ok())
        .collect();
    if nums.len() < 6 {
        return Err(Error::Parse(format!(
            "RINEX OBS epoch line has too few date fields: {line:?}"
        )));
    }
    let epoch = ObsEpochTime {
        year: nums[0] as i32,
        month: nums[1] as u8,
        day: nums[2] as u8,
        hour: nums[3] as u8,
        minute: nums[4] as u8,
        second: nums[5],
    };
    let flag = field(line, 31, 32).trim().parse::<u8>().unwrap_or(0);
    let numsat = field(line, 32, 35).trim().parse::<usize>().map_err(|_| {
        Error::Parse(format!(
            "RINEX OBS epoch satellite count unparsable: {line:?}"
        ))
    })?;
    Ok((epoch, flag, numsat))
}

/// Map a RINEX time-system label onto the core [`TimeScale`]. An empty/unknown
/// label defaults to GPS time, which is the scale a multi-GNSS observation file
/// uses in practice.
fn time_scale_from_label(label: &str) -> TimeScale {
    match label.trim() {
        "GPS" => TimeScale::Gpst,
        "GAL" => TimeScale::Gst,
        "BDT" => TimeScale::Bdt,
        "UTC" => TimeScale::Utc,
        "TAI" => TimeScale::Tai,
        _ => TimeScale::Gpst,
    }
}

/// Parse a 3-char SV token (e.g. `G01`, `C30`) into a [`GnssSatelliteId`].
fn parse_sv_token(token: &str) -> Option<GnssSatelliteId> {
    let token = token.trim();
    let first = token.chars().next()?;
    let system = GnssSystem::from_letter(first)?;
    let prn = token[first.len_utf8()..].trim().parse::<u8>().ok()?;
    Some(GnssSatelliteId::new(system, prn))
}

/// Read a single decimal digit at byte `col`, or `None` if it is blank /
/// non-digit / past end of line.
fn digit_at(line: &str, col: usize) -> Option<u8> {
    line.as_bytes()
        .get(col)
        .filter(|b| b.is_ascii_digit())
        .map(|b| b - b'0')
}

#[cfg(test)]
mod tests;