proj-core 0.10.0

use std::f64::consts::FRAC_PI_2;

use crate::ellipsoid::Ellipsoid;
use crate::error::{Error, Result};
use crate::projection::{
    ensure_finite_lon_lat, ensure_finite_xy, normalize_longitude, validate_angle,
    validate_latitude_param, validate_lon_lat, validate_offset, validate_projected, validate_scale,
};

/// Polar Stereographic projection.
///
/// Supports both Variant A (specified by `lat_ts` = latitude of true scale)
/// and Variant B (specified by `k0` = scale factor at the pole).
///
/// Used by EPSG:3413 (NSIDC Sea Ice Polar Stereographic North),
/// EPSG:3031 (Antarctic Polar Stereographic), and others.
#[derive(Clone)]
pub(crate) struct PolarStereographic {
    /// Central meridian / straight vertical longitude (radians).
    lon0: f64,
    /// True if the projection is for the north pole, false for south.
    is_north: bool,
    /// False easting (meters).
    false_easting: f64,
    /// False northing (meters).
    false_northing: f64,
    // Precomputed constants
    e: f64,
    two_a_k0: f64,
    /// t evaluated at the pole, `((1-e)/(1+e))^(e/2)`.
    t_polar: f64,
}

impl PolarStereographic {
    /// Create from parameters. `lat_ts` in radians determines hemisphere.
    /// If `|lat_ts| < 90°`, k0 is computed from lat_ts (Variant A).
    /// If `|lat_ts| ≈ 90°`, the provided k0 is used directly (Variant B).
    pub(crate) fn new(
        ellipsoid: Ellipsoid,
        lon0: f64,
        lat_ts: f64,
        k0_input: f64,
        false_easting: f64,
        false_northing: f64,
    ) -> Result<Self> {
        validate_angle("central meridian", lon0)?;
        validate_latitude_param("latitude of true scale", lat_ts)?;
        validate_scale("scale factor", k0_input)?;
        validate_offset("false easting", false_easting)?;
        validate_offset("false northing", false_northing)?;

        let is_north = lat_ts >= 0.0;
        let e = ellipsoid.e();
        let e2 = ellipsoid.e2();
        let a = ellipsoid.semi_major_axis();

        // Determine k0:
        // If lat_ts is at the pole (±90°), use provided k0.
        // Otherwise compute k0 from lat_ts (Variant A).
        let k0 = if (lat_ts.abs() - FRAC_PI_2).abs() < 1e-10 {
            k0_input
        } else {
            // Compute scale factor from latitude of true scale
            let abs_lat_ts = lat_ts.abs();
            let sin_ts = abs_lat_ts.sin();
            let cos_ts = abs_lat_ts.cos();
            let e_sin_ts = e * sin_ts;

            let m_ts = cos_ts / (1.0 - e2 * sin_ts * sin_ts).sqrt();
            let t_ts = ((FRAC_PI_2 - abs_lat_ts) / 2.0).tan()
                / ((1.0 - e_sin_ts) / (1.0 + e_sin_ts)).powf(e / 2.0);
            let t_90 = compute_t_polar(e);

            m_ts / (2.0 * t_ts / t_90)
        };
        if !k0.is_finite() || k0 <= 0.0 {
            return Err(Error::InvalidDefinition(
                "Polar Stereographic scale factor must resolve to a finite positive number".into(),
            ));
        }

        let two_a_k0 = 2.0 * a * k0;

        Ok(Self {
            lon0,
            is_north,
            false_easting,
            false_northing,
            e,
            two_a_k0,
            t_polar: compute_t_polar(e),
        })
    }

    /// Variant C (EPSG method 9830): the cone geometry of the standard-
    /// parallel form, but with the false origin on the standard parallel
    /// instead of at the pole. Equivalent to the lat_ts form with the false
    /// northing offset by ρF = a·mF, the radius of the standard parallel.
    pub(crate) fn new_variant_c(
        ellipsoid: Ellipsoid,
        lon0: f64,
        lat_ts: f64,
        easting_false_origin: f64,
        northing_false_origin: f64,
    ) -> Result<Self> {
        validate_latitude_param("latitude of standard parallel", lat_ts)?;
        validate_offset("northing at false origin", northing_false_origin)?;
        if (lat_ts.abs() - FRAC_PI_2).abs() < 1e-10 {
            return Err(Error::InvalidDefinition(
                "Polar Stereographic variant C requires a standard parallel off the pole".into(),
            ));
        }
        let e2 = ellipsoid.e2();
        let sin_ts = lat_ts.sin();
        let rho_f =
            ellipsoid.semi_major_axis() * lat_ts.cos() / (1.0 - e2 * sin_ts * sin_ts).sqrt();
        let effective_false_northing = if lat_ts >= 0.0 {
            northing_false_origin + rho_f
        } else {
            northing_false_origin - rho_f
        };
        Self::new(
            ellipsoid,
            lon0,
            lat_ts,
            1.0,
            easting_false_origin,
            effective_false_northing,
        )
    }
}

