sidereon_core/

terrain.rs

//! DTED tile reader and bilinear terrain lookup.

use std::collections::HashMap;
use std::fs;
use std::path::{Path, PathBuf};

use crate::Error;

pub(crate) const UHL_SIZE: usize = 80;
pub(crate) const DSI_SIZE: usize = 648;
pub(crate) const ACC_SIZE: usize = 2700;
pub(crate) const DATA_OFFSET: usize = UHL_SIZE + DSI_SIZE + ACC_SIZE;
pub(crate) const DATA_SENTINEL: u8 = 0xAA;
pub(crate) const DTED_SUFFIX: &str = concat!("_1arc_v3.d", "t", "2");
const MIN_LOOKUP_LATITUDE_DEG: f64 = -90.0;
const MAX_LOOKUP_LATITUDE_DEG: f64 = 90.0;
const MIN_LOOKUP_LONGITUDE_DEG: f64 = -180.0;
const MAX_LOOKUP_LONGITUDE_DEG: f64 = 180.0;

/// Interpolation mode for DTED terrain lookups.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum DtedInterpolation {
    /// Return the nearest DTED posting as an orthometric height in metres.
    NearestPosting,
    /// Bilinearly interpolate the four surrounding DTED postings as an
    /// orthometric height in metres.
    Bilinear,
}

/// Lookup options for DTED terrain queries.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct DtedLookupOptions {
    /// Interpolation mode used for each orthometric height query.
    pub interpolation: DtedInterpolation,
}

impl Default for DtedLookupOptions {
    fn default() -> Self {
        Self {
            interpolation: DtedInterpolation::Bilinear,
        }
    }
}

/// Lazy DTED terrain reader backed by raw `.dt2` tile bytes.
///
/// Heights returned by this reader are orthometric metres, `H`, above the
/// EGM96 mean sea level geoid used by DTED/SRTM terrain products. They are not
/// ellipsoidal heights above the WGS84 reference ellipsoid.
#[derive(Debug)]
pub struct DtedTerrain {
    root: PathBuf,
    tiles: HashMap<(i32, i32), DtedTile>,
}

impl DtedTerrain {
    /// Build a terrain reader rooted at a directory containing DTED `.dt2`
    /// tiles, either directly or under the repository's block directories.
    #[must_use]
    pub fn new(root: impl Into<PathBuf>) -> Self {
        Self {
            root: root.into(),
            tiles: HashMap::new(),
        }
    }

    /// Return the bilinearly interpolated orthometric height `H` in metres at a
    /// longitude-first geodetic position in degrees.
    pub fn height_m(&mut self, longitude_deg: f64, latitude_deg: f64) -> crate::Result<f64> {
        self.height_m_with_options(longitude_deg, latitude_deg, DtedLookupOptions::default())
    }

    /// Return the orthometric height `H` in metres at a longitude-first
    /// geodetic position in degrees using explicit lookup options.
    pub fn height_m_with_options(
        &mut self,
        longitude_deg: f64,
        latitude_deg: f64,
        options: DtedLookupOptions,
    ) -> crate::Result<f64> {
        validate_lookup_coordinates(longitude_deg, latitude_deg)?;
        let Some(tile) = self.load_tile(longitude_deg, latitude_deg)? else {
            return Ok(0.0);
        };
        height_from_tile(tile, longitude_deg, latitude_deg, options)
    }

    /// Evaluate `(longitude_deg, latitude_deg)` points in order using one
    /// mutable borrow of the resident tile cache.
    ///
    /// The tuple order is intentionally longitude-first, matching
    /// [`Self::height_m_with_options`], even though geoid batch helpers use
    /// latitude-first points.
    pub fn height_batch(
        &mut self,
        points: &[(f64, f64)],
        options: DtedLookupOptions,
    ) -> Vec<crate::Result<f64>> {
        let mut out = Vec::with_capacity(points.len());
        let mut current = None;

        for &(longitude_deg, latitude_deg) in points {
            if let Err(err) = validate_lookup_coordinates(longitude_deg, latitude_deg) {
                out.push(Err(err));
                continue;
            }

            let primary_grid = terrain_grid(longitude_deg, latitude_deg);
            if current == Some(primary_grid) {
                if let Some(tile) = self.tiles.get(&primary_grid) {
                    if tile.contains(longitude_deg, latitude_deg) {
                        out.push(height_from_tile(tile, longitude_deg, latitude_deg, options));
                        continue;
                    }
                }
            }

            match self.resolve_grid(longitude_deg, latitude_deg) {
                Ok(Some(grid_idx)) => {
                    current = Some(grid_idx);
                    let tile = self
                        .tiles
                        .get(&grid_idx)
                        .expect("resolved DTED grid must be present in tile cache");
                    out.push(height_from_tile(tile, longitude_deg, latitude_deg, options));
                }
                Ok(None) => {
                    current = None;
                    out.push(Ok(0.0));
                }
                Err(err) => out.push(Err(err)),
            }
        }

        out
    }

    fn load_tile(&mut self, longitude: f64, latitude: f64) -> crate::Result<Option<&DtedTile>> {
        let Some(grid_idx) = self.resolve_grid(longitude, latitude)? else {
            return Ok(None);
        };
        Ok(self.tiles.get(&grid_idx))
    }

