oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! Deformation model (`defmodel`) — port of PROJ 9.8.0
//! `src/transformations/defmodel.cpp` and the header-only engine
//! `defmodel.hpp` / `defmodel_impl.hpp`.
//!
//! `+proj=defmodel +model=<model.json>` applies a time-dependent deformation
//! described by a JSON *master file*. The file lists one or more *components*,
//! each pairing a *spatial model* (a grid, referenced by filename) with a *time
//! function*. For an input coordinate `(lon, lat, z, t)`, every component
//! evaluates its time function at epoch `t` to obtain a scale factor, samples
//! its grid at `(lon, lat)` for the local displacement, scales it by the factor
//! and sums the contributions across components.
//!
//! Time functions ported exactly per PROJ: `constant`, `velocity`, `step`,
//! `reverse_step`, `piecewise`, `exponential`. (PROJ 9.8.0 defines no
//! `logarithmic` function.)
//!
//! Pure Rust, zero additional dependencies: the JSON master file is parsed by a
//! small hand-rolled recursive-descent parser (mirroring the tokenizer approach
//! of `oxiproj-grids`'s `tinshift`). Component grids are obtained through the
//! [`GridRegistry`](crate::GridRegistry) dependency inversion already used by
//! `gridshift`/`deformation` — they are supplied by the engine, not read here.
//!
//! Supported, matching PROJ's common case: a geographic definition CRS (radians
//! in/out); horizontal offsets in `degree` or `metre`; `horizontal_offset_method`
//! `addition` or `geocentric`; interpolation `bilinear` or
//! `geocentric_bilinear`; displacement types `none`, `horizontal`, `vertical`,
//! `3d`. CRS resolution is not available within this crate, so the definition
//! CRS is treated as geographic (all real-world deformation models are).

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjError, ProjResult, DEG_TO_RAD, M_TWOPI};
use oxiproj_grids::{read_geotiff, read_gtx, read_ntv2, GridSet};

// =====================================================================
// Hand-rolled recursive-descent JSON parser (object/array/string/number/
// bool/null). Zero dependencies.
// =====================================================================

#[derive(Debug, Clone)]
enum Json {
    Null,
    Bool,
    Num(f64),
    Str(String),
    Arr(Vec<Json>),
    Obj(Vec<(String, Json)>),
}

impl Json {
    fn get(&self, key: &str) -> Option<&Json> {
        match self {
            Json::Obj(entries) => entries.iter().find(|(k, _)| k == key).map(|(_, v)| v),
            _ => None,
        }
    }
    fn as_str(&self) -> Option<&str> {
        match self {
            Json::Str(s) => Some(s.as_str()),
            _ => None,
        }
    }
    fn as_f64(&self) -> Option<f64> {
        match self {
            Json::Num(n) => Some(*n),
            _ => None,
        }
    }
    fn as_arr(&self) -> Option<&[Json]> {
        match self {
            Json::Arr(a) => Some(a.as_slice()),
            _ => None,
        }
    }
    fn is_object(&self) -> bool {
        matches!(self, Json::Obj(_))
    }
}

/// Maximum nesting depth (objects/arrays) the master-file JSON parser will
/// descend into before giving up. Guards against a crafted `+model=` file
/// with deeply nested arrays/objects driving unbounded native-stack
/// recursion (stack-overflow DoS), mirroring the depth cap used by the WKT
/// parser in `oxiproj-crs`.
const MAX_JSON_DEPTH: usize = 64;

struct JsonParser<'a> {
    bytes: &'a [u8],
    pos: usize,
    depth: usize,
}

impl<'a> JsonParser<'a> {
    fn parse(s: &str) -> ProjResult<Json> {
        let mut p = JsonParser {
            bytes: s.as_bytes(),
            pos: 0,
            depth: 0,
        };
        p.skip_ws();
        let v = p.parse_value()?;
        p.skip_ws();
        if p.pos != p.bytes.len() {
            return Err(ProjError::IllegalArgValue);
        }
        Ok(v)
    }

    fn skip_ws(&mut self) {
        while let Some(&b) = self.bytes.get(self.pos) {
            match b {
                b' ' | b'\t' | b'\n' | b'\r' => self.pos += 1,
                _ => break,
            }
        }
    }

    fn peek(&self) -> Option<u8> {
        self.bytes.get(self.pos).copied()
    }

    fn expect(&mut self, b: u8) -> ProjResult<()> {
        if self.peek() == Some(b) {
            self.pos += 1;
            Ok(())
        } else {
            Err(ProjError::IllegalArgValue)
        }
    }

    fn parse_value(&mut self) -> ProjResult<Json> {
        self.skip_ws();
        match self.peek().ok_or(ProjError::IllegalArgValue)? {
            b'{' => self.parse_object(),
            b'[' => self.parse_array(),
            b'"' => Ok(Json::Str(self.parse_string()?)),
            b't' | b'f' => self.parse_bool(),
            b'n' => self.parse_null(),
            b'-' | b'0'..=b'9' => self.parse_number(),
            _ => Err(ProjError::IllegalArgValue),
        }
    }

    fn parse_object(&mut self) -> ProjResult<Json> {
        self.expect(b'{')?;
        self.enter_nesting()?;
        let result = self.parse_object_body();
        self.depth -= 1;
        result
    }

    fn parse_object_body(&mut self) -> ProjResult<Json> {
        let mut entries = Vec::new();
        self.skip_ws();
        if self.peek() == Some(b'}') {
            self.pos += 1;
            return Ok(Json::Obj(entries));
        }
        loop {
            self.skip_ws();
            let key = self.parse_string()?;
            self.skip_ws();
            self.expect(b':')?;
            let val = self.parse_value()?;
            entries.push((key, val));
            self.skip_ws();
            match self.peek().ok_or(ProjError::IllegalArgValue)? {
                b',' => self.pos += 1,
                b'}' => {
                    self.pos += 1;
                    break;
                }
                _ => return Err(ProjError::IllegalArgValue),
            }
        }
        Ok(Json::Obj(entries))
    }

    fn parse_array(&mut self) -> ProjResult<Json> {
        self.expect(b'[')?;
        self.enter_nesting()?;
        let result = self.parse_array_body();
        self.depth -= 1;
        result
    }

    fn parse_array_body(&mut self) -> ProjResult<Json> {
        let mut items = Vec::new();
        self.skip_ws();
        if self.peek() == Some(b']') {
            self.pos += 1;
            return Ok(Json::Arr(items));
        }
        loop {
            let val = self.parse_value()?;
            items.push(val);
            self.skip_ws();
            match self.peek().ok_or(ProjError::IllegalArgValue)? {
                b',' => self.pos += 1,
                b']' => {
                    self.pos += 1;
                    break;
                }
                _ => return Err(ProjError::IllegalArgValue),
            }
        }
        Ok(Json::Arr(items))
    }

    /// Enter one level of object/array nesting, rejecting master files that
    /// exceed [`MAX_JSON_DEPTH`] before any further recursive descent takes
    /// place (defends against stack-overflow DoS from crafted input).
    fn enter_nesting(&mut self) -> ProjResult<()> {
        if self.depth >= MAX_JSON_DEPTH {
            return Err(ProjError::IllegalArgValue);
        }
        self.depth += 1;
        Ok(())
    }

