astrodynamics-gnss 0.9.7

//! IONEX vertical-TEC grid parser.
//!
//! Parses an IONEX (IONosphere map EXchange) ASCII product into the vertical-TEC
//! grids and the geometry needed to interpolate them: the latitude and longitude
//! node axes, the per-map TEC and (optionally) RMS grids, and the map epochs on a
//! J2000-second axis. The float math that turns the grid into a slant delay lives
//! in [`super::slant`]; this module is the deterministic byte/record reader.
//!
//! IONEX stores TEC as a scaled integer field: the printed value times
//! `10^EXPONENT`. Latitude bands are written north-to-south (a negative `DLAT`),
//! longitude nodes west-to-east (a positive `DLON`), and the longitude span
//! includes the `+180` wrap-seam column. The node axes are reconstructed as
//! `v1 + i * step` so they match the producer's `arange`-style construction
//! bit-for-bit, and the TEC scaling is `field * 10^EXPONENT` formed as a single
//! multiply.

use crate::error::{Error, Result};

/// A parsed IONEX vertical-TEC product.
///
/// The grids are indexed `[map][i_lat][i_lon]`, with `lat_nodes_deg` descending
/// (north-to-south) and `lon_nodes_deg` ascending (west-to-east). TEC and RMS are
/// in TECU after the `10^EXPONENT` scaling. Map epochs are J2000 seconds.
#[derive(Debug, Clone, PartialEq)]
pub struct Ionex {
    /// Latitude node values in degrees, descending (north-to-south).
    lat_nodes_deg: Vec<f64>,
    /// Longitude node values in degrees, ascending (west-to-east).
    lon_nodes_deg: Vec<f64>,
    /// Signed latitude step in degrees (negative for the standard ordering).
    dlat_deg: f64,
    /// Signed longitude step in degrees (positive for the standard ordering).
    dlon_deg: f64,
    /// Single-layer shell height in kilometers.
    shell_height_km: f64,
    /// Mean earth radius used by the geometry, in kilometers.
    base_radius_km: f64,
    /// The integer `EXPONENT` header field (TEC scale is `10^EXPONENT`).
    exponent: i32,
    /// Map epochs on the J2000-second axis, ascending.
    map_epochs_s: Vec<i64>,
    /// Per-map vertical-TEC grids, indexed `[map][i_lat][i_lon]` (TECU).
    tec_maps: Vec<Vec<Vec<f64>>>,
    /// Per-map RMS grids, indexed `[map][i_lat][i_lon]` (TECU); empty if absent.
    rms_maps: Vec<Vec<Vec<f64>>>,
}

impl Ionex {
    /// Latitude node values in degrees (descending, north-to-south).
    pub fn lat_nodes_deg(&self) -> &[f64] {
        &self.lat_nodes_deg
    }

    /// Longitude node values in degrees (ascending, west-to-east).
    pub fn lon_nodes_deg(&self) -> &[f64] {
        &self.lon_nodes_deg
    }

    /// Signed latitude step in degrees (negative for the standard ordering).
    pub fn dlat_deg(&self) -> f64 {
        self.dlat_deg
    }

    /// Signed longitude step in degrees.
    pub fn dlon_deg(&self) -> f64 {
        self.dlon_deg
    }

    /// Single-layer shell height in kilometers.
    pub fn shell_height_km(&self) -> f64 {
        self.shell_height_km
    }

    /// Mean earth radius used by the geometry, in kilometers.
    pub fn base_radius_km(&self) -> f64 {
        self.base_radius_km
    }

    /// The IONEX `EXPONENT` header field; the TEC scale is `10^EXPONENT`.
    pub fn exponent(&self) -> i32 {
        self.exponent
    }

    /// Map epochs on the J2000-second axis (ascending).
    pub fn map_epochs_s(&self) -> &[i64] {
        &self.map_epochs_s
    }

    /// Per-map vertical-TEC grids, indexed `[map][i_lat][i_lon]` (TECU).
    pub fn tec_maps(&self) -> &[Vec<Vec<f64>>] {
        &self.tec_maps
    }