    fn resolve_grid(&mut self, longitude: f64, latitude: f64) -> crate::Result<Option<(i32, i32)>> {
        for grid_idx in terrain_grid_candidates(longitude, latitude) {
            if !self.tiles.contains_key(&grid_idx) {
                let Some(path) = self.terrain_path_for_grid(grid_idx.0, grid_idx.1) else {
                    continue;
                };
                if !path.is_file() {
                    continue;
                }
                let tile = DtedTile::from_path(path).map_err(Error::Parse)?;
                self.tiles.insert(grid_idx, tile);
            }
            if let Some(tile) = self.tiles.get(&grid_idx) {
                if tile.contains(longitude, latitude) {
                    return Ok(Some(grid_idx));
                }
            }
        }
        Ok(None)
    }

    fn terrain_path_for_grid(&self, latitude_index: i32, longitude_index: i32) -> Option<PathBuf> {
        let tile_name = format!(
            "{}_{}{}",
            format_lat(latitude_index),
            format_lon(longitude_index),
            DTED_SUFFIX
        );

        let direct = self.root.join(&tile_name);
        if direct.is_file() {
            return Some(direct);
        }

        let block_dir = terrain_block_dir(latitude_index, longitude_index);
        let nested = self.root.join(&block_dir).join(&tile_name);
        if nested.is_file() {
            return Some(nested);
        }

        let sibling = self.root.parent()?.join(&block_dir).join(&tile_name);
        sibling.is_file().then_some(sibling)
    }
}

fn height_from_tile(
    tile: &DtedTile,
    longitude_deg: f64,
    latitude_deg: f64,
    options: DtedLookupOptions,
) -> crate::Result<f64> {
    if options.interpolation == DtedInterpolation::NearestPosting {
        return tile
            .get_elevation(longitude_deg, latitude_deg)
            .map(|v| v as f64)
            .map_err(Error::Parse);
    }

    let postings_per_deg_lon = tile.lon_count - 1;
    let postings_per_deg_lat = tile.lat_count - 1;

    let lon = scaled_cell_fraction(longitude_deg - tile.origin_longitude, postings_per_deg_lon);
    let lat = scaled_cell_fraction(latitude_deg - tile.origin_latitude, postings_per_deg_lat);
    let lon_lo = lon.cell;
    let lat_lo = lat.cell;
    let fx = lon.fraction;
    let fy = lat.fraction;

    let mut z = 0.0;
    for (di, wx) in [(0i64, 1.0 - fx), (1i64, fx)] {
        for (dj, wy) in [(0i64, 1.0 - fy), (1i64, fy)] {
            let w = wx * wy;
            if w == 0.0 {
                continue;
            }
            let posting_lon =
                tile.origin_longitude + (lon_lo + di) as f64 / postings_per_deg_lon as f64;
            let posting_lat =
                tile.origin_latitude + (lat_lo + dj) as f64 / postings_per_deg_lat as f64;
            z += w * f64::from(
                tile.get_elevation(posting_lon, posting_lat)
                    .map_err(Error::Parse)?,
            );
        }
    }
    Ok(z)
}

pub(crate) fn validate_lookup_coordinates(
    longitude_deg: f64,
    latitude_deg: f64,
) -> crate::Result<()> {
    if !longitude_deg.is_finite() {
        return Err(Error::InvalidInput(
            "longitude_deg must be finite".to_string(),
        ));
    }
    if !latitude_deg.is_finite() {
        return Err(Error::InvalidInput(
            "latitude_deg must be finite".to_string(),
        ));
    }
    if !(MIN_LOOKUP_LONGITUDE_DEG..=MAX_LOOKUP_LONGITUDE_DEG).contains(&longitude_deg) {
        return Err(Error::InvalidInput(
            "longitude_deg must be within [-180, 180]".to_string(),
        ));
    }
    if !(MIN_LOOKUP_LATITUDE_DEG..=MAX_LOOKUP_LATITUDE_DEG).contains(&latitude_deg) {
        return Err(Error::InvalidInput(
            "latitude_deg must be within [-90, 90]".to_string(),
        ));
    }
    Ok(())
}

/// Parsed DTED tile backed by raw `.dt2` bytes.
///
/// Posting values are decoded lazily from DTED signed-magnitude samples.
/// Returned heights are orthometric metres, `H`, above the EGM96 mean sea level
/// geoid.
#[derive(Debug)]
pub struct DtedTile {
    origin_latitude: f64,
    origin_longitude: f64,
    lon_count: usize,
    lat_count: usize,
    data_block_length: usize,
    bytes: Vec<u8>,
}

impl DtedTile {
    /// Read and parse a DTED `.dt2` tile from disk.
    pub fn from_path(path: impl AsRef<Path>) -> Result<Self, String> {
        let bytes =
            fs::read(path.as_ref()).map_err(|e| format!("{}: {e}", path.as_ref().display()))?;
        if bytes.len() < DATA_OFFSET {
            return Err(format!(
                "{} is too short for DTED headers",
                path.as_ref().display()
            ));
        }
        if &bytes[0..4] != b"UHL1" {
            return Err(format!("{} missing UHL1 header", path.as_ref().display()));
        }

        let origin_longitude =
            parse_dted_coord(std::str::from_utf8(&bytes[4..12]).map_err(|e| e.to_string())?)?;
        let origin_latitude =
            parse_dted_coord(std::str::from_utf8(&bytes[12..20]).map_err(|e| e.to_string())?)?;
        let lon_count = parse_ascii_usize(&bytes[47..51])?;
        let lat_count = parse_ascii_usize(&bytes[51..55])?;
        if lon_count < 2 || lat_count < 2 {
            return Err(format!(
                "{} has invalid DTED dimensions lon_count={} lat_count={}; both must be at least 2",
                path.as_ref().display(),
                lon_count,
                lat_count
            ));
        }
        let data_block_length = 12 + 2 * lat_count;
        let expected_len = DATA_OFFSET + lon_count * data_block_length;
        if bytes.len() < expected_len {
            return Err(format!(
                "{} has {} bytes but expected at least {}",
                path.as_ref().display(),
                bytes.len(),
                expected_len
            ));
        }

        Ok(Self {
            origin_latitude,
            origin_longitude,
            lon_count,
            lat_count,
            data_block_length,
            bytes,
        })
    }