    fn parse_string(&mut self) -> ProjResult<String> {
        self.expect(b'"')?;
        let mut out = String::new();
        loop {
            let c = self.peek().ok_or(ProjError::IllegalArgValue)?;
            self.pos += 1;
            match c {
                b'"' => break,
                b'\\' => {
                    let e = self.peek().ok_or(ProjError::IllegalArgValue)?;
                    self.pos += 1;
                    match e {
                        b'"' => out.push('"'),
                        b'\\' => out.push('\\'),
                        b'/' => out.push('/'),
                        b'b' => out.push('\u{0008}'),
                        b'f' => out.push('\u{000C}'),
                        b'n' => out.push('\n'),
                        b'r' => out.push('\r'),
                        b't' => out.push('\t'),
                        b'u' => {
                            let cp = self.parse_hex4()?;
                            if (0xD800..=0xDBFF).contains(&cp) {
                                self.expect(b'\\')?;
                                self.expect(b'u')?;
                                let lo = self.parse_hex4()?;
                                if !(0xDC00..=0xDFFF).contains(&lo) {
                                    return Err(ProjError::IllegalArgValue);
                                }
                                let combined = 0x10000u32
                                    + (((cp - 0xD800) as u32) << 10)
                                    + (lo - 0xDC00) as u32;
                                out.push(
                                    char::from_u32(combined).ok_or(ProjError::IllegalArgValue)?,
                                );
                            } else {
                                out.push(
                                    char::from_u32(cp as u32).ok_or(ProjError::IllegalArgValue)?,
                                );
                            }
                        }
                        _ => return Err(ProjError::IllegalArgValue),
                    }
                }
                _ if c < 0x80 => out.push(c as char),
                _ => {
                    // Multi-byte UTF-8: copy the whole sequence.
                    let start = self.pos - 1;
                    let n = utf8_len(c);
                    let end = start + n;
                    if end > self.bytes.len() {
                        return Err(ProjError::IllegalArgValue);
                    }
                    let s = core::str::from_utf8(&self.bytes[start..end])
                        .map_err(|_| ProjError::IllegalArgValue)?;
                    out.push_str(s);
                    self.pos = end;
                }
            }
        }
        Ok(out)
    }

    fn parse_hex4(&mut self) -> ProjResult<u16> {
        let mut v: u16 = 0;
        for _ in 0..4 {
            let c = self.peek().ok_or(ProjError::IllegalArgValue)?;
            self.pos += 1;
            let d = match c {
                b'0'..=b'9' => (c - b'0') as u16,
                b'a'..=b'f' => (c - b'a' + 10) as u16,
                b'A'..=b'F' => (c - b'A' + 10) as u16,
                _ => return Err(ProjError::IllegalArgValue),
            };
            v = v * 16 + d;
        }
        Ok(v)
    }

    fn parse_bool(&mut self) -> ProjResult<Json> {
        if self.bytes[self.pos..].starts_with(b"true") {
            self.pos += 4;
            Ok(Json::Bool)
        } else if self.bytes[self.pos..].starts_with(b"false") {
            self.pos += 5;
            Ok(Json::Bool)
        } else {
            Err(ProjError::IllegalArgValue)
        }
    }

    fn parse_null(&mut self) -> ProjResult<Json> {
        if self.bytes[self.pos..].starts_with(b"null") {
            self.pos += 4;
            Ok(Json::Null)
        } else {
            Err(ProjError::IllegalArgValue)
        }
    }

    fn parse_number(&mut self) -> ProjResult<Json> {
        let start = self.pos;
        if self.peek() == Some(b'-') {
            self.pos += 1;
        }
        while matches!(self.peek(), Some(b'0'..=b'9')) {
            self.pos += 1;
        }
        if self.peek() == Some(b'.') {
            self.pos += 1;
            while matches!(self.peek(), Some(b'0'..=b'9')) {
                self.pos += 1;
            }
        }
        if matches!(self.peek(), Some(b'e') | Some(b'E')) {
            self.pos += 1;
            if matches!(self.peek(), Some(b'+') | Some(b'-')) {
                self.pos += 1;
            }
            while matches!(self.peek(), Some(b'0'..=b'9')) {
                self.pos += 1;
            }
        }
        let s = core::str::from_utf8(&self.bytes[start..self.pos])
            .map_err(|_| ProjError::IllegalArgValue)?;
        let v: f64 = s.parse().map_err(|_| ProjError::IllegalArgValue)?;
        Ok(Json::Num(v))
    }
}

fn utf8_len(b: u8) -> usize {
    if b >> 5 == 0b110 {
        2
    } else if b >> 4 == 0b1110 {
        3
    } else if b >> 3 == 0b11110 {
        4
    } else {
        1
    }
}

// JSON field accessors ------------------------------------------------

fn req_str(j: &Json, key: &str) -> ProjResult<String> {
    j.get(key)
        .and_then(|v| v.as_str())
        .map(|s| s.to_owned())
        .ok_or(ProjError::IllegalArgValue)
}
fn opt_str(j: &Json, key: &str) -> String {
    j.get(key).and_then(|v| v.as_str()).unwrap_or("").to_owned()
}
fn req_f64(j: &Json, key: &str) -> ProjResult<f64> {
    j.get(key)
        .and_then(|v| v.as_f64())
        .ok_or(ProjError::IllegalArgValue)
}
fn req_obj<'a>(j: &'a Json, key: &str) -> ProjResult<&'a Json> {
    j.get(key)
        .filter(|v| v.is_object())
        .ok_or(ProjError::IllegalArgValue)
}
fn req_arr<'a>(j: &'a Json, key: &str) -> ProjResult<&'a [Json]> {
    j.get(key)
        .and_then(|v| v.as_arr())
        .ok_or(ProjError::IllegalArgValue)
}

// =====================================================================
// ISO 8601 date-time -> decimal year (ported from defmodel_impl.hpp).
// =====================================================================

