peteksim 0.1.1

//! `Project` — the ingest front-door. Walks a Petrel-export tree, routes each
//! file to the owning petekio loader (extension-dispatch does the per-file
//! work), and records a loud [`Inventory`] of what loaded and what was
//! skipped-with-reason. No parsing lives here; classification + routing only.

use super::spec::TrendSpec;
use crate::units::SrsError;
use petekio::{
    BBox, GeoData, GridGeometry, GridMethod, Interval, LogView, NameMap, Point3, Sidetrack,
    Surface, Unit,
};
use std::collections::BTreeSet;
use std::path::{Path, PathBuf};

mod helpers;
use helpers::*;

/// What a [`Project::load`] found — every loaded artifact by kind, plus the
/// files it skipped and why (never silent).
#[derive(Debug, Clone, Default)]
pub struct Inventory {
    pub surfaces: Vec<String>,
    pub polygons: Vec<String>,
    pub points: Vec<String>,
    /// Loaded wells as **bore-level** ids — one entry per positioned bore under
    /// its bore-qualified id (matching [`Project::wells`]), not the parent well id.
    pub wells: Vec<String>,
    /// Named tops-pick surfaces available across the loaded wells (e.g. "GOC").
    pub tops: Vec<String>,
    /// Core LAS files folded into a bore by petekio's well loader (the LAS 3.0
    /// core-file merge): `(filename, bore_id)`. Without this the merge leaves no
    /// inventory trace — the core log is absorbed into its bore and the file
    /// vanishes (F8). `bore_id` names the bore it merged into (`"<id> <bore>"`,
    /// or the well id for a main-bore merge).
    pub merged: Vec<(String, String)>,
    /// `(path, reason)` for every file the loader could not place.
    pub skipped: Vec<(String, String)>,
}

/// A named tops surface's per-well picks and a representative level.
#[derive(Debug, Clone)]
pub struct TopsPick {
    pub name: String,
    /// `(well_id, depth_m)` — subsea TVD (positive-down) at each well's pick.
    pub picks: Vec<(String, f64)>,
    /// Mean pick depth \[m\] — the representative contact level.
    pub level_m: f64,
    /// Sample standard deviation of the picks \[m\] (0 for a single pick).
    pub spread_m: f64,
}

/// One **bore** of a loaded well, presented as an independent positioned "well"
/// for modelling (the R-a wiring). petekio makes every bore of a multi-sidetrack
/// well a first-class [`Sidetrack`] with its own trajectory + logs + tops; this
/// wrapper pairs a bore with its **bore-qualified id** (`"<id> <bore>"`, or the
/// plain well id for the main bore) — the same key petekio's `model_inputs`
/// emits — so the framework tie + property upscaling consume each bore's *own*
/// positions. We never call the well-level accessors on a multi-bore well (they
/// raise / resolve-empty by design — that raise is the regression signal, not a
/// thing to swallow); we go straight to the bore.
pub struct BoreWell<'a> {
    /// Bore-qualified id: `"<well id> <bore>"` (or the plain well id for the
    /// main bore) — matches `model_inputs().spatial.well_curves[..].well_id`.
    pub id: String,
    bore: &'a Sidetrack,
}

impl BoreWell<'_> {
    /// The interval named by top `name` on this bore (its own picks), or `None`.
    pub fn top(&self, name: &str) -> Option<Interval<'_>> {
        self.bore.top(name)
    }

    /// Interpolated world position at `md` on this bore's own trajectory
    /// (`z` = negative-down elevation), or `None` if unpositioned.
    pub fn xyz(&self, md: f64) -> Option<Point3> {
        self.bore.xyz(md)
    }

    /// The bore's positioned path as world `[x, y, tvd]` stations (tvd
    /// **positive-down**, `tvd == -xyz.z`), sampled evenly along its active
    /// survey's MD range. Empty when the bore carries no trajectory. Feeds the
    /// viewer's well markers + along-bore section (`intersection_bundle(well=..)`).
    pub fn trajectory(&self) -> Vec<[f64; 3]> {
        if self.bore.trajectories().is_empty() {
            return Vec::new();
        }
        let (md0, md1) = self.bore.active().md_range();
        if md1 <= md0 {
            return Vec::new();
        }
        const N: usize = 48;
        (0..=N)
            .filter_map(|s| {
                let md = md0 + (md1 - md0) * (s as f64 / N as f64);
                self.bore.xyz(md).map(|p| [p.x, p.y, -p.z])
            })
            .collect()
    }

    /// A full-curve view of log `mnemonic` on this bore, or `None`.
    pub fn log(&self, mnemonic: &str) -> Option<LogView<'_>> {
        self.bore.log(mnemonic)
    }

    /// This bore's stratigraphic intervals (name + `[top_md, base_md)`), in depth
    /// order — used to colour a crossplot sample by the zone it falls in.
    pub fn zones(&self) -> Vec<Interval<'_>> {
        self.bore.zones()
    }
}