    /// Return the nearest orthometric posting value in metres for a
    /// longitude-first geodetic position in degrees.
    pub fn get_elevation(&self, longitude: f64, latitude: f64) -> Result<i16, String> {
        if !self.contains(longitude, latitude) {
            return Err(format!(
                "point ({longitude},{latitude}) is outside DTED tile ({},{})",
                self.origin_longitude, self.origin_latitude
            ));
        }

        let latitude_index =
            nearest_posting_index(latitude - self.origin_latitude, self.lat_count - 1)?;
        let longitude_index =
            nearest_posting_index(longitude - self.origin_longitude, self.lon_count - 1)?;
        if latitude_index >= self.lat_count || longitude_index >= self.lon_count {
            return Err(format!(
                "posting index out of bounds lon={longitude_index} lat={latitude_index}"
            ));
        }

        let block = self.validated_block(longitude_index)?;

        let sample_start = 8 + latitude_index * 2;
        let raw = i16::from_be_bytes([block[sample_start], block[sample_start + 1]]);
        Ok(convert_signed_magnitude(raw))
    }

    pub(crate) fn origin_latitude(&self) -> f64 {
        self.origin_latitude
    }

    pub(crate) fn origin_longitude(&self) -> f64 {
        self.origin_longitude
    }

    pub(crate) fn lon_count(&self) -> usize {
        self.lon_count
    }

    pub(crate) fn lat_count(&self) -> usize {
        self.lat_count
    }

    pub(crate) fn decoded_postings_lon_major(&self) -> Result<Vec<i16>, String> {
        let mut out = Vec::with_capacity(self.lon_count * self.lat_count);
        for longitude_index in 0..self.lon_count {
            let block = self.validated_block(longitude_index)?;
            for latitude_index in 0..self.lat_count {
                let sample_start = 8 + latitude_index * 2;
                let raw = i16::from_be_bytes([block[sample_start], block[sample_start + 1]]);
                out.push(convert_signed_magnitude(raw));
            }
        }
        Ok(out)
    }

    fn contains(&self, longitude: f64, latitude: f64) -> bool {
        latitude >= self.origin_latitude
            && latitude <= self.origin_latitude + 1.0
            && longitude >= self.origin_longitude
            && longitude <= self.origin_longitude + 1.0
    }

    fn validated_block(&self, longitude_index: usize) -> Result<&[u8], String> {
        let block_start = DATA_OFFSET + longitude_index * self.data_block_length;
        let block_end = block_start + self.data_block_length;
        let block = &self.bytes[block_start..block_end];
        if block[0] != DATA_SENTINEL {
            return Err(format!(
                "DTED block {longitude_index} missing data sentinel"
            ));
        }
        let checksum = i32::from_be_bytes([
            block[block.len() - 4],
            block[block.len() - 3],
            block[block.len() - 2],
            block[block.len() - 1],
        ]);
        let sum = block[..block.len() - 4]
            .iter()
            .fold(0i32, |acc, b| acc + i32::from(*b));
        if sum != checksum {
            return Err(format!(
                "DTED checksum failed for block {longitude_index}: expected {checksum}, found {sum}"
            ));
        }
        Ok(block)
    }
}

pub(crate) fn terrain_grid(longitude: f64, latitude: f64) -> (i32, i32) {
    (latitude.floor() as i32, longitude.floor() as i32)
}

pub(crate) fn terrain_grid_candidates(longitude: f64, latitude: f64) -> Vec<(i32, i32)> {
    let (lat, lon) = terrain_grid(longitude, latitude);
    let mut out = vec![(lat, lon)];
    let on_lat_edge = latitude == latitude.floor();
    let on_lon_edge = longitude == longitude.floor();
    if on_lat_edge {
        out.push((lat - 1, lon));
    }
    if on_lon_edge {
        out.push((lat, lon - 1));
    }
    if on_lat_edge && on_lon_edge {
        out.push((lat - 1, lon - 1));
    }
    out
}

pub(crate) fn format_lat(latitude_index: i32) -> String {
    if latitude_index >= 0 {
        format!("n{latitude_index:02}")
    } else {
        format!("s{:02}", -latitude_index)
    }
}

pub(crate) fn format_lon(longitude_index: i32) -> String {
    if longitude_index >= 0 {
        format!("e{longitude_index:03}")
    } else {
        format!("w{:03}", -longitude_index)
    }
}

pub(crate) fn terrain_block_dir(latitude_index: i32, longitude_index: i32) -> String {
    format!(
        "{}_{}",
        format_block_lat(latitude_index),
        format_block_lon(longitude_index)
    )
}

fn format_block_lat(latitude_index: i32) -> String {
    let origin = block_origin(latitude_index);
    if latitude_index >= 0 {
        format!("n{origin:02}")
    } else {
        format!("s{origin:02}")
    }
}

fn format_block_lon(longitude_index: i32) -> String {
    let origin = block_origin(longitude_index);
    if longitude_index >= 0 {
        format!("e{origin:03}")
    } else {
        format!("w{origin:03}")
    }
}

pub(crate) fn block_origin(index: i32) -> u32 {
    (index.unsigned_abs() / 10) * 10
}

