rustial_engine/math/
globe.rs

1//! Globe / geocentric projection.
2//!
3//! This projection maps WGS-84 geographic coordinates directly onto a 3D globe
4//! centered at the origin. It is useful as the geometric foundation for globe
5//! rendering and globe-aware camera helpers.
6
7use crate::bounds::GeoBounds;
8use crate::coord::{GeoCoord, WorldCoord};
9use crate::ellipsoid::Ellipsoid;
10use crate::projection::Projection;
11
12/// Geocentric globe projection using the WGS-84 ellipsoid.
13///
14/// Output coordinates are Earth-centered Earth-fixed style Cartesian meters:
15///
16/// - `x`: intersects the equator at longitude 0
17/// - `y`: intersects the equator at longitude 90E
18/// - `z`: north pole axis
19pub struct Globe;
20
21impl Globe {
22    /// Project geographic coordinates to globe-centered Cartesian meters.
23    #[inline]
24    pub fn project(geo: &GeoCoord) -> WorldCoord {
25        <Self as Projection>::project(&Self, geo)
26    }
27
28    /// Inverse-project globe-centered Cartesian meters to geographic coordinates.
29    #[inline]
30    pub fn unproject(world: &WorldCoord) -> GeoCoord {
31        <Self as Projection>::unproject(&Self, world)
32    }
33
34    /// Mean Earth radius used for simple scale queries.
35    #[inline]
36    pub fn radius() -> f64 {
37        Ellipsoid::WGS84.a
38    }
39}
40
41impl Projection for Globe {
42    fn project(&self, geo: &GeoCoord) -> WorldCoord {
43        let ellipsoid = Ellipsoid::WGS84;
44        let lat = geo.lat.to_radians();
45        let lon = geo.lon.to_radians();
46        let (sin_lat, cos_lat) = lat.sin_cos();
47        let (sin_lon, cos_lon) = lon.sin_cos();
48
49        let a = ellipsoid.a;
50        let e2 = ellipsoid.e2;
51        let n = a / (1.0_f64 - e2 * sin_lat * sin_lat).sqrt();
52        let alt = geo.alt;
53
54        let x = (n + alt) * cos_lat * cos_lon;
55        let y = (n + alt) * cos_lat * sin_lon;
56        let z = (n * (1.0 - e2) + alt) * sin_lat;
57        WorldCoord::new(x, y, z)
58    }
59
60    fn unproject(&self, world: &WorldCoord) -> GeoCoord {
61        let ellipsoid = Ellipsoid::WGS84;
62        let a = ellipsoid.a;
63        let b = ellipsoid.b;
64        let e2 = ellipsoid.e2;
65        let ep2 = (a * a - b * b) / (b * b);
66
67        let x = world.position.x;
68        let y = world.position.y;
69        let z = world.position.z;
70
71        let lon = y.atan2(x);
72        let p = (x * x + y * y).sqrt();
73
74        if p <= 1e-12 {
75            let lat = if z >= 0.0 {
76                std::f64::consts::FRAC_PI_2
77            } else {
78                -std::f64::consts::FRAC_PI_2
79            };
80            let alt = z.abs() - b;
81            return GeoCoord::new(lat.to_degrees(), 0.0, alt);
82        }
83
84        let theta = (z * a).atan2(p * b);
85        let (sin_theta, cos_theta) = theta.sin_cos();
86        let lat = (z + ep2 * b * sin_theta.powi(3)).atan2(p - e2 * a * cos_theta.powi(3));
87        let sin_lat = lat.sin();
88        let n = a / (1.0_f64 - e2 * sin_lat * sin_lat).sqrt();
89        let alt = p / lat.cos() - n;
90
91        GeoCoord::new(lat.to_degrees(), lon.to_degrees(), alt)
92    }
93
94    fn scale_factor(&self, _geo: &GeoCoord) -> f64 {
95        1.0
96    }
97
98    fn projection_bounds(&self) -> GeoBounds {
99        GeoBounds::new(
100            GeoCoord::from_lat_lon(-90.0, -180.0),
101            GeoCoord::from_lat_lon(90.0, 180.0),
102        )
103    }
104}
105
106#[cfg(test)]
107mod tests {
108    use super::*;
109
110    #[test]
111    fn roundtrip_origin() {
112        let geo = GeoCoord::from_lat_lon(0.0, 0.0);
113        let world = Globe::project(&geo);
114        let back = Globe::unproject(&world);
115        assert!((back.lat - geo.lat).abs() < 1e-8);
116        assert!((back.lon - geo.lon).abs() < 1e-8);
117    }
118
119    #[test]
120    fn roundtrip_nonzero() {
121        let geo = GeoCoord::new(51.1, 17.0, 1200.0);
122        let world = Globe::project(&geo);
123        let back = Globe::unproject(&world);
124        assert!((back.lat - geo.lat).abs() < 1e-6);
125        assert!((back.lon - geo.lon).abs() < 1e-6);
126        assert!((back.alt - geo.alt).abs() < 1e-3);
127    }
128
129    #[test]
130    fn equator_axes_match_expected_orientation() {
131        let zero = Globe::project(&GeoCoord::from_lat_lon(0.0, 0.0));
132        let east = Globe::project(&GeoCoord::from_lat_lon(0.0, 90.0));
133        let north = Globe::project(&GeoCoord::from_lat_lon(90.0, 0.0));
134
135        assert!(zero.position.x > 6_000_000.0);
136        assert!(east.position.y > 6_000_000.0);
137        assert!(north.position.z > 6_000_000.0);
138    }
139}
rustial_engine/math/globe.rs

rustial_engine/math/
globe.rs