/// A loaded project: the petekio data substrate + the load inventory.
pub struct Project {
    geo: GeoData,
    inventory: Inventory,
    crs: Option<String>,
    /// Fluid-contact picks (`Type == "Other"`: GOC/FWL/OWC/…) parsed facade-side
    /// from the Petrel tops files — petekio's stratigraphic tops loader keeps
    /// only `Type == "Horizon"`, so these would otherwise be dropped (F3). Each
    /// is `(surface_name, well, depth_m)` with depth = subsea TVD (positive-down).
    contacts: Vec<(String, String, f64)>,
}

const SURFACE_EXTS: &[&str] = &["irap", "gri", "cps3grid"];
const POLYGON_EXTS: &[&str] = &["pol", "shp", "cps3lines"];
// `.irapclassicpoints` (scattered horizon picks) + `.earthvisiongrid` are
// scattered point-sets petekio loads natively; classified here as points so
// they load (framework grids a point-set horizon onto the build lattice — F2).
const POINT_EXTS: &[&str] = &["xyz", "dat", "irapclassicpoints", "earthvisiongrid"];

impl Project {
    /// Walk `root`, classify + load every recognised file, and return the
    /// populated project. `crs` is recorded as a provenance label (petekio does
    /// not reproject); `aliases` canonicalise log mnemonics at load.
    ///
    /// # Errors
    /// [`SrsError`] only on an unreadable `root`; per-file failures are recorded
    /// on the [`Inventory`] as skips, never aborting the walk.
    pub fn load(
        root: impl AsRef<Path>,
        crs: Option<String>,
        aliases: Option<Vec<(String, String)>>,
    ) -> Result<Self, SrsError> {
        let root = root.as_ref();
        if !root.is_dir() {
            return Err(SrsError::InvalidInput(format!(
                "project root '{}' is not a directory",
                root.display()
            )));
        }
        let mut geo = GeoData::new(Unit::Metres);
        if let Some(pairs) = aliases {
            geo.set_curve_aliases(NameMap::from_pairs(pairs));
        }
        let mut inv = Inventory::default();

        // 1) Wells first (so tops can attach to already-loaded wells). Discover
        //    each well id from the survey/log filenames under `wells/`, then route
        //    it through petekio's directory-walking `load_well` (F1 + F6):
        //
        //    - `load_well` walks the tree **recursively** and keeps only this
        //      well's files, so the real Petrel split `Wells/Paths/` +
        //      `Wells/Logs/` layout Just Works (no more one-bogus-well-per-subdir).
        //    - A single `.wellpath` is the well's one (main) bore, so its logs
        //      position on the **real deviated trajectory** — no synthesized
        //      vertical (F6) — and the `.wellpath` header supplies the
        //      authoritative wellhead datum (we pass a zero placeholder).
        let mut tops_files: Vec<PathBuf> = Vec::new();
        for wells_dir in child_dirs_named(root, "wells") {
            let ids = discover_well_ids(&wells_dir);
            for id in &ids {
                // Discovery is by well id; the bore-level inventory is derived from
                // the loaded substrate at the end (so `inventory().wells` matches
                // `proj.wells()` exactly). A load failure is a loud skip.
                if let Err(e) = geo.load_well(id, (0.0, 0.0), 0.0, &wells_dir) {
                    inv.skipped.push((
                        format!("{} [well '{id}']", wells_dir.display()),
                        e.to_string(),
                    ));
                }
            }
            // F8: LAS 3.0 **core** files (stem carries "core") are folded into a
            // bore by `load_well` above and would otherwise leave no trace. Record
            // each against its now-loaded well, naming the bore it merged into
            // (mirroring petekio's filename→bore routing).
            record_core_merges(&geo, &wells_dir, &ids, &mut inv.merged);
        }

        // 2) Everything else by extension / parent-dir hint (recursive), holding
        //    `.tops` files for step 3.
        let mut files = Vec::new();
        collect_files(root, &mut files);
        for path in &files {
            if under_dir(path, "wells") {
                continue; // already ingested well-side
            }
            let ext = ext_of(path);
            let parent = parent_name(path);
            if is_tops_file(path, &ext) {
                // `.tops`/`.petrel`, or an extension-less Petrel tops export whose
                // name carries "tops" (e.g. `FieldTops`) — the real raw layout (F3).
                tops_files.push(path.clone());
            } else if SURFACE_EXTS.contains(&ext.as_str()) {
                load_named(&mut geo, path, &mut inv, Kind::Surface);
            } else if ext == "geojson" || ext == "json" {
                // geojson is polygons or points — disambiguate by parent dir.
                if dir_is(&parent, &["point", "scatter"]) {
                    load_named(&mut geo, path, &mut inv, Kind::Points);
                } else {
                    load_named(&mut geo, path, &mut inv, Kind::Polygons);
                }
            } else if POLYGON_EXTS.contains(&ext.as_str()) {
                load_named(&mut geo, path, &mut inv, Kind::Polygons);
            } else if ext == "csv" || POINT_EXTS.contains(&ext.as_str()) {
                load_named(&mut geo, path, &mut inv, Kind::Points);
            } else if ext == "las" || ext == "wellpath" {
                // Loose well files outside a wells/ dir — not classifiable to a
                // well here; record rather than guess.
                inv.skipped.push((
                    path.display().to_string(),
                    "well file outside a wells/<id>/ directory".into(),
                ));
            } else {
                inv.skipped.push((
                    path.display().to_string(),
                    format!("unrecognised extension '.{ext}'"),
                ));
            }
        }

        // 3) Tops (multi-well Petrel picks) after the wells exist. petekio
        //    distributes the `Type == "Horizon"` (stratigraphic) picks; the
        //    facade separately parses the `Type == "Other"` fluid contacts
        //    (GOC/FWL/OWC), which petekio's tops loader drops (F3).
        let mut contacts: Vec<(String, String, f64)> = Vec::new();
        for tp in tops_files {
            if let Err(e) = geo.load_well_tops(&tp) {
                inv.skipped.push((tp.display().to_string(), e.to_string()));
            }
            contacts.extend(parse_other_contacts(&tp));
        }
        inv.tops = discover_tops(&geo, &contacts);
        inv.surfaces.sort();
        inv.surfaces.dedup();

        let mut this = Self {
            geo,
            inventory: inv,
            crs,
            contacts,
        };
        // Inventory wells are **bore-level**, derived from the same positioned
        // bores `proj.wells()` yields — one entry per positioned bore under its
        // bore-qualified id (`"99/9-1 A"`; a single-bore well keeps its plain id).
        // So the inventory advertises the exact ids the modelling seams consume,
        // not the parent well id.
        this.inventory.wells = this.wells().iter().map(|bw| bw.id.clone()).collect();
        Ok(this)
    }