/// Compute t at the pole: t(90°) = exp(-e * atanh(e)) = ((1-e)/(1+e))^(e/2)
fn compute_t_polar(e: f64) -> f64 {
    ((1.0 - e) / (1.0 + e)).powf(e / 2.0)
}

/// Compute the conformal latitude parameter t for a given latitude.
///
/// Evaluated at the signed latitude: the formula extends continuously across
/// the equator, so opposite-hemisphere inputs map to large-radius points
/// (matching C PROJ) instead of mirroring into the projection's hemisphere.
fn compute_t(lat: f64, e: f64) -> f64 {
    let sin_lat = lat.sin();
    let e_sin = e * sin_lat;
    ((FRAC_PI_2 - lat) / 2.0).tan() / ((1.0 - e_sin) / (1.0 + e_sin)).powf(e / 2.0)
}

impl super::ProjectionImpl for PolarStereographic {
    fn forward(&self, lon: f64, lat: f64) -> Result<(f64, f64)> {
        validate_lon_lat(lon, lat)?;
        let e = self.e;
        let t_polar = self.t_polar;

        // For north polar: use latitude directly.
        // For south polar: negate latitude to work with formulas designed for north.
        let (effective_lat, sign) = if self.is_north {
            (lat, 1.0)
        } else {
            (-lat, -1.0)
        };

        let t = compute_t(effective_lat, e);
        let rho = self.two_a_k0 * t / t_polar;

        let d_lon = normalize_longitude(lon - self.lon0);

        let x = self.false_easting + rho * d_lon.sin();
        let y = self.false_northing - sign * rho * d_lon.cos();

        ensure_finite_xy("Polar Stereographic", x, y)
    }

