oxiproj-engine 0.1.1

//! The `Pj` operation wrapper: prepare/finalize around an inner operation.
//!
//! Ported from PROJ 9.8.0 `src/fwd.cpp` (`fwd_prepare`/`fwd_finalize`) and
//! `src/inv.cpp` (`inv_prepare`/`inv_finalize`).

use std::sync::Arc;

use oxiproj_core::{
    Coord, Ellipsoid, IoUnits, Lp, Operation, ProjError, ProjResult, ProjectGenericBox, Xy,
};

/// Latitude tolerance for the prepare-stage domain check (radians).
const PJ_EPS_LAT: f64 = 1e-12;

/// A prepared coordinate operation with input/output unit handling.
///
/// Wraps an inner [`Operation`] and applies the unit conversions, prime-meridian
/// and central-meridian shifts, and axis bookkeeping that PROJ performs around
/// the raw projection math.
#[derive(Debug)]
pub struct Pj {
    /// The inner operation (projection or conversion) driven via 4D.
    pub operation: Box<dyn Operation>,
    /// The reference ellipsoid.
    pub ellipsoid: Ellipsoid,
    /// Central meridian (radians).
    pub lam0: f64,
    /// Latitude of origin (radians).
    pub phi0: f64,
    /// False easting (meters, pre-scale).
    pub x0: f64,
    /// False northing (meters, pre-scale).
    pub y0: f64,
    /// False up (meters, pre-scale).
    pub z0: f64,
    /// Scale factor at the central meridian.
    pub k0: f64,
    /// Horizontal input-unit-to-meter factor.
    pub to_meter: f64,
    /// Horizontal meter-to-output-unit factor.
    pub fr_meter: f64,
    /// Vertical input-unit-to-meter factor.
    pub vto_meter: f64,
    /// Vertical meter-to-output-unit factor.
    pub vfr_meter: f64,
    /// Prime-meridian offset from Greenwich (radians).
    pub from_greenwich: f64,
    /// Allow longitudes outside the normal range (`+over`).
    pub over: bool,
    /// Geocentric latitude flag (`+geoc`).
    pub geoc: bool,
    /// Whether this is a pass-through geographic (lat/long) operation.
    pub is_latlong: bool,
    /// Input units (forward direction).
    pub left: IoUnits,
    /// Output units (forward direction).
    pub right: IoUnits,
    /// Whether the operation runs inverted in its containing pipeline.
    pub inverted: bool,
    /// Skip prepare/finalize entirely (used by the top-level pipeline driver).
    pub bypass_prepare_finalize: bool,
    /// Skip this step when the pipeline runs in the forward direction.
    pub omit_fwd: bool,
    /// Skip this step when the pipeline runs in the inverse direction.
    pub omit_inv: bool,
    /// Optional AD-capable projection for exact Jacobian computation (T7.1).
    ///
    /// When `Some`, `factors_exact` uses this to drive forward-mode AD via
    /// `Dual1<2>`. Set to `None` for projections that do not implement
    /// `ProjectGeneric` or whose math is not differentiable.
    pub ad_proj: Option<Arc<dyn ProjectGenericBox>>,
    /// The operation name (e.g. `"helmert"`, `"axisswap"`, `"noop"`).
    ///
    /// Populated from the `+proj=` key at construction time. Used by the
    /// pipeline optimizer (`pipeline_opt`) to detect and eliminate
    /// redundant steps. Empty for the outer pipeline wrapper itself.
    pub op_name: String,
}