    /// The load inventory (loaded artifacts + skipped-with-reason).
    pub fn inventory(&self) -> &Inventory {
        &self.inventory
    }

    /// The recorded CRS provenance label, if any.
    pub fn crs(&self) -> Option<&str> {
        self.crs.as_deref()
    }

    /// The loaded wells as **bore-level** entries — one [`BoreWell`] per
    /// positioned bore (the R-a wiring). A single-bore well yields one entry
    /// under its plain id; a multi-sidetrack well (`99/9-1` with bores A/B/ST2)
    /// yields one entry *per bore* under `"99/9-1 A"` etc., each positioned by
    /// that bore's own trajectory. Unpositioned bores (e.g. the empty main bore
    /// of a multi-bore well) carry no penetration and are skipped, so the
    /// framework tie + property upscaling see exactly the real, positioned bores.
    pub fn wells(&self) -> Vec<BoreWell<'_>> {
        let mut out = Vec::new();
        for well in self.geo.wells().iter() {
            for bore in well.sidetracks() {
                if bore.trajectories().is_empty() {
                    continue; // no trajectory → not an independent positioned bore
                }
                out.push(BoreWell {
                    id: well.bore_id(&bore.label),
                    bore,
                });
            }
        }
        out
    }

    /// A loaded surface by name (grid or gridded scattered set), or `None`.
    pub fn surface(&self, name: &str) -> Option<&Surface> {
        self.geo.surface(name)
    }

    /// The raw stored value of a loaded surface at world `(x, y)` — the nearest
    /// node on the surface's own axis-aligned lattice (`NaN` for a degenerate,
    /// zero-spacing geometry). The value is returned as stored: petekio's
    /// negative-down elevation for a depth surface, or the value itself for a
    /// value grid (a depositional-trend map). The same node→world mapping the
    /// collocated trend reads, so a value sampled here matches what steers a
    /// collocated population.
    ///
    /// # Errors
    /// [`SrsError::InvalidInput`] if `name` names no loaded surface.
    pub fn surface_value_at(&self, name: &str, x: f64, y: f64) -> Result<f64, SrsError> {
        let s = self.require_surface(name)?;
        Ok(nearest_node_value(s, x, y))
    }

    /// Batched [`surface_value_at`](Self::surface_value_at): sample the surface at
    /// each `(x, y)` in `points` with one name→surface resolution (the chatty
    /// per-point lookup collapsed to a single crossing + loop).
    ///
    /// # Errors
    /// [`SrsError::InvalidInput`] if `name` names no loaded surface.
    pub fn surface_values_at(&self, name: &str, points: &[[f64; 2]]) -> Result<Vec<f64>, SrsError> {
        let s = self.require_surface(name)?;
        Ok(points
            .iter()
            .map(|p| nearest_node_value(s, p[0], p[1]))
            .collect())
    }

    /// Build a collocated-cokriging [`TrendSpec`] from a loaded surface.
    ///
    /// **A trend is not a depth.** With `as_depth = false` the surface is read as
    /// a trend in its own units (a net-sand fraction, an amplitude, …), so `corr`
    /// reads against the surface value directly. With `as_depth = true` a
    /// structural surface's petekio negative-down elevation is flipped to
    /// positive-down depth (the framework's Z-flip) so `corr > 0` means the
    /// property increases with depth. Either way the field is shifted up to be
    /// non-negative when needed (the seam's trend is a multiplier domain);
    /// because the collocated secondary is standardized internally (normal
    /// scores), the additive offset is steering-neutral — it only satisfies the
    /// constructor.
    ///
    /// # Errors
    /// [`SrsError::InvalidInput`] if `name` names no loaded surface; the trend
    /// constructor's error otherwise.
    pub fn collocated_trend(
        &self,
        name: &str,
        corr: f64,
        as_depth: bool,
    ) -> Result<TrendSpec, SrsError> {
        let s = self.require_surface(name)?;
        let g = &s.geom;
        let (ncol, nrow) = (g.ncol, g.nrow);
        let vals = s.values();
        let mut values = vec![f64::NAN; ncol * nrow];
        let mut min_finite = f64::INFINITY;
        for j in 0..nrow {
            for i in 0..ncol {
                // as_depth: structural elevation → positive-down depth; else the
                // trend value as-is (its own units).
                let v = if as_depth {
                    -vals[[i, j]]
                } else {
                    vals[[i, j]]
                };
                values[j * ncol + i] = v;
                if v.is_finite() {
                    min_finite = min_finite.min(v);
                }
            }
        }
        // Shift any negative field up into the multiplier domain (steering-neutral).
        if min_finite.is_finite() && min_finite < 0.0 {
            for v in values.iter_mut() {
                if v.is_finite() {
                    *v -= min_finite;
                }
            }
        }
        TrendSpec::collocated(ncol, nrow, values, g.xori, g.yori, g.xinc, g.yinc, corr)
    }

    /// Resolve a loaded surface by name or raise the shared "not loaded" error.
    fn require_surface(&self, name: &str) -> Result<&Surface, SrsError> {
        self.geo
            .surface(name)
            .ok_or_else(|| SrsError::InvalidInput(format!("surface '{name}' not loaded")))
    }

    /// A named horizon as an owned gridded [`Surface`] (F2): a loaded **grid**
    /// surface is returned as-is; a loaded **scattered point-set** (real seismic
    /// horizons arrive as `.IrapClassicPoints`) is gridded onto `lattice` with
    /// Briggs minimum-curvature — the same kernel the Rust stack uses — via
    /// petekio's `PointSet::to_surface`. `lattice = None` derives a build lattice
    /// from the point-set's bounding box (the first horizon establishes the shared
    /// framework lattice; later horizons pass `Some(top_geom)` to grid onto it).
    ///
    /// `cell_size_m` sets the derived lattice's cell size when `lattice = None`
    /// (the first, lattice-establishing horizon): the node counts are `ceil(extent
    /// / cell_size) + 1` per axis, capped loudly at [`MAX_CELLS_PER_AXIS`]. `None`
    /// falls back to the legacy fixed [`DEFAULT_TARGET_CELLS`]²; ignored once a
    /// `lattice` is supplied.
    ///
    /// # Errors
    /// [`SrsError::InvalidInput`] if `name` names neither a loaded surface nor a
    /// loaded point-set (or the point-set is empty); the gridder's error otherwise.
    pub fn horizon_surface(
        &self,
        name: &str,
        lattice: Option<&GridGeometry>,
        cell_size_m: Option<f64>,
    ) -> Result<Surface, SrsError> {
        if let Some(s) = self.geo.surface(name) {
            return Ok(s.clone());
        }
        if let Some(ps) = self.geo.points(name) {
            if ps.is_empty() {
                return Err(SrsError::InvalidInput(format!(
                    "horizon point-set '{name}' has no points to grid"
                )));
            }
            let geom = lattice
                .cloned()
                .unwrap_or_else(|| lattice_from_bbox(&ps.bbox(), cell_size_m));
            return ps
                .to_surface(geom, GridMethod::MinimumCurvature)
                .map_err(SrsError::from);
        }
        Err(SrsError::InvalidInput(format!(
            "horizon '{name}' is neither a loaded surface nor a loaded point-set"
        )))
    }

    /// The horizon's build **lattice geometry only** — cheaply, WITHOUT gridding.
    ///
    /// A loaded grid surface returns its own geometry; a scattered point-set returns
    /// the lattice derived from its bounding box ([`lattice_from_bbox`], identical to
    /// the lattice [`Self::horizon_surface`] would grid it onto) — but skips the
    /// expensive minimum-curvature solve. Used by the framework to fix the shared
    /// areal lattice (`top_geom`) for the **raw-scatter stack** path, where the top's
    /// gridded VALUES are never used (the engine conditions the scatter itself) — only
    /// its lattice is needed. `cell_size_m` sizes a derived point-set lattice (ignored
    /// for a loaded surface); see [`Self::horizon_surface`].
    ///
    /// # Errors
    /// [`SrsError::InvalidInput`] if `name` names neither a loaded surface nor a loaded
    /// (non-empty) point-set.
    pub fn horizon_geom(
        &self,
        name: &str,
        cell_size_m: Option<f64>,
    ) -> Result<GridGeometry, SrsError> {
        if let Some(s) = self.geo.surface(name) {
            return Ok(s.geom.clone());
        }
        if let Some(ps) = self.geo.points(name) {
            if ps.is_empty() {
                return Err(SrsError::InvalidInput(format!(
                    "horizon point-set '{name}' has no points to grid"
                )));
            }
            return Ok(lattice_from_bbox(&ps.bbox(), cell_size_m));
        }
        Err(SrsError::InvalidInput(format!(
            "horizon '{name}' is neither a loaded surface nor a loaded point-set"
        )))
    }

    /// The underlying petekio substrate (for the framework/property seams).
    pub fn geo(&self) -> &GeoData {
        &self.geo
    }

    /// Per-well picks of a named tops surface (e.g. "GOC"/"FWL") + a
    /// representative level with spread. `None` if no well carries the pick.
    ///
    /// `wells` optionally restricts the aggregation to a subset of bores/wells: a
    /// pick is kept when its bore-qualified id is listed, or a listed **well** id
    /// prefixes it (`"99/9-1"` selects both `"99/9-1 A"` and `"99/9-1 B"`). `None`
    /// (or an empty list) aggregates across every well.
    pub fn pick(&self, name: &str, wells: Option<&[String]>) -> Option<TopsPick> {
        let keep = |id: &str| match wells {
            None | Some([]) => true,
            Some(ws) => ws
                .iter()
                .any(|w| id == w || id.starts_with(&format!("{w} "))),
        };
        let mut picks: Vec<(String, f64)> = Vec::new();
        // Stratigraphic (Type="Horizon") picks distributed onto the wells — read
        // **per bore** (a multi-sidetrack well's picks live on its bores; the
        // well-level accessor would resolve-empty on a multi-bore well). Depth is
        // subsea TVD (`-z`) at the pick MD on the bore's own trajectory.
        for bw in self.wells() {
            if !keep(&bw.id) {
                continue;
            }
            if let Some(iv) = bw.top(name) {
                if let Some(d) = bw.xyz(iv.top_md).map(|p| -p.z) {
                    picks.push((bw.id.clone(), d));
                }
            }
        }
        // Fluid-contact (Type="Other") picks parsed facade-side (F3): GOC/FWL/OWC
        // are exactly what the two-contact `grid.model(...)` API needs.
        for (surf, well, depth) in &self.contacts {
            if surf.eq_ignore_ascii_case(name) && keep(well) {
                picks.push((well.clone(), *depth));
            }
        }
        if picks.is_empty() {
            return None;
        }
        let n = picks.len() as f64;
        let level = picks.iter().map(|(_, d)| *d).sum::<f64>() / n;
        let spread = if picks.len() < 2 {
            0.0
        } else {
            let var = picks.iter().map(|(_, d)| (d - level).powi(2)).sum::<f64>() / (n - 1.0);
            var.sqrt()
        };
        Some(TopsPick {
            name: name.to_string(),
            picks,
            level_m: level,
            spread_m: spread,
        })
    }
}