    fn inverse(&self, x: f64, y: f64) -> Result<(f64, f64)> {
        validate_projected(x, y)?;
        let e = self.e;
        let t_polar = self.t_polar;

        let dx = x - self.false_easting;
        let dy = y - self.false_northing;

        let (dy_eff, sign) = if self.is_north {
            (-dy, 1.0)
        } else {
            (dy, -1.0)
        };

        let rho = (dx * dx + dy_eff * dy_eff).sqrt();

        if rho < 1e-10 {
            // At the pole
            let lat = if self.is_north { FRAC_PI_2 } else { -FRAC_PI_2 };
            return Ok((normalize_longitude(self.lon0), lat));
        }

        let t = rho * t_polar / self.two_a_k0;
        let lat_unsigned =
            super::latitude_from_conformal_t("Polar Stereographic inverse latitude", t, e)?;
        let lat = sign * lat_unsigned;

        let lon = self.lon0 + dx.atan2(dy_eff);

        ensure_finite_lon_lat("Polar Stereographic", lon, lat)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::ellipsoid;
    use crate::projection::ProjectionImpl;

    /// EPSG:3413 — WGS 84 / NSIDC Sea Ice Polar Stereographic North
    /// lat_ts = 70°N, lon0 = -45°, k0 derived from lat_ts
    fn epsg_3413() -> PolarStereographic {
        PolarStereographic::new(
            ellipsoid::WGS84,
            (-45.0_f64).to_radians(),
            70.0_f64.to_radians(),
            1.0,
            0.0,
            0.0,
        )
        .unwrap()
    }

    /// EPSG:3031 — WGS 84 / Antarctic Polar Stereographic
    /// lat_ts = -71°S, lon0 = 0°
    fn epsg_3031() -> PolarStereographic {
        PolarStereographic::new(
            ellipsoid::WGS84,
            0.0_f64.to_radians(),
            (-71.0_f64).to_radians(),
            1.0,
            0.0,
            0.0,
        )
        .unwrap()
    }

    #[test]
    fn north_pole_at_origin() {
        let proj = epsg_3413();
        let (x, y) = proj
            .forward((-45.0_f64).to_radians(), 90.0_f64.to_radians())
            .unwrap();
        assert!(x.abs() < 0.01, "x = {x}");
        assert!(y.abs() < 0.01, "y = {y}");
    }

    /// The EPSG Guidance Note 7-2 worked example for variant C
    /// (EPSG:2985, Petrels 1972 / Terre Adelie Polar Stereographic).
    /// C PROJ does not implement method 9830, so the guidance note is the
    /// authoritative reference.
    #[test]
    fn variant_c_matches_epsg_worked_example() {
        let proj = PolarStereographic::new_variant_c(
            ellipsoid::INTL1924,
            140.0_f64.to_radians(),
            (-67.0_f64).to_radians(),
            300_000.0,
            200_000.0,
        )
        .unwrap();
        // 66°36'18.820"S, 140°04'17.040"E
        let lat: f64 = -(66.0 + 36.0 / 60.0 + 18.820 / 3600.0);
        let lon: f64 = 140.0 + 4.0 / 60.0 + 17.040 / 3600.0;
        let (x, y) = proj.forward(lon.to_radians(), lat.to_radians()).unwrap();
        assert!((x - 303_169.52).abs() < 1e-2, "x = {x}");
        assert!((y - 244_055.72).abs() < 1e-2, "y = {y}");
        let (lon2, lat2) = proj.inverse(303_169.522, 244_055.721).unwrap();
        assert!(
            (lon2.to_degrees() - lon).abs() < 1e-7,
            "{}",
            lon2.to_degrees()
        );
        assert!(
            (lat2.to_degrees() - lat).abs() < 1e-7,
            "{}",
            lat2.to_degrees()
        );
    }

    #[test]
    fn variant_c_rejects_polar_standard_parallel() {
        assert!(PolarStereographic::new_variant_c(
            ellipsoid::WGS84,
            0.0,
            (-90.0_f64).to_radians(),
            0.0,
            0.0
        )
        .is_err());
    }

    #[test]
    fn roundtrip_north() {
        let proj = epsg_3413();
        let lon = (-45.0_f64).to_radians();
        let lat = 75.0_f64.to_radians();
        let (x, y) = proj.forward(lon, lat).unwrap();
        let (lon_back, lat_back) = proj.inverse(x, y).unwrap();
        assert!((lon_back - lon).abs() < 1e-8, "lon: {lon_back} vs {lon}");
        assert!((lat_back - lat).abs() < 1e-8, "lat: {lat_back} vs {lat}");
    }

    #[test]
    fn forward_wraps_longitude_delta() {
        let proj = PolarStereographic::new(
            ellipsoid::WGS84,
            179.0_f64.to_radians(),
            70.0_f64.to_radians(),
            1.0,
            0.0,
            0.0,
        )
        .unwrap();
        let lat = 80.0_f64.to_radians();

        let wrapped = proj.forward((-181.0_f64).to_radians(), lat).unwrap();
        let canonical = proj.forward(179.0_f64.to_radians(), lat).unwrap();

        assert!((wrapped.0 - canonical.0).abs() < 1e-8);
        assert!((wrapped.1 - canonical.1).abs() < 1e-8);
    }

    #[test]
    fn roundtrip_south() {
        let proj = epsg_3031();
        let lon = 45.0_f64.to_radians();
        let lat = (-75.0_f64).to_radians();
        let (x, y) = proj.forward(lon, lat).unwrap();
        let (lon_back, lat_back) = proj.inverse(x, y).unwrap();
        assert!((lon_back - lon).abs() < 1e-8, "lon: {lon_back} vs {lon}");
        assert!((lat_back - lat).abs() < 1e-8, "lat: {lat_back} vs {lat}");
    }

    #[test]
    fn known_value_3413() {
        // Point: 75°N, 0°E in EPSG:3413
        // Verified by roundtrip; exact C PROJ values to be added in reference test suite
        let proj = epsg_3413();
        let lon = 0.0_f64.to_radians();
        let lat = 75.0_f64.to_radians();
        let (x, y) = proj.forward(lon, lat).unwrap();

        // Roundtrip verification
        let (lon_back, lat_back) = proj.inverse(x, y).unwrap();
        assert!(
            (lon_back - lon).abs() < 1e-8,
            "lon roundtrip: {lon_back} vs {lon}"
        );
        assert!(
            (lat_back - lat).abs() < 1e-8,
            "lat roundtrip: {lat_back} vs {lat}"
        );
    }

    #[test]
    fn wrong_hemisphere_extends_continuously() {
        // EPSG:3413-style north polar stereographic (lat_ts 70N, lon0 -45).
        let proj = PolarStereographic::new(
            ellipsoid::WGS84,
            (-45.0_f64).to_radians(),
            70.0_f64.to_radians(),
            1.0,
            0.0,
            0.0,
        )
        .unwrap();

        // The conformal formula is continuous across the equator: a southern
        // point maps to a large-radius coordinate, not a mirrored northern
        // one, and it roundtrips.
        let lon = (-45.0_f64).to_radians();
        let lat = (-30.0_f64).to_radians();
        let (x_south, y_south) = proj.forward(lon, lat).unwrap();
        let (_, y_north) = proj.forward(lon, -lat).unwrap();
        assert!(
            y_south < y_north && y_south < -2.0e7,
            "southern point must not mirror: y_south={y_south}, y_north={y_north}"
        );
        let (lon2, lat2) = proj.inverse(x_south, y_south).unwrap();
        assert!((lon2 - lon).abs() < 1e-12, "lon: {lon2}");
        assert!((lat2 - lat).abs() < 1e-12, "lat: {lat2}");

        // Continuity across the equator itself.
        let eps = 1e-9;
        let (_, y_plus) = proj.forward(lon, eps).unwrap();
        let (_, y_minus) = proj.forward(lon, -eps).unwrap();
        assert!(
            (y_plus - y_minus).abs() < 1.0,
            "equator crossing must be continuous: {y_plus} vs {y_minus}"
        );
    }
}