    /// Per-map RMS grids, indexed `[map][i_lat][i_lon]` (TECU); empty if the
    /// product carries no RMS maps.
    pub fn rms_maps(&self) -> &[Vec<Vec<f64>>] {
        &self.rms_maps
    }

    /// Parse an IONEX product from its bytes.
    pub fn parse(bytes: &[u8]) -> Result<Self> {
        let text = core::str::from_utf8(bytes)
            .map_err(|_| Error::Parse("IONEX is not valid UTF-8".into()))?;
        Self::parse_str(text)
    }

    /// Parse an IONEX product from its text.
    pub fn parse_str(text: &str) -> Result<Self> {
        let mut lines = text.lines();

        // ---- Header ----
        let mut lat1 = None;
        let mut lat2 = None;
        let mut dlat = None;
        let mut lon1 = None;
        let mut lon2 = None;
        let mut dlon = None;
        let mut shell_height_km = None;
        let mut base_radius_km = None;
        let mut exponent: i32 = -1; // IONEX default when no EXPONENT record.

        for line in lines.by_ref() {
            let label = label_of(line);
            match label {
                "LAT1 / LAT2 / DLAT" => {
                    let (a, b, c) = three_fields(line, "LAT1 / LAT2 / DLAT")?;
                    lat1 = Some(a);
                    lat2 = Some(b);
                    dlat = Some(c);
                }
                "LON1 / LON2 / DLON" => {
                    let (a, b, c) = three_fields(line, "LON1 / LON2 / DLON")?;
                    lon1 = Some(a);
                    lon2 = Some(b);
                    dlon = Some(c);
                }
                "HGT1 / HGT2 / DHGT" => {
                    let (a, _b, _c) = three_fields(line, "HGT1 / HGT2 / DHGT")?;
                    shell_height_km = Some(a);
                }
                "BASE RADIUS" => {
                    base_radius_km = Some(first_field(line, "BASE RADIUS")?);
                }
                "EXPONENT" => {
                    exponent = first_field_int(line, "EXPONENT")?;
                }
                "END OF HEADER" => break,
                _ => {}
            }
        }

        let lat1 = lat1.ok_or_else(|| Error::Parse("IONEX missing LAT1 / LAT2 / DLAT".into()))?;
        let lat2 = lat2.unwrap();
        let dlat = dlat.unwrap();
        let lon1 = lon1.ok_or_else(|| Error::Parse("IONEX missing LON1 / LON2 / DLON".into()))?;
        let lon2 = lon2.unwrap();
        let dlon = dlon.unwrap();
        let shell_height_km = shell_height_km
            .ok_or_else(|| Error::Parse("IONEX missing HGT1 / HGT2 / DHGT".into()))?;
        let base_radius_km =
            base_radius_km.ok_or_else(|| Error::Parse("IONEX missing BASE RADIUS".into()))?;

        let lat_nodes_deg = node_axis(lat1, lat2, dlat)?;
        let lon_nodes_deg = node_axis(lon1, lon2, dlon)?;
        let nlat = lat_nodes_deg.len();
        let nlon = lon_nodes_deg.len();

        // TEC scale `10^EXPONENT`, formed as a single multiply per field.
        let scale = 10f64.powi(exponent);

        // ---- Body ----
        let mut map_epochs_s = Vec::new();
        let mut tec_maps: Vec<Vec<Vec<f64>>> = Vec::new();
        let mut rms_maps: Vec<Vec<Vec<f64>>> = Vec::new();

        // The reader is line-driven: a START record opens a map, an EPOCH record
        // sets its time, each LAT/LON1/LON2/DLON/H record opens a latitude band
        // whose TEC values are read from the following continuation lines.
        let mut cur_kind: Option<MapKind> = None;
        let mut cur_grid: Vec<Vec<f64>> = Vec::new();
        let mut cur_epoch: Option<i64> = None;
        let mut band_vals: Vec<f64> = Vec::new();
        let mut reading_band = false;

        for line in lines {
            let label = label_of(line);
            match label {
                "START OF TEC MAP" => {
                    cur_kind = Some(MapKind::Tec);
                    cur_grid = Vec::with_capacity(nlat);
                    cur_epoch = None;
                }
                "START OF RMS MAP" => {
                    cur_kind = Some(MapKind::Rms);
                    cur_grid = Vec::with_capacity(nlat);
                    cur_epoch = None;
                }
                "EPOCH OF CURRENT MAP" => {
                    cur_epoch = Some(parse_epoch_j2000_s(line)?);
                }
                "LAT/LON1/LON2/DLON/H" => {
                    if reading_band {
                        finish_band(&mut cur_grid, &mut band_vals, nlon)?;
                    }
                    reading_band = true;
                    band_vals = Vec::with_capacity(nlon);
                }
                "END OF TEC MAP" | "END OF RMS MAP" => {
                    if reading_band {
                        finish_band(&mut cur_grid, &mut band_vals, nlon)?;
                        reading_band = false;
                    }
                    if cur_grid.len() != nlat {
                        return Err(Error::Parse(format!(
                            "IONEX map has {} latitude bands, expected {nlat}",
                            cur_grid.len()
                        )));
                    }
                    match cur_kind {
                        Some(MapKind::Tec) => {
                            let ep = cur_epoch.ok_or_else(|| {
                                Error::Parse("IONEX TEC map missing EPOCH OF CURRENT MAP".into())
                            })?;
                            map_epochs_s.push(ep);
                            tec_maps.push(core::mem::take(&mut cur_grid));
                        }
                        Some(MapKind::Rms) => {
                            rms_maps.push(core::mem::take(&mut cur_grid));
                        }
                        None => {
                            return Err(Error::Parse("IONEX END OF MAP without START".into()));
                        }
                    }
                    cur_kind = None;
                }
                _ => {
                    // A continuation line holding the TEC/RMS integer fields of
                    // the current latitude band, or an unrecognized header line.
                    if reading_band {
                        for tok in line.split_whitespace() {
                            let v: i64 = tok.parse().map_err(|_| {
                                Error::Parse(format!("IONEX TEC field unparsable: {tok:?}"))
                            })?;
                            band_vals.push(v as f64 * scale);
                        }
                    }
                }
            }
        }

        if tec_maps.is_empty() {
            return Err(Error::Parse("IONEX has no TEC maps".into()));
        }
        // Bilinear interpolation brackets a cell with `node[i+1]` / `node[j+1]`,
        // so each axis needs at least two nodes. Reject a degenerate grid here
        // rather than letting evaluation index past the end.
        if lat_nodes_deg.len() < 2 || lon_nodes_deg.len() < 2 {
            return Err(Error::Parse(format!(
                "IONEX grid needs >= 2 nodes per axis (got {} lat, {} lon)",
                lat_nodes_deg.len(),
                lon_nodes_deg.len()
            )));
        }
        if !rms_maps.is_empty() && rms_maps.len() != tec_maps.len() {
            return Err(Error::Parse(
                "IONEX RMS map count does not match TEC map count".into(),
            ));
        }

        Ok(Self {
            lat_nodes_deg,
            lon_nodes_deg,
            dlat_deg: dlat,
            dlon_deg: dlon,
            shell_height_km,
            base_radius_km,
            exponent,
            map_epochs_s,
            tec_maps,
            rms_maps,
        })
    }
}