#[cfg(test)]
mod tests {
    //! F3 regression: fluid-contact picks (Petrel `Type == "Other"`: GOC/FWL) are
    //! parsed facade-side, since petekio's tops loader keeps only `Type ==
    //! "Horizon"`. Synthesized to the Petrel well-tops shape (headered, quoted
    //! Surface/Well fields incl. a bore suffix with a space).
    use super::*;

    fn write_tmp(name: &str, body: &str) -> PathBuf {
        let p = std::env::temp_dir().join(format!("srscore_f3_{}_{name}", std::process::id()));
        std::fs::write(&p, body).unwrap();
        p
    }

    #[test]
    fn cell_size_derives_lattice_from_extent() {
        // A synthetic ~30 km × 20 km extent at 100 m cells → ~300 × 200 cells
        // (301 × 201 nodes), not the fixed 100² fallback. (coordinator addendum:
        // regional point-horizon frameworks need resolution control.)
        let b = BBox {
            xmin: 500_000.0,
            xmax: 530_000.0, // 30 km
            ymin: 6_700_000.0,
            ymax: 6_720_000.0, // 20 km
        };
        let g = lattice_from_bbox(&b, Some(100.0));
        assert_eq!(g.ncol, 301, "30 km / 100 m → 300 cells → 301 nodes");
        assert_eq!(g.nrow, 201, "20 km / 100 m → 200 cells → 201 nodes");
        assert!((g.xinc - 100.0).abs() < 1e-9 && (g.yinc - 100.0).abs() < 1e-9);
        // No cell_size_m → the legacy fixed 100²-node fallback, extent-independent.
        let f = lattice_from_bbox(&b, None);
        assert_eq!((f.ncol, f.nrow), (101, 101));
    }