impl Pj {
    /// Forward transform with prepare/finalize.
    ///
    /// Ported from `src/fwd.cpp` (`fwd_prepare` -> operation -> `fwd_finalize`).
    fn forward_impl(&self, c: Coord) -> ProjResult<Coord> {
        // --- fwd_prepare ---
        let v = c.v();
        if !v[0].is_finite() || !v[1].is_finite() || !v[2].is_finite() {
            return Ok(Coord::error());
        }

        let prepared = if self.left == IoUnits::Radians {
            let mut lam = v[0];
            let mut phi = v[1];
            let z = v[2];
            let t = v[3];

            let tt = phi.abs() - oxiproj_core::M_HALFPI;
            if tt > PJ_EPS_LAT {
                return Err(ProjError::InvalidCoord);
            }
            if !(-10.0..=10.0).contains(&lam) {
                return Err(ProjError::InvalidCoord);
            }
            // Clamp latitude to ±π/2 (PROJ src/fwd.cpp).
            phi = phi.clamp(-oxiproj_core::M_HALFPI, oxiproj_core::M_HALFPI);

            if self.geoc && self.ellipsoid.es != 0.0 {
                phi = (self.ellipsoid.one_es * phi.tan()).atan();
            }
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            // Vertical grid shift (vgridshift) in prepare is not yet applied here.
            lam = (lam - self.from_greenwich) - self.lam0;
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            Coord::new(lam, phi, z, t)
        } else {
            c
        };

        // --- operation ---
        let mid = self.operation.forward_4d(prepared)?;

        // --- fwd_finalize ---
        let mv = mid.v();
        let mut x = mv[0];
        let mut y = mv[1];
        let mut z = mv[2];
        let t = mv[3];
        match self.right {
            IoUnits::Classic | IoUnits::Projected => {
                if self.right == IoUnits::Classic {
                    x *= self.ellipsoid.a;
                    y *= self.ellipsoid.a;
                }
                x = self.fr_meter * (x + self.x0);
                y = self.fr_meter * (y + self.y0);
                z = self.vfr_meter * (z + self.z0);
            }
            IoUnits::Radians => {
                z = self.vfr_meter * (z + self.z0);
                if !self.over {
                    x = oxiproj_core::adjlon(x);
                }
            }
            IoUnits::Cartesian => {
                x *= self.fr_meter;
                y *= self.fr_meter;
                z *= self.fr_meter;
                // Note: is_geocent cartesian handling is not yet applied here.
            }
            IoUnits::Whatever | IoUnits::Degrees => {}
        }
        // Note: axisswap in finalize is not yet applied here.
        Ok(Coord::new(x, y, z, t))
    }

    /// Inverse transform with prepare/finalize.
    ///
    /// Ported from `src/inv.cpp` (`inv_prepare` -> operation -> `inv_finalize`).
    fn inverse_impl(&self, c: Coord) -> ProjResult<Coord> {
        // --- inv_prepare ---
        let v = c.v();
        if !v[0].is_finite() || !v[1].is_finite() || !v[2].is_finite() {
            return Err(ProjError::OutsideProjectionDomain);
        }
        // Note: axisswap in inv_prepare is not yet applied here.
        let t = v[3];
        let prepared = match self.right {
            IoUnits::Whatever | IoUnits::Degrees => c,
            IoUnits::Cartesian => {
                let x = self.to_meter * v[0];
                let y = self.to_meter * v[1];
                let z = self.to_meter * v[2];
                // Note: is_geocent cartesian handling is not yet applied here.
                Coord::new(x, y, z, t)
            }
            IoUnits::Projected | IoUnits::Classic => {
                let mut x = self.to_meter * v[0] - self.x0;
                let mut y = self.to_meter * v[1] - self.y0;
                let z = self.vto_meter * v[2] - self.z0;
                if self.right == IoUnits::Classic {
                    x *= self.ellipsoid.ra;
                    y *= self.ellipsoid.ra;
                }
                Coord::new(x, y, z, t)
            }
            IoUnits::Radians => {
                let z = self.vto_meter * v[2] - self.z0;
                Coord::new(v[0], v[1], z, t)
            }
        };

        // --- operation ---
        let mid = self.operation.inverse_4d(prepared)?;

        // --- inv_finalize ---
        if self.left == IoUnits::Radians {
            let mv = mid.v();
            let mut lam = mv[0];
            let phi = mv[1];
            let z = mv[2];
            let t = mv[3];
            lam = lam + self.from_greenwich + self.lam0;
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            let phi = if self.geoc && self.ellipsoid.es != 0.0 {
                let limit = oxiproj_core::M_HALFPI - 1e-9;
                if phi.abs() < limit {
                    (self.ellipsoid.rone_es * phi.tan()).atan()
                } else {
                    phi
                }
            } else {
                phi
            };
            // Note: vgridshift in inv_finalize is not yet applied here; datum shift handled via pipeline injection in create.rs.
            Ok(Coord::new(lam, phi, z, t))
        } else {
            Ok(mid)
        }
    }

    /// Run the operation in the forward direction, honoring `inverted`.
    pub fn forward(&self, c: Coord) -> ProjResult<Coord> {
        if self.bypass_prepare_finalize {
            if self.inverted {
                return self.operation.inverse_4d(c);
            } else {
                return self.operation.forward_4d(c);
            }
        }
        if self.inverted {
            self.inverse_impl(c)
        } else {
            self.forward_impl(c)
        }
    }