/// Whether the current map being read is a TEC or an RMS map.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum MapKind {
    Tec,
    Rms,
}

/// Move an accumulated latitude band into the current grid, validating width.
fn finish_band(grid: &mut Vec<Vec<f64>>, band: &mut Vec<f64>, nlon: usize) -> Result<()> {
    if band.len() != nlon {
        return Err(Error::Parse(format!(
            "IONEX latitude band has {} values, expected {nlon}",
            band.len()
        )));
    }
    grid.push(core::mem::take(band));
    Ok(())
}

/// The 20-character label field of an IONEX record (columns 60..80), trimmed.
fn label_of(line: &str) -> &str {
    if line.len() <= 60 {
        line.trim()
    } else {
        line[60..].trim_end().trim_start()
    }
}

/// Parse the first whitespace-delimited float of the data portion of a record.
fn first_field(line: &str, label: &str) -> Result<f64> {
    let data = data_of(line);
    data.split_whitespace()
        .next()
        .ok_or_else(|| Error::Parse(format!("IONEX {label} record empty")))?
        .parse()
        .map_err(|_| Error::Parse(format!("IONEX {label} field unparsable")))
}

/// Parse the first whitespace-delimited integer of the data portion.
fn first_field_int(line: &str, label: &str) -> Result<i32> {
    let data = data_of(line);
    data.split_whitespace()
        .next()
        .ok_or_else(|| Error::Parse(format!("IONEX {label} record empty")))?
        .parse()
        .map_err(|_| Error::Parse(format!("IONEX {label} field unparsable")))
}