    #[test]
    fn cell_size_caps_loudly_at_the_axis_maximum() {
        // A sub-metre cell over a basin extent would blow the budget → capped.
        let b = BBox {
            xmin: 0.0,
            xmax: 1_000_000.0,
            ymin: 0.0,
            ymax: 10.0,
        };
        let g = lattice_from_bbox(&b, Some(0.5));
        assert_eq!(g.ncol, MAX_CELLS_PER_AXIS + 1, "x-axis capped");
        assert_eq!(
            g.nrow, 21,
            "small y-axis uncapped (10 m / 0.5 m = 20 cells)"
        );
    }

    #[test]
    fn parse_other_contacts_surfaces_goc_fwl() {
        // GOC/FWL rows are Type="Other"; the strat tops are Type="Horizon".
        // Z is negative-down elevation → depth = -Z (positive-down subsea).
        let body = "\
# Petrel well tops
VERSION 2
BEGIN HEADER
X
Y
Z
MD
Type
Surface
Well
END HEADER
431600 6521600 -2000.0 2025.0 Horizon \"Top Reservoir\" \"99/9-1 A\"
431600 6521600 -2111.0 2248.0 Other \"GOC\" \"99/9-1 A\"
431600 6521600 -2190.0 2325.0 Other \"FWL\" \"99/9-1 A\"
";
        let p = write_tmp("tops.tops", body);
        let contacts = parse_other_contacts(&p);
        std::fs::remove_file(&p).ok();

        // Only the two Type="Other" rows; the Horizon row is left to petekio.
        assert_eq!(contacts.len(), 2, "{contacts:?}");
        let goc = contacts.iter().find(|(s, _, _)| s == "GOC").unwrap();
        assert!((goc.2 - 2111.0).abs() < 1e-9, "GOC depth {}", goc.2);
        assert_eq!(goc.1, "99/9-1 A"); // quoted well with a space + bore suffix
        let fwl = contacts.iter().find(|(s, _, _)| s == "FWL").unwrap();
        assert!((fwl.2 - 2190.0).abs() < 1e-9);
    }