    /// Run the operation in the inverse direction, honoring `inverted`.
    pub fn inverse(&self, c: Coord) -> ProjResult<Coord> {
        if self.bypass_prepare_finalize {
            if self.inverted {
                return self.operation.forward_4d(c);
            } else {
                return self.operation.inverse_4d(c);
            }
        }
        if self.inverted {
            self.forward_impl(c)
        } else {
            self.inverse_impl(c)
        }
    }

    /// Compute map-scale factors (Tissot indicatrix) at a geographic coordinate.
    ///
    /// `coord` must be in geographic input units — radians if `self.left == IoUnits::Radians`,
    /// degrees if `self.left == IoUnits::Degrees`. The coordinate layout is
    /// `Coord::new(longitude, latitude, z, t)`.
    ///
    /// Returns the [`oxiproj_core::Factors`] struct describing Tissot indicatrix parameters
    /// at the given point.
    ///
    /// Ported from PROJ 9.8.0 `src/factors.cpp` (`proj_factors` → `pj_factors` → `pj_deriv`).
    pub fn factors(&self, coord: Coord) -> ProjResult<oxiproj_core::Factors> {
        let v = coord.v();
        // Convert lon/lat to radians
        let (phi, lam_abs) = if self.left == IoUnits::Radians {
            (v[1], v[0])
        } else {
            // Assume degrees
            (v[1].to_radians(), v[0].to_radians())
        };

        // Longitude relative to central meridian (as done in fwd_prepare)
        let lam = lam_abs - self.from_greenwich - self.lam0;

        let es = self.ellipsoid.es;
        let one_es = self.ellipsoid.one_es;

        // The closure calls the raw inner operation's forward in normalized space.
        let op = &*self.operation;

        oxiproj_core::Factors::compute(phi, lam, es, one_es, |lam_f, phi_f| {
            let c = Coord::new(lam_f, phi_f, 0.0, 0.0);
            let r = op.forward_4d(c)?;
            let rv = r.v();
            Ok((rv[0], rv[1]))
        })
    }

    /// Compute exact distortion factors via automatic differentiation.
    ///
    /// Returns `Err(ProjError::UnsupportedOperation)` when:
    /// - `ad_proj` is `None` (no AD-capable projection registered), or
    /// - the projection's `project_fwd_generic` returns `UnsupportedOperation`.
    ///
    /// `coord` must be in the same geographic input units as for [`Pj::factors`].
    pub fn factors_exact(&self, coord: Coord) -> ProjResult<oxiproj_core::FactorsExact> {
        use oxiproj_core::autodiff::Dual1;

        let ad = match &self.ad_proj {
            None => return Err(ProjError::UnsupportedOperation),
            Some(p) => p,
        };

        let v = coord.v();
        let (phi, lam_abs) = if self.left == IoUnits::Radians {
            (v[1], v[0])
        } else {
            (v[1].to_radians(), v[0].to_radians())
        };
        let lam = lam_abs - self.from_greenwich - self.lam0;

        let lam_d = Dual1::<2>::variable(lam, 0);
        let phi_d = Dual1::<2>::variable(phi, 1);
        let (x, y) = ad.project_fwd_dual2(lam_d, phi_d)?;

        let a = self.ellipsoid.a;
        Ok(oxiproj_core::FactorsExact::from_jacobian(
            x.d[0], x.d[1], y.d[0], y.d[1], phi, a,
        ))
    }
}