fn parse_ascii_usize(bytes: &[u8]) -> Result<usize, String> {
    std::str::from_utf8(bytes)
        .map_err(|e| e.to_string())?
        .trim()
        .parse::<usize>()
        .map_err(|e| e.to_string())
}

fn parse_dted_coord(input: &str) -> Result<f64, String> {
    let hemi = input
        .chars()
        .last()
        .ok_or_else(|| "empty DTED coordinate".to_string())?;
    let sign = match hemi {
        'S' | 'W' => -1.0,
        'N' | 'E' => 1.0,
        _ => return Err(format!("invalid DTED hemisphere {hemi}")),
    };
    let coord = &input[..input.len() - 1];
    let seconds_index = if coord.as_bytes().get(coord.len().saturating_sub(2)) == Some(&b'.') {
        coord.len() - 4
    } else {
        coord.len() - 2
    };
    let minutes_index = seconds_index - 2;
    let degree = coord[..minutes_index]
        .parse::<i32>()
        .map_err(|e| e.to_string())?;
    let minute = coord[minutes_index..seconds_index]
        .parse::<i32>()
        .map_err(|e| e.to_string())?;
    let second = coord[seconds_index..]
        .parse::<f64>()
        .map_err(|e| e.to_string())?;
    Ok(sign * (degree as f64 + ((minute as f64 + second / 60.0) / 60.0)))
}

#[derive(Clone, Copy, Debug)]
pub(crate) struct ScaledCellFraction {
    pub(crate) cell: i64,
    pub(crate) fraction: f64,
    nearest: i64,
}

/// Scale an exact binary64 value by an integer without first rounding their
/// product to binary64.
pub(crate) fn scaled_cell_fraction(offset: f64, postings_per_degree: usize) -> ScaledCellFraction {
    debug_assert!(offset.is_finite());
    debug_assert!(postings_per_degree > 0);

    let bits = offset.to_bits();
    let negative = bits >> 63 != 0;
    let exponent_bits = ((bits >> 52) & 0x7ff) as i32;
    let stored_significand = bits & ((1_u64 << 52) - 1);
    let (significand, exponent) = if exponent_bits == 0 {
        (stored_significand, -1074)
    } else {
        (
            stored_significand | (1_u64 << 52),
            exponent_bits - 1023 - 52,
        )
    };
    if significand == 0 {
        return ScaledCellFraction {
            cell: 0,
            fraction: 0.0,
            nearest: 0,
        };
    }

    let numerator = u128::from(significand) * postings_per_degree as u128;
    if exponent >= 0 {
        let magnitude = numerator
            .checked_shl(exponent as u32)
            .expect("in-tile scaled coordinate must fit u128");
        let magnitude =
            i64::try_from(magnitude).expect("in-tile scaled coordinate must fit a posting index");
        let cell = if negative { -magnitude } else { magnitude };
        return ScaledCellFraction {
            cell,
            fraction: 0.0,
            nearest: cell,
        };
    }

    let denominator_exponent = (-exponent) as u32;
    if denominator_exponent >= 128 {
        if negative {
            return ScaledCellFraction {
                cell: -1,
                fraction: one_minus_dyadic(numerator, denominator_exponent),
                nearest: 0,
            };
        }
        return ScaledCellFraction {
            cell: 0,
            fraction: dyadic_to_f64(numerator, exponent),
            nearest: 0,
        };
    }

    let denominator = 1_u128 << denominator_exponent;
    let integer = numerator >> denominator_exponent;
    let remainder = numerator & (denominator - 1);
    let (cell, euclidean_remainder) = if negative {
        if remainder == 0 {
            (-(integer as i64), 0)
        } else {
            (-(integer as i64) - 1, denominator - remainder)
        }
    } else {
        (integer as i64, remainder)
    };
    let half = denominator >> 1;
    let nearest = if euclidean_remainder < half || (euclidean_remainder == half && cell % 2 == 0) {
        cell
    } else {
        cell + 1
    };
    ScaledCellFraction {
        cell,
        fraction: dyadic_to_f64(euclidean_remainder, exponent),
        nearest,
    }
}

pub(crate) fn nearest_posting_index(
    offset: f64,
    postings_per_degree: usize,
) -> Result<usize, String> {
    let scaled = scaled_cell_fraction(offset, postings_per_degree);
    usize::try_from(scaled.nearest)
        .map_err(|_| format!("cannot round negative posting index {}", scaled.nearest))
}

fn one_minus_dyadic(numerator: u128, denominator_exponent: u32) -> f64 {
    let deficit_units = round_shift_right(numerator, denominator_exponent - 53);
    1.0 - deficit_units as f64 * (f64::EPSILON / 2.0)
}

fn dyadic_to_f64(numerator: u128, exponent: i32) -> f64 {
    if numerator == 0 {
        return 0.0;
    }

    let bit_length = 128 - numerator.leading_zeros();
    let mut binary_exponent = bit_length as i32 - 1 + exponent;
    if binary_exponent >= -1022 {
        let mut significand = if bit_length <= 53 {
            numerator << (53 - bit_length)
        } else {
            round_shift_right(numerator, bit_length - 53)
        };
        if significand == 1_u128 << 53 {
            significand >>= 1;
            binary_exponent += 1;
        }
        let exponent_bits = u64::try_from(binary_exponent + 1023)
            .expect("normal binary64 exponent must be nonnegative");
        let fraction_bits = significand as u64 & ((1_u64 << 52) - 1);
        return f64::from_bits((exponent_bits << 52) | fraction_bits);
    }

    let subnormal_shift = exponent + 1074;
    let significand = if subnormal_shift >= 0 {
        numerator << subnormal_shift as u32
    } else {
        round_shift_right(numerator, (-subnormal_shift) as u32)
    };
    f64::from_bits(significand as u64)
}