    #[test]
    fn split_petrel_tokens_keeps_quoted_fields() {
        let toks = split_petrel_tokens("1 2 3 Other \"Field Main top\" \"99/9-1 ST2\"");
        assert_eq!(
            toks,
            ["1", "2", "3", "Other", "Field Main top", "99/9-1 ST2"]
        );
    }

    #[test]
    fn parse_other_contacts_decodes_latin1() {
        // R-b: a real Petrel tops export is often Latin-1, not UTF-8. A high byte
        // in a name (here `å` = 0xE5, `æ` = 0xE6 — "Blåbær") must decode rather
        // than abort the whole parse (the old `read_to_string` returned empty).
        let mut body: Vec<u8> =
            b"BEGIN HEADER\nX\nY\nZ\nMD\nType\nSurface\nWell\nEND HEADER\n".to_vec();
        // 431600 6521600 -2111.0 2248.0 Other "Blåbær" "15/9-A1"
        body.extend_from_slice(b"431600 6521600 -2111.0 2248.0 Other \"");
        body.extend_from_slice(&[b'B', b'l', 0xE5, b'b', 0xE6, b'r']); // Blåbær (Latin-1)
        body.extend_from_slice(b"\" \"15/9-A1\"\n");
        let p = std::env::temp_dir().join(format!("srscore_rb_{}.tops", std::process::id()));
        std::fs::write(&p, &body).unwrap();
        let contacts = parse_other_contacts(&p);
        std::fs::remove_file(&p).ok();

        assert_eq!(contacts.len(), 1, "{contacts:?}");
        assert_eq!(contacts[0].0, "Blåbær"); // decoded, not dropped
        assert!((contacts[0].2 - 2111.0).abs() < 1e-9);
    }
}