fn iso8601_to_decimal_year(dt: &str) -> ProjResult<f64> {
    let (date, time) = dt.split_once('T').ok_or(ProjError::IllegalArgValue)?;
    let time = time.strip_suffix('Z').unwrap_or(time);
    let mut di = date.split('-');
    let year: i64 = di
        .next()
        .and_then(|s| s.parse().ok())
        .ok_or(ProjError::IllegalArgValue)?;
    let month: i64 = di
        .next()
        .and_then(|s| s.parse().ok())
        .ok_or(ProjError::IllegalArgValue)?;
    let day: i64 = di
        .next()
        .and_then(|s| s.parse().ok())
        .ok_or(ProjError::IllegalArgValue)?;
    let mut ti = time.split(':');
    let hour: i64 = ti
        .next()
        .and_then(|s| s.parse().ok())
        .ok_or(ProjError::IllegalArgValue)?;
    let min: i64 = ti
        .next()
        .and_then(|s| s.parse().ok())
        .ok_or(ProjError::IllegalArgValue)?;
    let sec: i64 = ti
        .next()
        .and_then(|s| s.parse().ok())
        .ok_or(ProjError::IllegalArgValue)?;
    if year < 1582
        || !(1..=12).contains(&month)
        || !(1..=31).contains(&day)
        || !(0..24).contains(&hour)
        || !(0..60).contains(&min)
        || !(0..61).contains(&sec)
    {
        return Err(ProjError::IllegalArgValue);
    }
    let is_leap = (year % 4 == 0 && year % 100 != 0) || year % 400 == 0;
    let months: [[i64; 12]; 2] = [
        [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
        [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
    ];
    let li = if is_leap { 1usize } else { 0usize };
    let mut day_in_year = day - 1;
    for m in 1..month {
        day_in_year += months[li][(m - 1) as usize];
    }
    if day > months[li][(month - 1) as usize] {
        return Err(ProjError::IllegalArgValue);
    }
    let seconds = (day_in_year * 86400 + hour * 3600 + min * 60 + sec) as f64;
    let year_seconds = if is_leap {
        86400.0 * 366.0
    } else {
        86400.0 * 365.0
    };
    Ok(year as f64 + seconds / year_seconds)
}

// =====================================================================
// Time functions (ported exactly from defmodel_impl.hpp).
// =====================================================================

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum Clamp {
    Zero,
    Constant,
    Linear,
}

impl Clamp {
    fn parse(s: &str) -> ProjResult<Clamp> {
        match s {
            "zero" => Ok(Clamp::Zero),
            "constant" => Ok(Clamp::Constant),
            "linear" => Ok(Clamp::Linear),
            _ => Err(ProjError::IllegalArgValue),
        }
    }
}

#[derive(Debug, Clone)]
enum TimeFunction {
    Constant,
    Velocity {
        reference_epoch: f64,
    },
    Step {
        step_epoch: f64,
    },
    ReverseStep {
        step_epoch: f64,
    },
    Piecewise {
        before_first: Clamp,
        after_last: Clamp,
        model: Vec<(f64, f64)>,
    },
    Exponential {
        reference_epoch: f64,
        end_epoch: Option<f64>,
        relaxation_constant: f64,
        before_scale_factor: f64,
        initial_scale_factor: f64,
        final_scale_factor: f64,
    },
}

impl TimeFunction {
    fn evaluate_at(&self, dt: f64) -> f64 {
        match self {
            TimeFunction::Constant => 1.0,
            TimeFunction::Velocity { reference_epoch } => dt - reference_epoch,
            TimeFunction::Step { step_epoch } => {
                if dt < *step_epoch {
                    0.0
                } else {
                    1.0
                }
            }
            TimeFunction::ReverseStep { step_epoch } => {
                if dt < *step_epoch {
                    -1.0
                } else {
                    0.0
                }
            }
            TimeFunction::Piecewise {
                before_first,
                after_last,
                model,
            } => eval_piecewise(*before_first, *after_last, model, dt),
            TimeFunction::Exponential {
                reference_epoch,
                end_epoch,
                relaxation_constant,
                before_scale_factor,
                initial_scale_factor,
                final_scale_factor,
            } => {
                let t0 = *reference_epoch;
                if dt < t0 {
                    return *before_scale_factor;
                }
                let mut dt = dt;
                if let Some(end) = end_epoch {
                    dt = dt.min(*end);
                }
                initial_scale_factor
                    + (final_scale_factor - initial_scale_factor)
                        * (1.0 - (-(dt - t0) / relaxation_constant).exp())
            }
        }
    }
}

fn eval_piecewise(before_first: Clamp, after_last: Clamp, model: &[(f64, f64)], dt: f64) -> f64 {
    if model.is_empty() {
        return 0.0;
    }
    let (dt1, f1) = model[0];
    if dt < dt1 {
        match before_first {
            Clamp::Zero => return 0.0,
            Clamp::Constant => return f1,
            Clamp::Linear => {
                if model.len() == 1 {
                    return f1;
                }
                let (dt2, f2) = model[1];
                if dt1 == dt2 {
                    return f1;
                }
                return (f1 * (dt2 - dt) + f2 * (dt - dt1)) / (dt2 - dt1);
            }
        }
    }
    for w in model.windows(2) {
        let (dti, fi) = w[0];
        let (dtip1, fip1) = w[1];
        if dt < dtip1 {
            if dti == dtip1 {
                continue;
            }
            return (fi * (dtip1 - dt) + fip1 * (dt - dti)) / (dtip1 - dti);
        }
    }
    match after_last {
        Clamp::Zero => 0.0,
        Clamp::Constant => model[model.len() - 1].1,
        Clamp::Linear => {
            if model.len() == 1 {
                return model[model.len() - 1].1;
            }
            let (dt_prev, f_prev) = model[model.len() - 2];
            let (dt_last, f_last) = model[model.len() - 1];
            if dt_prev == dt_last {
                return f_last;
            }
            (f_prev * (dt_last - dt) + f_last * (dt - dt_prev)) / (dt_last - dt_prev)
        }
    }
}

// =====================================================================
// Model schema.
// =====================================================================

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum DisplacementType {
    None,
    Horizontal,
    Vertical,
    ThreeD,
}

impl DisplacementType {
    fn parse(s: &str) -> ProjResult<DisplacementType> {
        match s {
            "none" => Ok(DisplacementType::None),
            "horizontal" => Ok(DisplacementType::Horizontal),
            "vertical" => Ok(DisplacementType::Vertical),
            "3d" => Ok(DisplacementType::ThreeD),
            _ => Err(ProjError::IllegalArgValue),
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum InterpMethod {
    Bilinear,
    GeocentricBilinear,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum HUnit {
    Degree,
    Metre,
    Unspecified,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum HMethod {
    Addition,
    Geocentric,
    Unspecified,
}

/// Bounding box, stored in radians (geographic definition CRS).
#[derive(Debug, Clone, Copy)]
struct Extent {
    minx: f64,
    miny: f64,
    maxx: f64,
    maxy: f64,
}

impl Extent {
    fn parse(j: &Json) -> ProjResult<Extent> {
        if req_str(j, "type")? != "bbox" {
            return Err(ProjError::IllegalArgValue);
        }
        let params = req_obj(j, "parameters")?;
        let bbox = req_arr(params, "bbox")?;
        if bbox.len() != 4 {
            return Err(ProjError::IllegalArgValue);
        }
        let mut v = [0.0f64; 4];
        for (i, e) in bbox.iter().enumerate() {
            v[i] = e.as_f64().ok_or(ProjError::IllegalArgValue)?;
        }
        Ok(Extent {
            minx: v[0] * DEG_TO_RAD,
            miny: v[1] * DEG_TO_RAD,
            maxx: v[2] * DEG_TO_RAD,
            maxy: v[3] * DEG_TO_RAD,
        })
    }
}

/// Component specification parsed from JSON (grid not yet loaded).
struct CompSpec {
    extent: Extent,
    displacement_type: DisplacementType,
    interp: InterpMethod,
    filename: String,
    time_function: TimeFunction,
}

fn parse_time_function(j: &Json) -> ProjResult<TimeFunction> {
    let typ = req_str(j, "type")?;
    if typ == "constant" {
        return Ok(TimeFunction::Constant);
    }
    let p = req_obj(j, "parameters")?;
    match typ.as_str() {
        "velocity" => Ok(TimeFunction::Velocity {
            reference_epoch: iso8601_to_decimal_year(&req_str(p, "reference_epoch")?)?,
        }),
        "step" => Ok(TimeFunction::Step {
            step_epoch: iso8601_to_decimal_year(&req_str(p, "step_epoch")?)?,
        }),
        "reverse_step" => Ok(TimeFunction::ReverseStep {
            step_epoch: iso8601_to_decimal_year(&req_str(p, "step_epoch")?)?,
        }),
        "piecewise" => {
            let before_first = Clamp::parse(&req_str(p, "before_first")?)?;
            let after_last = Clamp::parse(&req_str(p, "after_last")?)?;
            let model_arr = req_arr(p, "model")?;
            let mut model = Vec::with_capacity(model_arr.len());
            for elt in model_arr {
                if !elt.is_object() {
                    return Err(ProjError::IllegalArgValue);
                }
                let epoch = iso8601_to_decimal_year(&req_str(elt, "epoch")?)?;
                let scale = req_f64(elt, "scale_factor")?;
                model.push((epoch, scale));
            }
            Ok(TimeFunction::Piecewise {
                before_first,
                after_last,
                model,
            })
        }
        "exponential" => {
            let reference_epoch = iso8601_to_decimal_year(&req_str(p, "reference_epoch")?)?;
            let end_s = opt_str(p, "end_epoch");
            let end_epoch = if end_s.is_empty() {
                None
            } else {
                Some(iso8601_to_decimal_year(&end_s)?)
            };
            let relaxation_constant = req_f64(p, "relaxation_constant")?;
            if relaxation_constant <= 0.0 {
                return Err(ProjError::IllegalArgValue);
            }
            Ok(TimeFunction::Exponential {
                reference_epoch,
                end_epoch,
                relaxation_constant,
                before_scale_factor: req_f64(p, "before_scale_factor")?,
                initial_scale_factor: req_f64(p, "initial_scale_factor")?,
                final_scale_factor: req_f64(p, "final_scale_factor")?,
            })
        }
        _ => Err(ProjError::IllegalArgValue),
    }
}

fn parse_component(j: &Json) -> ProjResult<CompSpec> {
    if !j.is_object() {
        return Err(ProjError::IllegalArgValue);
    }
    let extent = Extent::parse(req_obj(j, "extent")?)?;
    let displacement_type = DisplacementType::parse(&req_str(j, "displacement_type")?)?;
    req_str(j, "uncertainty_type")?; // required by schema; value unused
    let spatial = req_obj(j, "spatial_model")?;
    req_str(spatial, "type")?; // required by schema; value unused
    let interp = match req_str(spatial, "interpolation_method")?.as_str() {
        "bilinear" => InterpMethod::Bilinear,
        "geocentric_bilinear" => InterpMethod::GeocentricBilinear,
        _ => return Err(ProjError::IllegalArgValue),
    };
    let filename = req_str(spatial, "filename")?;
    let time_function = parse_time_function(req_obj(j, "time_function")?)?;
    Ok(CompSpec {
        extent,
        displacement_type,
        interp,
        filename,
        time_function,
    })
}

/// Parsed master file (grids not yet loaded).
struct ModelSpec {
    extent: Extent,
    time_first: f64,
    time_last: f64,
    h_unit: HUnit,
    h_method: HMethod,
    components: Vec<CompSpec>,
}

fn parse_model(text: &str) -> ProjResult<ModelSpec> {
    let j = JsonParser::parse(text)?;
    if !j.is_object() {
        return Err(ProjError::IllegalArgValue);
    }
    req_str(&j, "file_type")?;
    req_str(&j, "format_version")?;
    let source_crs = req_str(&j, "source_crs")?;
    req_str(&j, "target_crs")?;
    let definition_crs = req_str(&j, "definition_crs")?;
    if source_crs != definition_crs {
        // PROJ: "source_crs != definition_crs not currently supported"
        return Err(ProjError::IllegalArgValue);
    }
    let h_unit = match j.get("horizontal_offset_unit").and_then(|v| v.as_str()) {
        Some("metre") => HUnit::Metre,
        Some("degree") => HUnit::Degree,
        Some(_) => return Err(ProjError::IllegalArgValue),
        None => HUnit::Unspecified,
    };
    let v_unit = j.get("vertical_offset_unit").and_then(|v| v.as_str());
    if let Some(u) = v_unit {
        if u != "metre" {
            return Err(ProjError::IllegalArgValue);
        }
    }
    let v_unit_metre = matches!(v_unit, Some("metre"));
    let h_method = match j.get("horizontal_offset_method").and_then(|v| v.as_str()) {
        Some("addition") => HMethod::Addition,
        Some("geocentric") => HMethod::Geocentric,
        Some(_) => return Err(ProjError::IllegalArgValue),
        None => HMethod::Unspecified,
    };
    let extent = Extent::parse(req_obj(&j, "extent")?)?;
    let te = req_obj(&j, "time_extent")?;
    let time_first = iso8601_to_decimal_year(&req_str(te, "first")?)?;
    let time_last = iso8601_to_decimal_year(&req_str(te, "last")?)?;

    let comps_arr = req_arr(&j, "components")?;
    let mut components = Vec::with_capacity(comps_arr.len());
    for c in comps_arr {
        let spec = parse_component(c)?;
        if matches!(
            spec.displacement_type,
            DisplacementType::Horizontal | DisplacementType::ThreeD
        ) {
            if h_unit == HUnit::Unspecified {
                return Err(ProjError::IllegalArgValue);
            }
            if h_method == HMethod::Unspecified {
                return Err(ProjError::IllegalArgValue);
            }
        }
        if matches!(
            spec.displacement_type,
            DisplacementType::Vertical | DisplacementType::ThreeD
        ) && !v_unit_metre
        {
            return Err(ProjError::IllegalArgValue);
        }
        if h_unit == HUnit::Degree && spec.interp != InterpMethod::Bilinear {
            return Err(ProjError::IllegalArgValue);
        }
        components.push(spec);
    }
    if h_unit == HUnit::Degree && h_method != HMethod::Addition {
        return Err(ProjError::IllegalArgValue);
    }

    Ok(ModelSpec {
        extent,
        time_first,
        time_last,
        h_unit,
        h_method,
        components,
    })
}

// =====================================================================
// Evaluator.
// =====================================================================

#[derive(Debug)]
struct Comp {
    extent: Extent,
    displacement_type: DisplacementType,
    interp: InterpMethod,
    time_function: TimeFunction,
    grid: Option<GridSet>,
}

#[derive(Debug)]
struct Defmodel {
    a: f64,
    b: f64,
    es: f64,
    is_unit_degree: bool,
    is_addition: bool,
    extent: Extent,
    time_first: f64,
    time_last: f64,
    components: Vec<Comp>,
}

const EPS_GEO: f64 = 1e-10;

/// Port of PROJ `bboxCheck`. Returns whether `(x,y)` is acceptable; for the
/// inverse computation it clamps slightly-outside points to the edges.
fn bbox_check(
    x: &mut f64,
    y: &mut f64,
    for_inverse: bool,
    b: Extent,
    eps: f64,
    extra: f64,
) -> bool {
    if *x < b.minx - eps || *x > b.maxx + eps || *y < b.miny - eps || *y > b.maxy + eps {
        if !for_inverse {
            return false;
        }
        let mut x_ok = false;
        if *x >= b.minx - eps && *x <= b.maxx + eps {
            x_ok = true;
        } else if *x > b.minx - extra && *x < b.minx {
            *x = b.minx;
            x_ok = true;
        } else if *x < b.maxx + extra && *x > b.maxx {
            *x = b.maxx;
            x_ok = true;
        }
        let mut y_ok = false;
        if *y >= b.miny - eps && *y <= b.maxy + eps {
            y_ok = true;
        } else if *y > b.miny - extra && *y < b.miny {
            *y = b.miny;
            y_ok = true;
        } else if *y < b.maxy + extra && *y > b.maxy {
            *y = b.maxy;
            y_ok = true;
        }
        x_ok && y_ok
    } else {
        true
    }
}

/// Convert an easting/northing offset (metres) to a longitude/latitude offset
/// (radians). Port of PROJ `DeltaEastingNorthingToLongLat`.
fn delta_en_to_lonlat(cosphi: f64, de: f64, dn: f64, a: f64, b: f64, es: f64) -> (f64, f64) {
    let one_minus_x = es * (1.0 - cosphi * cosphi);
    let xx = 1.0 - one_minus_x;
    let sqrtx = xx.sqrt();
    let dlam = de * sqrtx / (a * cosphi);
    let dphi = dn * a * sqrtx * xx / (b * b);
    (dlam, dphi)
}

fn geographic_to_geocentric(lam: f64, phi: f64, h: f64, a: f64, es: f64) -> (f64, f64, f64) {
    let sphi = phi.sin();
    let cphi = phi.cos();
    let n = a / (1.0 - es * sphi * sphi).sqrt();
    let x = (n + h) * cphi * lam.cos();
    let y = (n + h) * cphi * lam.sin();
    let z = (n * (1.0 - es) + h) * sphi;
    (x, y, z)
}

fn geocentric_to_geographic(x: f64, y: f64, z: f64, a: f64, es: f64) -> (f64, f64, f64) {
    let lam = y.atan2(x);
    let p = (x * x + y * y).sqrt();
    let mut phi = z.atan2(p * (1.0 - es));
    let mut h = 0.0;
    for _ in 0..10 {
        let sphi = phi.sin();
        let n = a / (1.0 - es * sphi * sphi).sqrt();
        h = p / phi.cos() - n;
        let new_phi = z.atan2(p * (1.0 - es * n / (n + h)));
        let converged = (new_phi - phi).abs() < 1e-13;
        phi = new_phi;
        if converged {
            break;
        }
    }
    (lam, phi, h)
}

impl Defmodel {
    /// Port of PROJ `Evaluator::forward`. Coordinates are radians (geographic).
    fn forward(
        &self,
        x_in: f64,
        y_in: f64,
        z_in: f64,
        t: f64,
        for_inverse: bool,
    ) -> Option<(f64, f64, f64)> {
        let mut x = x_in;
        let mut y = y_in;
        let mut x_out = x_in;
        let mut y_out = y_in;
        let mut z_out = z_in;
        let eps = EPS_GEO;

        // Global spatial extent (with longitude wrapping).
        while x < self.extent.minx - eps {
            x += M_TWOPI;
        }
        while x > self.extent.maxx + eps {
            x -= M_TWOPI;
        }
        let extra = 0.1 * DEG_TO_RAD;
        if !bbox_check(&mut x, &mut y, for_inverse, self.extent, eps, extra) {
            return None;
        }

        // Global temporal extent.
        if t < self.time_first || t > self.time_last {
            return None;
        }

        let mut dlam = 0.0;
        let mut dphi = 0.0;
        let mut de = 0.0;
        let mut dn = 0.0;
        let mut dz = 0.0;

        for comp in &self.components {
            if comp.displacement_type == DisplacementType::None {
                continue;
            }
            let grid = match &comp.grid {
                Some(g) => g,
                None => continue,
            };
            let mut xfg = x;
            let mut yfg = y;
            if !bbox_check(&mut xfg, &mut yfg, for_inverse, comp.extent, eps, 0.0) {
                continue;
            }
            xfg = xfg.max(comp.extent.minx).min(comp.extent.maxx);
            yfg = yfg.max(comp.extent.miny).min(comp.extent.maxy);

            let tfactor = comp.time_function.evaluate_at(t);
            if tfactor == 0.0 {
                continue;
            }

            let first_band = grid.bands.first()?;
            let ext = &first_band.extent;
            let rows = ext.rows();
            let cols = ext.cols();
            if cols < 2 || rows < 2 {
                return None;
            }
            let gminx = ext.ll_lon * DEG_TO_RAD;
            let gminy = ext.ll_lat * DEG_TO_RAD;
            let resx = ext.lon_inc * DEG_TO_RAD;
            let resy = ext.lat_inc * DEG_TO_RAD;

            let ix_d = (xfg - gminx) / resx;
            let iy_d = (yfg - gminy) / resy;
            if ix_d < -eps
                || iy_d < -eps
                || ix_d + 1.0 >= cols as f64 + eps
                || iy_d + 1.0 >= rows as f64 + eps
            {
                continue;
            }
            let ix0 = (ix_d.floor() as usize).min(cols - 2);
            let iy0 = (iy_d.floor() as usize).min(rows - 2);
            let ix1 = ix0 + 1;
            let iy1 = iy0 + 1;
            let frct_x = ix_d - ix0 as f64;
            let frct_y = iy_d - iy0 as f64;
            let omfx = 1.0 - frct_x;
            let omfy = 1.0 - frct_y;
            let m00 = omfx * omfy; // SW (ix0, iy0)
            let m10 = frct_x * omfy; // SE (ix1, iy0)
            let m01 = omfx * frct_y; // NW (ix0, iy1)
            let m11 = frct_x * frct_y; // NE (ix1, iy1)

            // Read band `bi` at grid node (ix, iy) where iy is measured from the
            // south. Our `GridBand` rows run north (row 0) to south.
            let val = |bi: usize, ix: usize, iy: usize| -> Option<f64> {
                let band = grid.bands.get(bi)?;
                let row_north = rows - 1 - iy;
                band.get(row_north, ix).map(|v| v as f64)
            };

            match comp.displacement_type {
                DisplacementType::Vertical => {
                    let zidx = if grid.bands.len() == 1 { 0 } else { 2 };
                    let dz00 = val(zidx, ix0, iy0)?;
                    let dz10 = val(zidx, ix1, iy0)?;
                    let dz01 = val(zidx, ix0, iy1)?;
                    let dz11 = val(zidx, ix1, iy1)?;
                    dz += tfactor * (dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11);
                }
                _ if self.is_unit_degree => {
                    // Horizontal / 3D, offsets stored in degrees -> radians.
                    let dx00 = val(0, ix0, iy0)? * DEG_TO_RAD;
                    let dx10 = val(0, ix1, iy0)? * DEG_TO_RAD;
                    let dx01 = val(0, ix0, iy1)? * DEG_TO_RAD;
                    let dx11 = val(0, ix1, iy1)? * DEG_TO_RAD;
                    let dy00 = val(1, ix0, iy0)? * DEG_TO_RAD;
                    let dy10 = val(1, ix1, iy0)? * DEG_TO_RAD;
                    let dy01 = val(1, ix0, iy1)? * DEG_TO_RAD;
                    let dy11 = val(1, ix1, iy1)? * DEG_TO_RAD;
                    if comp.displacement_type == DisplacementType::ThreeD {
                        let dz00 = val(2, ix0, iy0)?;
                        let dz10 = val(2, ix1, iy0)?;
                        let dz01 = val(2, ix0, iy1)?;
                        let dz11 = val(2, ix1, iy1)?;
                        dz += tfactor * (dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11);
                    }
                    dlam += tfactor * (dx00 * m00 + dx01 * m01 + dx10 * m10 + dx11 * m11);
                    dphi += tfactor * (dy00 * m00 + dy01 * m01 + dy10 * m10 + dy11 * m11);
                }
                _ => {
                    // Horizontal / 3D, offsets stored in metres.
                    let de00 = val(0, ix0, iy0)?;
                    let de10 = val(0, ix1, iy0)?;
                    let de01 = val(0, ix0, iy1)?;
                    let de11 = val(0, ix1, iy1)?;
                    let dn00 = val(1, ix0, iy0)?;
                    let dn10 = val(1, ix1, iy0)?;
                    let dn01 = val(1, ix0, iy1)?;
                    let dn11 = val(1, ix1, iy1)?;
                    if comp.displacement_type == DisplacementType::ThreeD {
                        let dz00 = val(2, ix0, iy0)?;
                        let dz10 = val(2, ix1, iy0)?;
                        let dz01 = val(2, ix0, iy1)?;
                        let dz11 = val(2, ix1, iy1)?;
                        dz += tfactor * (dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11);
                    }
                    match comp.interp {
                        InterpMethod::Bilinear => {
                            de += tfactor * (de00 * m00 + de01 * m01 + de10 * m10 + de11 * m11);
                            dn += tfactor * (dn00 * m00 + dn01 * m01 + dn10 * m10 + dn11 * m11);
                        }
                        InterpMethod::GeocentricBilinear => {
                            // Port of PROJ GridEx::getBilinearGeocentric: convert
                            // each corner's E/N offset to a geocentric delta in a
                            // local frame centred on the cell, interpolate, then
                            // convert back to E/N at the query longitude.
                            let y0 = gminy + iy0 as f64 * resy;
                            let sinphi0 = y0.sin();
                            let cosphi0 = y0.cos();
                            let sinphi1 = (y0 + resy).sin();
                            let cosphi1 = (y0 + resy).cos();
                            let sinhalf = (resx / 2.0).sin();
                            let coshalf = (resx / 2.0).cos();

                            let s00 = dn00 * sinphi0;
                            let dx00 = de00 * sinhalf - s00 * coshalf;
                            let dy00 = de00 * coshalf + s00 * sinhalf;
                            let dz00 = dn00 * cosphi0;
                            let s01 = dn01 * sinphi1;
                            let dx01 = de01 * sinhalf - s01 * coshalf;
                            let dy01 = de01 * coshalf + s01 * sinhalf;
                            let dz01 = dn01 * cosphi1;
                            let s10 = dn10 * sinphi0;
                            let dx10 = -de10 * sinhalf - s10 * coshalf;
                            let dy10 = de10 * coshalf - s10 * sinhalf;
                            let dz10 = dn10 * cosphi0;
                            let s11 = dn11 * sinphi1;
                            let dx11 = -de11 * sinhalf - s11 * coshalf;
                            let dy11 = de11 * coshalf - s11 * sinhalf;
                            let dz11 = dn11 * cosphi1;

                            let dxg = m00 * dx00 + m01 * dx01 + m10 * dx10 + m11 * dx11;
                            let dyg = m00 * dy00 + m01 * dy01 + m10 * dy10 + m11 * dy11;
                            let dzg = m00 * dz00 + m01 * dz01 + m10 * dz10 + m11 * dz11;

                            let sinphi = y.sin();
                            let cosphi = y.cos();
                            let lam_rel = (frct_x - 0.5) * resx;
                            let sinlam = lam_rel.sin();
                            let coslam = lam_rel.cos();
                            let de_i = -dxg * sinlam + dyg * coslam;
                            let dn_i = (-dxg * coslam - dyg * sinlam) * sinphi + dzg * cosphi;
                            de += tfactor * de_i;
                            dn += tfactor * dn_i;
                        }
                    }
                }
            }
        }

        if self.is_unit_degree {
            x_out += dlam;
            y_out += dphi;
        } else if self.is_addition {
            let cosphi = y.cos();
            let (d_lam, d_phi) = delta_en_to_lonlat(cosphi, de, dn, self.a, self.b, self.es);
            x_out += d_lam;
            y_out += d_phi;
        } else {
            // Geocentric application of the offset.
            let sinphi = y.sin();
            let cosphi = y.cos();
            let sinlam = x.sin();
            let coslam = x.cos();
            let dnsinphi = dn * sinphi;
            let dxg = -de * sinlam - dnsinphi * coslam;
            let dyg = de * coslam - dnsinphi * sinlam;
            let dzg = dn * cosphi;
            let (mut xg, mut yg, mut zg) = geographic_to_geocentric(x, y, 0.0, self.a, self.es);
            xg += dxg;
            yg += dyg;
            zg += dzg;
            let (lam, phi, _h) = geocentric_to_geographic(xg, yg, zg, self.a, self.es);
            x_out = lam;
            y_out = phi;
        }
        z_out += dz;

        Some((x_out, y_out, z_out))
    }

    /// Port of PROJ `Evaluator::inverse`: fixed-point iteration on `forward`.
    fn inverse(&self, x: f64, y: f64, z: f64, t: f64) -> Option<(f64, f64, f64)> {
        let mut x_out = x;
        let mut y_out = y;
        let mut z_out = z;
        for _ in 0..10 {
            let (xn, yn, zn) = self.forward(x_out, y_out, z_out, t, true)?;
            let dx = xn - x;
            let dy = yn - y;
            let dz = zn - z;
            x_out -= dx;
            y_out -= dy;
            z_out -= dz;
            if dx.abs().max(dy.abs()) < 1e-12 && dz.abs() < 1e-3 {
                return Some((x_out, y_out, z_out));
            }
        }
        None
    }
}

impl Operation for Defmodel {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let t = v[3];
        if !t.is_finite() {
            return Err(ProjError::MissingTime);
        }
        let (x, y, z) = self
            .forward(v[0], v[1], v[2], t, false)
            .ok_or(ProjError::CoordTransfm)?;
        Ok(Coord::new(x, y, z, t))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let t = v[3];
        if !t.is_finite() {
            return Err(ProjError::MissingTime);
        }
        let (x, y, z) = self
            .inverse(v[0], v[1], v[2], t)
            .ok_or(ProjError::CoordTransfm)?;
        Ok(Coord::new(x, y, z, t))
    }

    fn has_inverse(&self) -> bool {
        true
    }
}

/// Parse a grid file from raw bytes, trying the formats a deformation model may
/// reference (GeoTIFF canonically, then NTv2, then GTX).
fn load_grid(data: &[u8], name: &str) -> ProjResult<GridSet> {
    if let Ok(gs) = read_geotiff(data, name) {
        return Ok(gs);
    }
    if let Ok(mut sets) = read_ntv2(data, name) {
        if let Some(gs) = sets.drain(..).next() {
            return Ok(gs);
        }
    }
    read_gtx(data, name)
}

/// Construct a `defmodel` transform.
///
/// Requires `+model=<name>`: the name of the deformation-model master JSON file
/// in the registry. Each component's `spatial_model.filename` is also resolved
/// through the registry.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let model_name = p.params.get_str("model").ok_or(ProjError::MissingArg)?;
    let registry = p.registry.ok_or(ProjError::FileNotFound)?;
    let bytes = registry
        .get_grid(model_name)
        .ok_or(ProjError::FileNotFound)?;
    let text = core::str::from_utf8(bytes).map_err(|_| ProjError::FileNotFound)?;
    let spec = parse_model(text)?;

    let is_unit_degree = spec.h_unit == HUnit::Degree;
    let is_addition = spec.h_method != HMethod::Geocentric;

    let mut components = Vec::with_capacity(spec.components.len());
    for cs in spec.components {
        let grid = if cs.displacement_type == DisplacementType::None {
            None
        } else {
            let gbytes = registry
                .get_grid(&cs.filename)
                .ok_or(ProjError::FileNotFound)?;
            Some(load_grid(gbytes, &cs.filename)?)
        };
        components.push(Comp {
            extent: cs.extent,
            displacement_type: cs.displacement_type,
            interp: cs.interp,
            time_function: cs.time_function,
            grid,
        });
    }

    let op = Defmodel {
        a: p.ellipsoid.a,
        b: p.ellipsoid.b,
        es: p.ellipsoid.es,
        is_unit_degree,
        is_addition,
        extent: spec.extent,
        time_first: spec.time_first,
        time_last: spec.time_last,
        components,
    };

    Ok(TransBuild::new(
        Box::new(op),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use oxiproj_core::Ellipsoid;
    use oxiproj_grids::{GridBand, GridExtent, GridSet, GridSource};

    fn wgs84() -> Ellipsoid {
        Ellipsoid::named("WGS84").unwrap()
    }

    // ---- JSON parser & time-function math (pure, no grids) ----

    fn parse_test_model() -> ModelSpec {
        let json = r#"{
          "file_type": "deformation_model_master_file",
          "format_version": "1.0",
          "name": "test model",
          "source_crs": "EPSG:4959",
          "target_crs": "EPSG:7907",
          "definition_crs": "EPSG:4959",
          "reference_epoch": "2000-01-01T00:00:00Z",
          "vertical_offset_unit": "metre",
          "extent": { "type": "bbox", "parameters": { "bbox": [0.0, 0.0, 1.0, 1.0] } },
          "time_extent": { "first": "1900-01-01T00:00:00Z", "last": "2100-01-01T00:00:00Z" },
          "components": [
            { "description": "constant", "extent": {"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
              "displacement_type": "vertical", "uncertainty_type": "none",
              "spatial_model": {"type":"GeoTIFF","interpolation_method":"bilinear","filename":"a.tif"},
              "time_function": {"type":"constant"} },
            { "description": "velocity", "extent": {"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
              "displacement_type": "vertical", "uncertainty_type": "none",
              "spatial_model": {"type":"GeoTIFF","interpolation_method":"bilinear","filename":"b.tif"},
              "time_function": {"type":"velocity","parameters":{"reference_epoch":"2000-01-01T00:00:00Z"}} },
            { "description": "step", "extent": {"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
              "displacement_type": "vertical", "uncertainty_type": "none",
              "spatial_model": {"type":"GeoTIFF","interpolation_method":"bilinear","filename":"c.tif"},
              "time_function": {"type":"step","parameters":{"step_epoch":"2005-01-01T00:00:00Z"}} },
            { "description": "reverse", "extent": {"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
              "displacement_type": "vertical", "uncertainty_type": "none",
              "spatial_model": {"type":"GeoTIFF","interpolation_method":"bilinear","filename":"d.tif"},
              "time_function": {"type":"reverse_step","parameters":{"step_epoch":"2005-01-01T00:00:00Z"}} },
            { "description": "piecewise", "extent": {"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
              "displacement_type": "vertical", "uncertainty_type": "none",
              "spatial_model": {"type":"GeoTIFF","interpolation_method":"bilinear","filename":"e.tif"},
              "time_function": {"type":"piecewise","parameters":{"before_first":"constant","after_last":"constant",
                "model":[{"epoch":"2000-01-01T00:00:00Z","scale_factor":0.0},
                         {"epoch":"2010-01-01T00:00:00Z","scale_factor":10.0}]}} },
            { "description": "exp", "extent": {"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
              "displacement_type": "vertical", "uncertainty_type": "none",
              "spatial_model": {"type":"GeoTIFF","interpolation_method":"bilinear","filename":"f.tif"},
              "time_function": {"type":"exponential","parameters":{
                "reference_epoch":"2000-01-01T00:00:00Z","relaxation_constant":10.0,
                "before_scale_factor":0.0,"initial_scale_factor":0.0,"final_scale_factor":1.0}} }
          ]
        }"#;
        parse_model(json).expect("model should parse")
    }

    #[test]
    fn test_parse_structure() {
        let m = parse_test_model();
        assert_eq!(m.components.len(), 6);
        assert!(m.h_unit == HUnit::Unspecified); // no horizontal components
        assert!((m.time_first - 1900.0).abs() < 1e-9);
        assert!((m.time_last - 2100.0).abs() < 1e-9);
        for c in &m.components {
            assert_eq!(c.displacement_type, DisplacementType::Vertical);
            assert_eq!(c.interp, InterpMethod::Bilinear);
        }
    }

    #[test]
    fn test_iso8601_decimal_year() {
        assert!((iso8601_to_decimal_year("2000-01-01T00:00:00Z").unwrap() - 2000.0).abs() < 1e-12);
        assert!((iso8601_to_decimal_year("2001-01-01T00:00:00Z").unwrap() - 2001.0).abs() < 1e-12);
        // 2000 is a leap year; 2 July is exactly the middle (day-in-year 183/366).
        assert!((iso8601_to_decimal_year("2000-07-02T00:00:00Z").unwrap() - 2000.5).abs() < 1e-12);
        assert!(iso8601_to_decimal_year("not-a-date").is_err());
    }

    #[test]
    fn test_time_functions() {
        let m = parse_test_model();
        // constant
        assert!((m.components[0].time_function.evaluate_at(2050.0) - 1.0).abs() < 1e-12);
        // velocity (ref 2000)
        assert!((m.components[1].time_function.evaluate_at(2010.0) - 10.0).abs() < 1e-9);
        assert!(m.components[1].time_function.evaluate_at(2000.0).abs() < 1e-9);
        // step (2005)
        assert!(m.components[2].time_function.evaluate_at(2004.0).abs() < 1e-12);
        assert!((m.components[2].time_function.evaluate_at(2005.0) - 1.0).abs() < 1e-12);
        // reverse_step (2005)
        assert!((m.components[3].time_function.evaluate_at(2004.0) + 1.0).abs() < 1e-12);
        assert!(m.components[3].time_function.evaluate_at(2005.0).abs() < 1e-12);
        // piecewise linear midpoint 2005 -> 5.0, clamped constant outside
        assert!((m.components[4].time_function.evaluate_at(2005.0) - 5.0).abs() < 1e-9);
        assert!(m.components[4].time_function.evaluate_at(1990.0).abs() < 1e-9);
        assert!((m.components[4].time_function.evaluate_at(2050.0) - 10.0).abs() < 1e-9);
        // exponential: at ref -> initial (0); at ref+relax -> 1-exp(-1)
        assert!(m.components[5].time_function.evaluate_at(2000.0).abs() < 1e-12);
        let expected = 1.0 - (-1.0f64).exp();
        assert!((m.components[5].time_function.evaluate_at(2010.0) - expected).abs() < 1e-9);
    }

    #[test]
    fn test_json_parser_rejects_deeply_nested_arrays() {
        // A crafted master file with far more nested arrays than
        // MAX_JSON_DEPTH must be rejected gracefully (IllegalArgValue)
        // rather than blowing the native stack.
        let depth = 100_000;
        let mut json = String::with_capacity(depth * 2);
        for _ in 0..depth {
            json.push('[');
        }
        for _ in 0..depth {
            json.push(']');
        }
        assert_eq!(
            JsonParser::parse(&json).err(),
            Some(ProjError::IllegalArgValue)
        );
    }

    #[test]
    fn test_json_parser_rejects_deeply_nested_objects() {
        let depth = 100_000;
        let mut json = String::with_capacity(depth * 10);
        for _ in 0..depth {
            json.push_str(r#"{"a":"#);
        }
        json.push_str("null");
        for _ in 0..depth {
            json.push('}');
        }
        assert_eq!(
            JsonParser::parse(&json).err(),
            Some(ProjError::IllegalArgValue)
        );
    }

    #[test]
    fn test_json_parser_accepts_nesting_within_limit() {
        // Sanity check: nesting up to (but not exceeding) MAX_JSON_DEPTH
        // still parses successfully so the cap doesn't break legitimate,
        // reasonably-structured master files.
        let depth = MAX_JSON_DEPTH;
        let mut json = String::new();
        for _ in 0..depth {
            json.push('[');
        }
        json.push('1');
        for _ in 0..depth {
            json.push(']');
        }
        assert!(JsonParser::parse(&json).is_ok());
    }

    #[test]
    fn test_parse_rejects_source_ne_definition() {
        let json = r#"{
          "file_type":"deformation_model_master_file","format_version":"1.0",
          "source_crs":"EPSG:1","target_crs":"EPSG:2","definition_crs":"EPSG:3",
          "extent":{"type":"bbox","parameters":{"bbox":[0,0,1,1]}},
          "time_extent":{"first":"2000-01-01T00:00:00Z","last":"2001-01-01T00:00:00Z"},
          "components":[]
        }"#;
        assert!(parse_model(json).is_err());
    }

    // ---- Grid helpers ----

    fn uniform_band(value: f32) -> GridBand {
        GridBand {
            extent: GridExtent {
                ll_lat: 0.0,
                ll_lon: 0.0,
                ur_lat: 1.0,
                ur_lon: 1.0,
                lat_inc: 1.0,
                lon_inc: 1.0,
            },
            values: vec![value; 4],
            semantics: Default::default(),
        }
    }

    fn full_extent() -> Extent {
        Extent {
            minx: 0.0,
            miny: 0.0,
            maxx: 1.0 * DEG_TO_RAD,
            maxy: 1.0 * DEG_TO_RAD,
        }
    }

    fn make_op(comp: Comp, is_unit_degree: bool, is_addition: bool) -> Defmodel {
        let ell = wgs84();
        Defmodel {
            a: ell.a,
            b: ell.b,
            es: ell.es,
            is_unit_degree,
            is_addition,
            extent: full_extent(),
            time_first: 1900.0,
            time_last: 2100.0,
            components: vec![comp],
        }
    }

    // ---- Forward / inverse round-trips (direct struct construction) ----

    #[test]
    fn test_horizontal_degree_forward_inverse() {
        // east offset 0.001 deg, north offset 0.002 deg, velocity ref 2000.
        let grid = GridSet {
            bands: vec![uniform_band(0.001), uniform_band(0.002)],
            source: GridSource::Memory,
        };
        let comp = Comp {
            extent: full_extent(),
            displacement_type: DisplacementType::Horizontal,
            interp: InterpMethod::Bilinear,
            time_function: TimeFunction::Velocity {
                reference_epoch: 2000.0,
            },
            grid: Some(grid),
        };
        let op = make_op(comp, true, true);
        let lon = 0.5 * DEG_TO_RAD;
        let lat = 0.5 * DEG_TO_RAD;
        let (x, y, _z) = op.forward(lon, lat, 0.0, 2010.0, false).unwrap();
        // dt = 10 -> dlon = 10*0.001 deg, dlat = 10*0.002 deg.
        // Tolerance accommodates the f32 storage of the grid offsets (0.001/0.002
        // are not exactly representable in f32, giving ~1e-11 absolute error).
        assert!((x - (0.5 + 0.01) * DEG_TO_RAD).abs() < 1e-9, "lon shift");
        assert!((y - (0.5 + 0.02) * DEG_TO_RAD).abs() < 1e-9, "lat shift");
        let (xi, yi, _zi) = op.inverse(x, y, 0.0, 2010.0).unwrap();
        assert!((xi - lon).abs() < 1e-12, "lon round-trip");
        assert!((yi - lat).abs() < 1e-12, "lat round-trip");
    }

    #[test]
    fn test_metre_addition_roundtrip() {
        let grid = GridSet {
            bands: vec![uniform_band(1.0), uniform_band(1.5)],
            source: GridSource::Memory,
        };
        let comp = Comp {
            extent: full_extent(),
            displacement_type: DisplacementType::Horizontal,
            interp: InterpMethod::Bilinear,
            time_function: TimeFunction::Velocity {
                reference_epoch: 2000.0,
            },
            grid: Some(grid),
        };
        let op = make_op(comp, false, true);
        let lon = 0.5 * DEG_TO_RAD;
        let lat = 0.5 * DEG_TO_RAD;
        let (x, y, _z) = op.forward(lon, lat, 0.0, 2010.0, false).unwrap();
        assert!((x - lon).abs() > 0.0, "should shift longitude");
        let (xi, yi, _zi) = op.inverse(x, y, 0.0, 2010.0).unwrap();
        assert!((xi - lon).abs() < 1e-12, "lon round-trip");
        assert!((yi - lat).abs() < 1e-12, "lat round-trip");
    }

    #[test]
    fn test_metre_geocentric_bilinear_roundtrip() {
        let grid = GridSet {
            bands: vec![uniform_band(2.0), uniform_band(-1.0)],
            source: GridSource::Memory,
        };
        let comp = Comp {
            extent: full_extent(),
            displacement_type: DisplacementType::Horizontal,
            interp: InterpMethod::GeocentricBilinear,
            time_function: TimeFunction::Velocity {
                reference_epoch: 2000.0,
            },
            grid: Some(grid),
        };
        let op = make_op(comp, false, true);
        let lon = 0.5 * DEG_TO_RAD;
        let lat = 0.5 * DEG_TO_RAD;
        let (x, y, _z) = op.forward(lon, lat, 0.0, 2010.0, false).unwrap();
        let (xi, yi, _zi) = op.inverse(x, y, 0.0, 2010.0).unwrap();
        assert!((xi - lon).abs() < 1e-11, "lon round-trip");
        assert!((yi - lat).abs() < 1e-11, "lat round-trip");
    }

    #[test]
    fn test_outside_time_extent_is_error() {
        let grid = GridSet {
            bands: vec![uniform_band(0.001), uniform_band(0.002)],
            source: GridSource::Memory,
        };
        let comp = Comp {
            extent: full_extent(),
            displacement_type: DisplacementType::Horizontal,
            interp: InterpMethod::Bilinear,
            time_function: TimeFunction::Constant,
            grid: Some(grid),
        };
        let op = make_op(comp, true, true);
        // t = 3000 is outside [1900, 2100].
        assert!(op
            .forward(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 3000.0, false)
            .is_none());
    }

    // ---- End-to-end through new() + GridRegistry (vertical, GTX grid) ----

    fn build_gtx(value: f32) -> Vec<u8> {
        let mut buf = Vec::new();
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // south_lat = 0
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // west_lon = 0
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lat_inc = 1
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lon_inc = 1
        buf.extend_from_slice(&2i32.to_be_bytes()); // rows = 2
        buf.extend_from_slice(&2i32.to_be_bytes()); // cols = 2
        for _ in 0..4 {
            buf.extend_from_slice(&value.to_be_bytes());
        }
        buf
    }

    struct InMemReg(std::collections::HashMap<String, Vec<u8>>);
    impl crate::GridRegistry for InMemReg {
        fn get_grid(&self, name: &str) -> Option<&[u8]> {
            self.0.get(name).map(|v| v.as_slice())
        }
    }

    struct ModelParam(String);
    impl crate::TransParamLookup for ModelParam {
        fn get_str(&self, key: &str) -> Option<&str> {
            if key == "model" {
                Some(&self.0)
            } else {
                None
            }
        }
        fn get_f64(&self, _: &str) -> Option<f64> {
            None
        }
        fn get_dms(&self, _: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _: &str) -> Option<i64> {
            None
        }
        fn get_bool(&self, _: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            key == "model"
        }
    }

    #[test]
    fn test_end_to_end_vertical_via_registry() {
        let model_json = r#"{
          "file_type":"deformation_model_master_file","format_version":"1.0",
          "source_crs":"EPSG:1","target_crs":"EPSG:1","definition_crs":"EPSG:1",
          "vertical_offset_unit":"metre",
          "extent":{"type":"bbox","parameters":{"bbox":[0.0,0.0,1.0,1.0]}},
          "time_extent":{"first":"1900-01-01T00:00:00Z","last":"2100-01-01T00:00:00Z"},
          "components":[
            { "extent":{"type":"bbox","parameters":{"bbox":[0.0,0.0,1.0,1.0]}},
              "displacement_type":"vertical","uncertainty_type":"none",
              "spatial_model":{"type":"GTX","interpolation_method":"bilinear","filename":"v.gtx"},
              "time_function":{"type":"velocity","parameters":{"reference_epoch":"2000-01-01T00:00:00Z"}} }
          ]
        }"#;
        let mut map = std::collections::HashMap::new();
        map.insert("model.json".to_string(), model_json.as_bytes().to_vec());
        map.insert("v.gtx".to_string(), build_gtx(2.0));
        let reg = InMemReg(map);
        let param = ModelParam("model.json".to_string());
        let ell = wgs84();
        let p = TransParams {
            ellipsoid: &ell,
            params: &param,
            registry: Some(&reg),
        };
        let tb = new(&p).expect("defmodel should build");

        // dt = 10 yr, vertical offset 2.0 m -> dz = 20 m.
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 100.0, 2010.0);
        let out = tb.operation.forward_4d(input).unwrap();
        assert!(
            (out.v()[2] - 120.0).abs() < 1e-6,
            "z forward: {}",
            out.v()[2]
        );
        // Horizontal untouched.
        assert!((out.v()[0] - input.v()[0]).abs() < 1e-12);
        assert!((out.v()[1] - input.v()[1]).abs() < 1e-12);

        let back = tb.operation.inverse_4d(out).unwrap();
        assert!(
            (back.v()[2] - 100.0).abs() < 1e-6,
            "z inverse: {}",
            back.v()[2]
        );
    }

    #[test]
    fn test_missing_model_param_errors() {
        let ell = wgs84();
        struct NoParams;
        impl crate::TransParamLookup for NoParams {
            fn get_str(&self, _: &str) -> Option<&str> {
                None
            }
            fn get_f64(&self, _: &str) -> Option<f64> {
                None
            }
            fn get_dms(&self, _: &str) -> Option<f64> {
                None
            }
            fn get_int(&self, _: &str) -> Option<i64> {
                None
            }
            fn get_bool(&self, _: &str) -> bool {
                false
            }
            fn exists(&self, _: &str) -> bool {
                false
            }
        }
        let np = NoParams;
        let p = TransParams {
            ellipsoid: &ell,
            params: &np,
            registry: None,
        };
        assert_eq!(new(&p).err(), Some(ProjError::MissingArg));
    }
}