impl Operation for Pj {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        self.forward(c)
    }
    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        self.inverse(c)
    }
    fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
        let c = self.forward(Coord::new(lp.lam, lp.phi, 0.0, 0.0))?;
        Ok(c.xy())
    }
    fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
        let c = self.inverse(Coord::new(xy.x, xy.y, 0.0, 0.0))?;
        Ok(c.lp())
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[derive(Debug)]
    struct Id;
    impl Operation for Id {
        fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
            Ok(Xy::new(lp.lam, lp.phi))
        }
        fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
            Ok(Lp::new(xy.x, xy.y))
        }
    }

    fn unit_latlong() -> Pj {
        Pj {
            operation: Box::new(Id),
            ellipsoid: Ellipsoid::named("WGS84").unwrap(),
            lam0: 0.0,
            phi0: 0.0,
            x0: 0.0,
            y0: 0.0,
            z0: 0.0,
            k0: 1.0,
            to_meter: 1.0,
            fr_meter: 1.0,
            vto_meter: 1.0,
            vfr_meter: 1.0,
            from_greenwich: 0.0,
            over: false,
            geoc: false,
            is_latlong: true,
            left: IoUnits::Radians,
            right: IoUnits::Radians,
            inverted: false,
            bypass_prepare_finalize: false,
            omit_fwd: false,
            omit_inv: false,
            ad_proj: None,
            op_name: String::new(),
        }
    }

    #[test]
    fn latlong_round_trip() {
        let pj = unit_latlong();
        let lam = 0.2;
        let phi = 0.5;
        let fwd = pj.forward(Coord::new(lam, phi, 0.0, 0.0)).unwrap();
        let inv = pj.inverse(fwd).unwrap();
        let r = inv.v();
        assert!((r[0] - lam).abs() < 1e-12);
        assert!((r[1] - phi).abs() < 1e-12);
    }

    #[test]
    fn non_finite_forward_yields_error_coord() {
        let pj = unit_latlong();
        let out = pj.forward(Coord::new(f64::NAN, 0.0, 0.0, 0.0)).unwrap();
        assert!(out.is_error());
    }

    #[test]
    fn out_of_domain_latitude_rejected() {
        let pj = unit_latlong();
        let res = pj.forward(Coord::new(0.0, 2.0, 0.0, 0.0));
        assert!(res.is_err());
    }

    #[test]
    fn factors_exact_spherical_mercator_analytic() {
        use oxiproj_core::autodiff::Dual1;

        // Spherical Mercator Jacobian at lam=0.3, phi=0.5:
        //   x = k0*lam  => dx/dlam = k0, dx/dphi = 0
        //   y = k0*ln(tan(pi/4 + phi/2)) = k0*asinh(tan(phi))
        //      => dy/dlam = 0, dy/dphi = k0/cos(phi) = k0*sec(phi)
        let phi0 = 0.5_f64;
        let lam_d = Dual1::<2>::variable(0.3, 0);
        let phi_d = Dual1::<2>::variable(phi0, 1);

        let k0 = 1.0_f64;
        // x = k0 * lam
        let x = lam_d * Dual1::constant(k0);
        // y = k0 * asinh(tan(phi)) = k0 * ln(tan(phi) + sqrt(tan^2(phi) + 1))
        let t = phi_d.tan();
        let y = (t + (t * t + Dual1::constant(1.0)).sqrt()).ln() * Dual1::constant(k0);

        // Verify Jacobian entries
        assert!((x.d[0] - k0).abs() < 1e-12, "dx/dlam = k0");
        assert!(x.d[1].abs() < 1e-12, "dx/dphi = 0");
        assert!(y.d[0].abs() < 1e-12, "dy/dlam = 0");
        let sec_phi = 1.0 / phi0.cos();
        assert!(
            (y.d[1] - sec_phi * k0).abs() < 1e-9,
            "dy/dphi = sec(phi)*k0"
        );

        // FactorsExact: for spherical Mercator, h = k = sec(phi), omega = 0.
        // Mercator is conformal (h=k, omega=0) but NOT equal-area.
        // s = |det J| / (a^2 * cos phi) = sec(phi) / cos(phi) = sec^2(phi).
        let fe =
            oxiproj_core::FactorsExact::from_jacobian(x.d[0], x.d[1], y.d[0], y.d[1], phi0, 1.0);
        assert!(
            (fe.meridian_scale - sec_phi).abs() < 1e-9,
            "h={}",
            fe.meridian_scale
        );
        assert!(
            (fe.parallel_scale - sec_phi).abs() < 1e-9,
            "k={}",
            fe.parallel_scale
        );
        let expected_s = sec_phi * sec_phi;
        assert!(
            (fe.areal_scale - expected_s).abs() < 1e-9,
            "s={}",
            fe.areal_scale
        );
        assert!(
            fe.angular_distortion.abs() < 1e-9,
            "omega={}",
            fe.angular_distortion
        );
        assert!(fe.is_exact);
    }

    #[test]
    fn factors_exact_no_ad_proj_returns_error() {
        let pj = unit_latlong();
        let res = pj.factors_exact(Coord::new(0.1, 0.3, 0.0, 0.0));
        assert!(
            matches!(res, Err(oxiproj_core::ProjError::UnsupportedOperation)),
            "expected UnsupportedOperation, got {res:?}"
        );
    }
}