/// Parse the first three whitespace-delimited floats of the data portion.
fn three_fields(line: &str, label: &str) -> Result<(f64, f64, f64)> {
    let data = data_of(line);
    let mut it = data.split_whitespace();
    let a = next_float(&mut it, label)?;
    let b = next_float(&mut it, label)?;
    let c = next_float(&mut it, label)?;
    Ok((a, b, c))
}

fn next_float<'a>(it: &mut impl Iterator<Item = &'a str>, label: &str) -> Result<f64> {
    it.next()
        .ok_or_else(|| Error::Parse(format!("IONEX {label} record short")))?
        .parse()
        .map_err(|_| Error::Parse(format!("IONEX {label} field unparsable")))
}

/// The data portion of a record (columns 0..60), or the whole short line.
fn data_of(line: &str) -> &str {
    if line.len() <= 60 {
        line
    } else {
        &line[..60]
    }
}

/// Build a node axis `v1 + i * step`, with the count taken from the inclusive
/// `[v1, v2]` span (a half-step guard on the end matches the producer's
/// `arange`-style construction).
fn node_axis(v1: f64, v2: f64, step: f64) -> Result<Vec<f64>> {
    if step == 0.0 {
        return Err(Error::Parse("IONEX grid step is zero".into()));
    }
    let guard = if step > 0.0 { 0.1 } else { -0.1 };
    let span = (v2 + guard - v1) / step;
    if span < 0.0 {
        return Err(Error::Parse("IONEX grid span has the wrong sign".into()));
    }
    let n = span.floor() as usize + 1;
    Ok((0..n).map(|i| v1 + (i as f64) * step).collect())
}

/// Parse an `EPOCH OF CURRENT MAP` record into J2000 seconds.
///
/// The record carries `year month day hour minute second`; the result is the
/// integer number of seconds from the J2000 epoch (2000-01-01 12:00:00).
fn parse_epoch_j2000_s(line: &str) -> Result<i64> {
    let data = data_of(line);
    let mut it = data.split_whitespace();
    let mut next_int = |what: &str| -> Result<i64> {
        it.next()
            .ok_or_else(|| Error::Parse(format!("IONEX epoch missing {what}")))?
            .parse::<i64>()
            .map_err(|_| Error::Parse(format!("IONEX epoch {what} unparsable")))
    };
    let year = next_int("year")?;
    let month = next_int("month")?;
    let day = next_int("day")?;
    let hour = next_int("hour")?;
    let minute = next_int("minute")?;
    let second = next_int("second")?;

    let days = days_from_civil(year, month, day) - days_from_civil(2000, 1, 1);
    // J2000 origin is noon; subtract 12 h from the midnight-based day count.
    Ok(days * 86_400 - 43_200 + hour * 3_600 + minute * 60 + second)
}

/// Days from the civil date to 1970-01-01 (Howard Hinnant's algorithm), used
/// only as a difference so the absolute origin cancels.
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 mp = m + if m > 2 { -3 } else { 9 };
    let doy = (153 * mp + 2) / 5 + d - 1;
    let doe = yoe * 365 + yoe / 4 - yoe / 100 + doy;
    era * 146_097 + doe - 719_468
}