rustial-engine 0.0.1

//! Equirectangular (Plate Carree) projection (EPSG:4326 scaled to meters).
//!
//! # Overview
//!
//! The Equirectangular projection is the simplest cylindrical projection:
//! longitude and latitude are mapped linearly to X and Y, then scaled by
//! the WGS-84 semi-major axis so that output units are **meters** at the
//! equator.
//!
//! ```text
//! x = R * lon_rad
//! y = R * lat_rad
//! ```
//!
//! # Distortion
//!
//! - **Equator:** no distortion (scale factor = 1.0).
//! - **Away from equator:** X distances are progressively stretched by
//!   `sec(lat)`.  Y distances remain undistorted at all latitudes.
//! - **Poles:** the projection is valid up to +/-90 degrees latitude,
//!   but X distortion is infinite at the poles.
//!
//! # When to use
//!
//! Equirectangular is useful as a simple baseline, for data interchange
//! (many raster datasets are stored in EPSG:4326 pixel grids), and for
//! server-side tile pre-processing where rendering fidelity is not
//! critical.  For map display, prefer [`WebMercator`](crate::WebMercator).
//!
//! # Altitude
//!
//! Altitude (`alt` / `z`) is passed through unchanged in both directions.

use crate::coord::{GeoCoord, WorldCoord};
use crate::ellipsoid::Ellipsoid;
use crate::projection::Projection;

/// Equirectangular projection: trivial linear mapping of lon/lat to meters.
///
/// Uses the WGS-84 semi-major axis (`R = 6 378 137 m`) as the scaling
/// factor so that output units are meters at the equator.  This is a
/// zero-size type -- construction is free.
pub struct Equirectangular;

impl Projection for Equirectangular {
    /// Forward projection: geographic degrees to world meters.
    ///
    /// `x = R * lon_rad`, `y = R * lat_rad`, `z = alt` (passthrough).
    fn project(&self, geo: &GeoCoord) -> WorldCoord {
        let a = Ellipsoid::WGS84.a;
        WorldCoord::new(a * geo.lon.to_radians(), a * geo.lat.to_radians(), geo.alt)
    }

    /// Inverse projection: world meters back to geographic degrees.
    ///
    /// `lat = degrees(y / R)`, `lon = degrees(x / R)`, `alt = z`.
    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
        let a = Ellipsoid::WGS84.a;
        let lat = (world.position.y / a).to_degrees().clamp(-90.0, 90.0);
        let lon = ((world.position.x / a).to_degrees() + 180.0).rem_euclid(360.0) - 180.0;
        GeoCoord::new(lat, lon, world.position.z)
    }

    /// Scale factor for Equirectangular: `sec(lat) = 1 / cos(lat)`.
    ///
    /// The projection maps longitude linearly, so X distances are
    /// stretched relative to true distances by `sec(lat)`, the same
    /// as the meridian convergence factor.  Along Y the scale is 1.0,
    /// so this returns the *maximum* linear distortion at the point.
    fn scale_factor(&self, geo: &GeoCoord) -> f64 {
        1.0 / geo.lat.to_radians().cos()
    }

    // `projection_bounds` is not overridden: the default full-globe
    // range (+/-90 lat, +/-180 lon) is correct for Equirectangular.
}

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

    // -- Roundtrip ---------------------------------------------------------

    #[test]
    fn roundtrip_origin() {
        let geo = GeoCoord::from_lat_lon(0.0, 0.0);
        let world = Equirectangular.project(&geo);
        let back = Equirectangular.unproject(&world);
        assert!((back.lat - geo.lat).abs() < 1e-10);
        assert!((back.lon - geo.lon).abs() < 1e-10);
    }

    #[test]
    fn roundtrip_nonzero() {
        let geo = GeoCoord::from_lat_lon(45.0, 90.0);
        let world = Equirectangular.project(&geo);
        let back = Equirectangular.unproject(&world);
        assert!((back.lat - geo.lat).abs() < 1e-8);
        assert!((back.lon - geo.lon).abs() < 1e-8);
    }

    #[test]
    fn poles_roundtrip() {
        let north = GeoCoord::from_lat_lon(90.0, 0.0);
        let north_w = Equirectangular.project(&north);
        let north_back = Equirectangular.unproject(&north_w);
        assert!((north_back.lat - 90.0).abs() < 1e-9);

        let south = GeoCoord::from_lat_lon(-90.0, 0.0);
        let south_w = Equirectangular.project(&south);
        let south_back = Equirectangular.unproject(&south_w);
        assert!((south_back.lat + 90.0).abs() < 1e-9);
    }

    #[test]
    fn anti_meridian_roundtrip() {
        // +180 and -180 represent the same meridian; the normalisation
        // in `unproject` canonicalises both to the (-180, 180] range.
        let east = GeoCoord::from_lat_lon(10.0, 180.0);
        let east_back = Equirectangular.unproject(&Equirectangular.project(&east));
        assert!((east_back.lon.abs() - 180.0).abs() < 1e-9);

        let west = GeoCoord::from_lat_lon(10.0, -180.0);
        let west_back = Equirectangular.unproject(&Equirectangular.project(&west));
        assert!((west_back.lon.abs() - 180.0).abs() < 1e-9);
    }

    #[test]
    fn extreme_latitudes_stable() {
        for lat in [89.999, -89.999, 89.9, -89.9] {
            let geo = GeoCoord::from_lat_lon(lat, 45.0);
            let world = Equirectangular.project(&geo);
            let back = Equirectangular.unproject(&world);
            assert!(
                (back.lat - lat).abs() < 1e-6,
                "lat roundtrip failed for {lat}"
            );
        }
    }

    // -- Metric output at the equator -------------------------------------

    #[test]
    fn equator_scale() {
        let a = GeoCoord::from_lat_lon(0.0, 0.0);
        let b = GeoCoord::from_lat_lon(0.0, 1.0);
        let wa = Equirectangular.project(&a);
        let wb = Equirectangular.project(&b);
        let dx = wb.position.x - wa.position.x;
        // 1 degree of longitude at the equator ~ 111,319.49 m.
        assert!((dx - 111_319.49).abs() < 1.0);
    }

    // -- Scale factor ------------------------------------------------------

    #[test]
    fn scale_factor_equator() {
        let sf = Equirectangular.scale_factor(&GeoCoord::from_lat_lon(0.0, 0.0));
        assert!((sf - 1.0).abs() < 1e-10);
    }

    #[test]
    fn scale_factor_45_degrees() {
        // sec(45 deg) = sqrt(2)
        let sf = Equirectangular.scale_factor(&GeoCoord::from_lat_lon(45.0, 0.0));
        assert!((sf - std::f64::consts::SQRT_2).abs() < 1e-10);
    }

    // -- Altitude passthrough ----------------------------------------------

    #[test]
    fn altitude_passthrough() {
        let geo = GeoCoord::new(30.0, 60.0, 1234.5);
        let world = Equirectangular.project(&geo);
        let back = Equirectangular.unproject(&world);
        assert!((world.position.z - 1234.5).abs() < 1e-12);
        assert!((back.alt - 1234.5).abs() < 1e-12);
    }

    // -- Projection bounds default to full globe --------------------------

    #[test]
    fn projection_bounds_full_globe() {
        let bounds = Equirectangular.projection_bounds();
        assert!((bounds.sw().lat - (-90.0)).abs() < 1e-10);
        assert!((bounds.ne().lat - 90.0).abs() < 1e-10);
        assert!((bounds.sw().lon - (-180.0)).abs() < 1e-10);
        assert!((bounds.ne().lon - 180.0).abs() < 1e-10);
    }
}