fn round_shift_right(value: u128, shift: u32) -> u128 {
    if shift == 0 {
        return value;
    }
    if shift > 128 {
        return 0;
    }

    let quotient = if shift == 128 { 0 } else { value >> shift };
    let remainder = if shift == 128 {
        value
    } else {
        value & ((1_u128 << shift) - 1)
    };
    let half = 1_u128 << (shift - 1);
    if remainder > half || (remainder == half && !quotient.is_multiple_of(2)) {
        quotient + 1
    } else {
        quotient
    }
}

fn convert_signed_magnitude(raw: i16) -> i16 {
    if raw < 0 {
        (-32768i32 - i32::from(raw)) as i16
    } else {
        raw
    }
}

#[cfg(all(test, sidereon_repo_tests))]
mod tests {
    //! DTED batch fixture provenance: adjacent synthetic tiles under
    //! `tests/fixtures/dted/tiles` are written by
    //! `crates/sidereon-core/fixtures-generators/generate_dted_points.py` using
    //! the public DTED UHL/DSI/ACC/data-record layout. Tests compare
    //! `f64::to_bits` exactly, never tolerances.

    use std::fs;
    use std::path::Path;
    use std::path::PathBuf;
    use std::time::{SystemTime, UNIX_EPOCH};

    use serde_json::Value;

    use crate::test_parity::f64_from_hex;
    use crate::Error;

    use super::{
        nearest_posting_index, scaled_cell_fraction, terrain_block_dir, DtedInterpolation,
        DtedLookupOptions, DtedTerrain, DtedTile, DATA_OFFSET, DATA_SENTINEL,
    };

    #[test]
    fn exact_scaling_preserves_the_split_point_fraction() {
        let tile_origin = -107.0;
        let coordinate = -106.265_141_029_846_36;
        let offset = coordinate - tile_origin;
        assert_eq!(offset, 0.734_858_970_153_638_3);

        let scaled = scaled_cell_fraction(offset, 3600);
        let naive = offset * 3600.0;
        let naive_fraction = naive - naive.floor();
        assert_eq!(scaled.cell, 2645);
        assert_eq!(scaled.fraction.to_bits(), 0x3fdf_81b8_9fe7_b000);
        assert_eq!(naive.floor() as i64, 2645);
        assert_eq!(naive_fraction.to_bits(), 0x3fdf_81b8_9fe7_c000);
        assert_eq!(naive_fraction.to_bits() - scaled.fraction.to_bits(), 4096);

        let nondiscriminating_offset = -0.265_141_029_846_361_7;
        let nondiscriminating = scaled_cell_fraction(nondiscriminating_offset, 3600);
        let naive = nondiscriminating_offset * 3600.0;
        assert_eq!(
            nondiscriminating.fraction.to_bits(),
            (naive - naive.floor()).to_bits()
        );
    }

    #[test]
    fn exact_scaling_keeps_a_coordinate_below_the_posting_in_the_lower_cell() {
        let posting = 3.0_f64 / 3600.0;
        let coordinate = f64::from_bits(posting.to_bits() - 1);
        assert!(coordinate < posting);
        assert_eq!(coordinate * 3600.0, 3.0);

        let scaled = scaled_cell_fraction(coordinate, 3600);
        assert_eq!(scaled.cell, 2);
        assert_eq!(scaled.fraction.to_bits(), 0x3fef_ffff_ffff_fffe);
        assert_eq!(scaled.nearest, 3);
    }

    #[test]
    fn nearest_posting_rounds_the_exact_product_instead_of_the_binary64_product() {
        let half_posting = 1.5_f64 / 3600.0;
        let coordinate = f64::from_bits(half_posting.to_bits() - 1);
        assert!(coordinate < half_posting);
        assert_eq!(coordinate * 3600.0, 1.5);

        assert_eq!(nearest_posting_index(coordinate, 3600), Ok(1));
    }

    #[test]
    fn exact_scaled_value_tracks_the_binary64_product_over_deterministic_offsets() {
        let mut state = 0x764e_279d_9f41_2c03_u64;
        for _ in 0..10_000 {
            state = state.wrapping_add(0x9e37_79b9_7f4a_7c15);
            let mut random = state;
            random = (random ^ (random >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
            random = (random ^ (random >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
            random ^= random >> 31;
            let offset = f64::from_bits(0x3ff0_0000_0000_0000 | (random >> 12)) - 1.0;

            let exact = scaled_cell_fraction(offset, 3600);
            let reconstructed = exact.cell as f64 + exact.fraction;
            let naive = offset * 3600.0;
            assert!(
                reconstructed.to_bits().abs_diff(naive.to_bits()) <= 1,
                "offset={offset} exact={reconstructed} naive={naive}"
            );
            assert!(exact.cell <= naive.floor() as i64, "offset={offset}");
        }
    }

    #[test]
    fn terrain_block_dir_matches_reference_bucket_names() {
        assert_eq!(terrain_block_dir(36, -107), "n30_w100");
        assert_eq!(terrain_block_dir(32, -117), "n30_w110");
        assert_eq!(terrain_block_dir(43, -112), "n40_w110");
        assert_eq!(terrain_block_dir(20, -103), "n20_w100");
        assert_eq!(terrain_block_dir(36, 107), "n30_e100");
        assert_eq!(terrain_block_dir(-1, -1), "s00_w000");
        assert_eq!(terrain_block_dir(1, 1), "n00_e000");
        assert_eq!(terrain_block_dir(-1, 1), "s00_e000");
        assert_eq!(terrain_block_dir(32, -110), "n30_w110");
        assert_eq!(terrain_block_dir(32, -111), "n30_w110");
        assert_eq!(terrain_block_dir(32, -1), "n30_w000");
        assert_eq!(terrain_block_dir(32, -10), "n30_w010");
    }

    #[test]
    fn negative_tile_indices_resolve_to_negative_block_dir() {
        let nonce = SystemTime::now()
            .duration_since(UNIX_EPOCH)
            .expect("system time after epoch")
            .as_nanos();
        let root = std::env::temp_dir().join(format!(
            "sidereon-dted-negative-block-{}-{nonce}",
            std::process::id()
        ));
        let tile_dir = root.join("s00_w000");
        let tile_path = tile_dir.join("s01_w001_1arc_v3.dt2");
        fs::create_dir_all(&tile_dir).expect("create nested DTED block dir");
        fs::write(&tile_path, []).expect("create nested DTED tile path");

        let terrain = DtedTerrain::new(&root);
        let got = terrain
            .terrain_path_for_grid(-1, -1)
            .expect("negative nested tile path");
        assert_eq!(got, tile_path);

        fs::remove_dir_all(root).expect("remove temp DTED block dir");
    }

    fn fixture_path(name: &str) -> PathBuf {
        PathBuf::from(env!("CARGO_MANIFEST_DIR"))
            .join("tests")
            .join("fixtures")
            .join("dted")
            .join(name)
    }

    fn bits(v: &Value) -> f64 {
        f64_from_hex(v.as_str().expect("hex-bit string")).expect("valid f64 bits")
    }

    fn temp_path(name: &str) -> PathBuf {
        let nonce = SystemTime::now()
            .duration_since(UNIX_EPOCH)
            .expect("system time after epoch")
            .as_nanos();
        std::env::temp_dir().join(format!("sidereon-{name}-{}-{nonce}", std::process::id()))
    }

    fn scalar_loop(
        root: &Path,
        points: &[(f64, f64)],
        options: DtedLookupOptions,
    ) -> Vec<crate::Result<f64>> {
        let mut terrain = DtedTerrain::new(root);
        points
            .iter()
            .map(|&(lon, lat)| terrain.height_m_with_options(lon, lat, options))
            .collect()
    }

    fn assert_height_results_match(
        got: &[crate::Result<f64>],
        want: &[crate::Result<f64>],
        context: &str,
    ) {
        assert_eq!(got.len(), want.len(), "{context} result length");
        for (idx, (got, want)) in got.iter().zip(want).enumerate() {
            match (got, want) {
                (Ok(got), Ok(want)) => assert_eq!(
                    got.to_bits(),
                    want.to_bits(),
                    "{context} index {idx} height bits"
                ),
                (Err(got), Err(want)) => {
                    assert_eq!(got, want, "{context} index {idx} error")
                }
                (got, want) => panic!("{context} index {idx} mismatch: {got:?} != {want:?}"),
            }
        }
    }

    fn copy_fixture_tile(root: &Path, tile_name: &str) {
        fs::copy(
            fixture_path(&format!("tiles/{tile_name}")),
            root.join(tile_name),
        )
        .expect("copy DTED fixture tile");
    }

    fn copy_primary_fixture_root(name: &str) -> PathBuf {
        let root = temp_path(name);
        fs::create_dir_all(&root).expect("create temp DTED dir");
        copy_fixture_tile(&root, "n36_w107_1arc_v3.dt2");
        root
    }

    fn write_synthetic_dted_tile(
        path: &Path,
        lon_count: usize,
        lat_count: usize,
        sample: impl Fn(usize, usize) -> i16,
    ) {
        let data_block_length = 12 + 2 * lat_count;
        let mut bytes = vec![b' '; DATA_OFFSET];
        bytes[0..4].copy_from_slice(b"UHL1");
        bytes[4..12].copy_from_slice(b"1070000W");
        bytes[12..20].copy_from_slice(b"0360000N");
        bytes[47..51].copy_from_slice(format!("{lon_count:04}").as_bytes());
        bytes[51..55].copy_from_slice(format!("{lat_count:04}").as_bytes());

        for lon_index in 0..lon_count {
            let mut block = vec![0u8; data_block_length];
            block[0] = DATA_SENTINEL;
            for lat_index in 0..lat_count {
                let sample_start = 8 + lat_index * 2;
                block[sample_start..sample_start + 2]
                    .copy_from_slice(&sample(lon_index, lat_index).to_be_bytes());
            }
            let checksum = block[..block.len() - 4]
                .iter()
                .fold(0i32, |acc, b| acc + i32::from(*b));
            let checksum_start = block.len() - 4;
            block[checksum_start..].copy_from_slice(&checksum.to_be_bytes());
            bytes.extend(block);
        }

        fs::write(path, bytes).expect("write synthetic DTED tile");
    }

    #[test]
    fn dted_rejects_degenerate_header_counts() {
        let root = temp_path("dted-degenerate-counts");
        fs::create_dir_all(&root).expect("create temp DTED dir");

        for (lon_count, lat_count) in [(0, 2), (1, 2), (2, 0), (2, 1)] {
            let tile_path = root.join(format!("tile-{lon_count}-{lat_count}.dt2"));
            write_synthetic_dted_tile(&tile_path, lon_count, lat_count, |_, _| 0);

            let err = DtedTile::from_path(&tile_path).expect_err("degenerate counts must error");
            assert!(
                err.contains("invalid DTED dimensions"),
                "unexpected error for lon_count={lon_count} lat_count={lat_count}: {err}"
            );
        }

        fs::remove_dir_all(root).expect("remove temp DTED dir");
    }

    #[test]
    fn dted_lookup_rejects_nonfinite_coordinates() {
        let root = temp_path("dted-nonfinite-coordinates");
        let mut terrain = DtedTerrain::new(&root);

        for (lon, lat, field) in [
            (f64::NAN, 36.5, "longitude_deg"),
            (f64::INFINITY, 36.5, "longitude_deg"),
            (f64::NEG_INFINITY, 36.5, "longitude_deg"),
            (-106.5, f64::NAN, "latitude_deg"),
            (-106.5, f64::INFINITY, "latitude_deg"),
            (-106.5, f64::NEG_INFINITY, "latitude_deg"),
        ] {
            assert_eq!(
                terrain
                    .height_m_with_options(lon, lat, DtedLookupOptions::default())
                    .expect_err("non-finite DTED coordinate must error"),
                Error::InvalidInput(format!("{field} must be finite"))
            );
        }

        assert_eq!(
            terrain
                .height_m(f64::NAN, 36.5)
                .expect_err("height_m must also reject non-finite coordinates"),
            Error::InvalidInput("longitude_deg must be finite".to_string())
        );
    }

    #[test]
    fn dted_lookup_rejects_out_of_range_coordinates() {
        let root = temp_path("dted-out-of-range-coordinates");
        let mut terrain = DtedTerrain::new(&root);

        for (lon, lat, error) in [
            (
                -106.5,
                91.0,
                Error::InvalidInput("latitude_deg must be within [-90, 90]".to_string()),
            ),
            (
                -106.5,
                -90.5,
                Error::InvalidInput("latitude_deg must be within [-90, 90]".to_string()),
            ),
            (
                200.0,
                36.5,
                Error::InvalidInput("longitude_deg must be within [-180, 180]".to_string()),
            ),
            (
                -180.5,
                36.5,
                Error::InvalidInput("longitude_deg must be within [-180, 180]".to_string()),
            ),
        ] {
            assert_eq!(
                terrain
                    .height_m_with_options(lon, lat, DtedLookupOptions::default())
                    .expect_err("out-of-range DTED coordinate must error"),
                error
            );
        }

        assert_eq!(
            terrain
                .height_m(-106.5, 36.5)
                .expect("missing in-range tile keeps sea-level fallback"),
            0.0
        );
    }

    #[test]
    fn dted_valid_minimum_tile_parses_and_interpolates() {
        let root = temp_path("dted-valid-minimum");
        fs::create_dir_all(&root).expect("create temp DTED dir");
        let tile_path = root.join("n36_w107_1arc_v3.dt2");
        write_synthetic_dted_tile(&tile_path, 2, 2, |lon_index, lat_index| {
            match (lon_index, lat_index) {
                (0, 0) => 10,
                (0, 1) => 30,
                (1, 0) => 50,
                (1, 1) => 70,
                _ => unreachable!("2x2 synthetic tile"),
            }
        });

        DtedTile::from_path(&tile_path).expect("valid 2x2 DTED tile");
        let mut terrain = DtedTerrain::new(&root);
        assert_eq!(
            terrain
                .height_m_with_options(
                    -106.5,
                    36.5,
                    DtedLookupOptions {
                        interpolation: DtedInterpolation::Bilinear,
                    },
                )
                .expect("bilinear height"),
            40.0
        );

        fs::remove_dir_all(root).expect("remove temp DTED dir");
    }

    // Fixture provenance: `tests/fixtures/dted/tiles/n36_w107_1arc_v3.dt2` is a
    // synthetic public-format DTED tile written by the committed generator
    // `crates/sidereon-core/fixtures-generators/generate_dted_points.py` using the
    // DTED UHL/DSI/ACC/data-record layout (tile id `n36_w107`, elevation formula
    // `z_m = -20 + 7*lon_i - 5*lat_i + lon_i*lat_i`); no external terrain payload is
    // copied. `tests/fixtures/dted/dted_points.json` holds nearest-posting and
    // bilinear lookup cases generated from that tile. Floating-point fixture
    // values are serialized as f64 hex-bit strings and must be compared with
    // `f64::to_bits`, never tolerances.
    #[test]
    fn dted_lookup_matches_generated_fixture_bits() {
        let raw =
            std::fs::read_to_string(fixture_path("dted_points.json")).expect("read dted fixture");
        let doc: Value = serde_json::from_str(&raw).expect("parse dted fixture");
        assert_eq!(doc["schema"], "gnss-dted-points-v1");

        let root = copy_primary_fixture_root("dted-fixture-single-scalar");
        let mut terrain = DtedTerrain::new(&root);
        let nearest = DtedLookupOptions {
            interpolation: DtedInterpolation::NearestPosting,
        };
        let bilinear = DtedLookupOptions {
            interpolation: DtedInterpolation::Bilinear,
        };

        let mut checked = 0usize;
        for case in doc["nearest_cases"].as_array().expect("nearest_cases") {
            let lon = bits(&case["longitude_bits"]);
            let lat = bits(&case["latitude_bits"]);
            let got = terrain
                .height_m_with_options(lon, lat, nearest)
                .expect("nearest DTED height");
            let want = bits(&case["elevation_bits"]);
            assert_eq!(
                got.to_bits(),
                want.to_bits(),
                "nearest DTED {},{}",
                lon,
                lat
            );
            checked += 1;
        }

        for case in doc["bilinear_cases"].as_array().expect("bilinear_cases") {
            let lon = bits(&case["longitude_bits"]);
            let lat = bits(&case["latitude_bits"]);
            let got = terrain
                .height_m_with_options(lon, lat, bilinear)
                .expect("bilinear DTED height");
            let want = bits(&case["elevation_bits"]);
            assert_eq!(
                got.to_bits(),
                want.to_bits(),
                "bilinear DTED {},{}",
                lon,
                lat
            );
            checked += 1;
        }
        assert!(checked > 0, "empty DTED fixture");
        fs::remove_dir_all(root).expect("remove temp DTED dir");
    }

    #[test]
    fn height_batch_matches_scalar_loop_on_fixture_bits() {
        let raw =
            std::fs::read_to_string(fixture_path("dted_points.json")).expect("read dted fixture");
        let doc: Value = serde_json::from_str(&raw).expect("parse dted fixture");
        assert_eq!(doc["schema"], "gnss-dted-points-v1");

        let points: Vec<(f64, f64)> = ["nearest_cases", "bilinear_cases"]
            .into_iter()
            .flat_map(|cases_key| {
                doc[cases_key]
                    .as_array()
                    .expect(cases_key)
                    .iter()
                    .map(|case| (bits(&case["longitude_bits"]), bits(&case["latitude_bits"])))
            })
            .collect();

        for options in [
            DtedLookupOptions {
                interpolation: DtedInterpolation::NearestPosting,
            },
            DtedLookupOptions {
                interpolation: DtedInterpolation::Bilinear,
            },
        ] {
            let root = copy_primary_fixture_root("dted-fixture-single-batch");
            let want = scalar_loop(&root, &points, options);
            let mut terrain = DtedTerrain::new(&root);
            let got = terrain.height_batch(&points, options);
            assert_height_results_match(&got, &want, "single-tile fixture batch");
            fs::remove_dir_all(root).expect("remove temp DTED dir");
        }
    }

    #[test]
    fn height_batch_matches_scalar_loop_across_adjacent_tiles_bits() {
        let root = fixture_path("tiles");
        let options = DtedLookupOptions {
            interpolation: DtedInterpolation::Bilinear,
        };
        let raw =
            std::fs::read_to_string(fixture_path("dted_points.json")).expect("read dted fixture");
        let doc: Value = serde_json::from_str(&raw).expect("parse dted fixture");
        for case in doc["multi_tile_cases"]
            .as_array()
            .expect("multi_tile_cases")
        {
            let lon = bits(&case["longitude_bits"]);
            let lat = bits(&case["latitude_bits"]);
            let expected = bits(&case["bilinear_bits"]);
            let mut terrain = DtedTerrain::new(&root);
            let got = terrain
                .height_m_with_options(lon, lat, options)
                .expect("multi-tile generated bilinear height");
            assert_eq!(
                got.to_bits(),
                expected.to_bits(),
                "multi-tile generated case {}",
                case["case_id"].as_str().expect("case_id")
            );
        }

        let sequences = [
            (
                "all_in_a_then_all_in_b",
                vec![
                    (-106.875, 36.125),
                    (-106.625, 36.375),
                    (-105.875, 36.125),
                    (-105.625, 36.375),
                ],
            ),
            (
                "interleaved_a_b_a_b",
                vec![
                    (-106.875, 36.625),
                    (-105.875, 36.625),
                    (-106.625, 36.125),
                    (-105.625, 36.125),
                ],
            ),
            (
                "boundary_after_a_then_missing",
                vec![
                    (-106.875, 36.5),
                    (-106.0, 36.5),
                    (-104.5, 36.5),
                    (-105.875, 36.5),
                ],
            ),
        ];

        for (name, points) in sequences {
            let want = scalar_loop(&root, &points, options);
            let mut terrain = DtedTerrain::new(&root);
            let got = terrain.height_batch(&points, options);
            assert_height_results_match(&got, &want, name);
        }

        let mut terrain = DtedTerrain::new(&root);
        let missing = terrain.height_batch(&[(-104.5, 36.5)], options);
        assert_eq!(
            missing[0].as_ref().map(|v| v.to_bits()),
            Ok(0.0f64.to_bits())
        );
    }

    #[test]
    fn height_batch_places_errors_at_input_indices() {
        let root = temp_path("dted-batch-errors");
        fs::create_dir_all(&root).expect("create temp DTED dir");
        copy_fixture_tile(&root, "n36_w107_1arc_v3.dt2");
        copy_fixture_tile(&root, "n36_w106_1arc_v3.dt2");
        fs::write(root.join("n37_w107_1arc_v3.dt2"), b"not a DTED tile")
            .expect("write corrupt DTED tile");

        let points = [
            (-106.875, 36.125),
            (-106.5, f64::NAN),
            (-105.875, 36.125),
            (-106.5, 37.5),
            (-106.625, 36.375),
        ];
        let options = DtedLookupOptions {
            interpolation: DtedInterpolation::Bilinear,
        };
        let want = scalar_loop(&root, &points, options);
        let mut terrain = DtedTerrain::new(&root);
        let got = terrain.height_batch(&points, options);
        assert_height_results_match(&got, &want, "batch error placement");

        assert!(got[0].is_ok(), "index 0 remains valid");
        assert_eq!(
            got[1],
            Err(Error::InvalidInput(
                "latitude_deg must be finite".to_string()
            ))
        );
        assert!(got[2].is_ok(), "index 2 remains valid");
        assert!(
            matches!(&got[3], Err(Error::Parse(msg)) if msg.contains("too short")),
            "index 3 is the corrupt-tile error: {:?}",
            got[3]
        );
        assert!(got[4].is_ok(), "index 4 remains valid");

        fs::remove_dir_all(root).expect("remove temp DTED dir");
    }
}