sidereon_core/
geodesic.rs

1//! WGS84 geodesic direct and inverse helpers.
2//!
3//! The computations use Karney's geodesic algorithm on the WGS84 ellipsoid.
4//! Public angles are degrees and distances are meters.
5
6use oxiproj_geodesic::Geodesic;
7
8/// Error returned when geodesic inputs are outside the accepted domain.
9#[derive(Debug, Clone, Copy, PartialEq, Eq, thiserror::Error)]
10pub enum GeodesicError {
11    /// A geodesic input was non-finite or outside its documented range.
12    #[error("invalid geodesic input {field}: {reason}")]
13    InvalidInput {
14        /// Name of the invalid field.
15        field: &'static str,
16        /// Reason the field was rejected.
17        reason: &'static str,
18    },
19}
20
21fn invalid_input(field: &'static str, reason: &'static str) -> GeodesicError {
22    GeodesicError::InvalidInput { field, reason }
23}
24
25fn validate_latitude(field: &'static str, value: f64) -> Result<(), GeodesicError> {
26    if !value.is_finite() {
27        return Err(invalid_input(field, "must be finite"));
28    }
29    if !(-90.0..=90.0).contains(&value) {
30        return Err(invalid_input(field, "must be in [-90, 90] degrees"));
31    }
32    Ok(())
33}
34
35fn validate_finite(field: &'static str, value: f64) -> Result<(), GeodesicError> {
36    if value.is_finite() {
37        Ok(())
38    } else {
39        Err(invalid_input(field, "must be finite"))
40    }
41}
42
43/// Solve the WGS84 inverse geodesic problem.
44///
45/// Inputs are point 1 latitude and longitude followed by point 2 latitude and
46/// longitude, all in degrees. Longitudes may be outside `[-180, 180]`; they are
47/// reduced by the geodesic solver. The returned tuple is `(s12_m, azi1_deg,
48/// azi2_deg)`, where `s12_m` is the geodesic distance in meters, `azi1_deg` is
49/// the forward azimuth at point 1, and `azi2_deg` is the forward azimuth at
50/// point 2.
51pub fn geodesic_inverse(
52    lat1_deg: f64,
53    lon1_deg: f64,
54    lat2_deg: f64,
55    lon2_deg: f64,
56) -> Result<(f64, f64, f64), GeodesicError> {
57    validate_latitude("lat1_deg", lat1_deg)?;
58    validate_finite("lon1_deg", lon1_deg)?;
59    validate_latitude("lat2_deg", lat2_deg)?;
60    validate_finite("lon2_deg", lon2_deg)?;
61
62    let geodesic = Geodesic::wgs84();
63    let result = geodesic.inverse(lat1_deg, lon1_deg, lat2_deg, lon2_deg);
64    Ok((result.s12, result.azi1, result.azi2))
65}
66
67/// Solve the WGS84 direct geodesic problem.
68///
69/// Inputs are point 1 latitude, longitude, forward azimuth, and geodesic
70/// distance. Angles are degrees and `s12_m` is meters. Longitudes and azimuths
71/// may be outside one turn; they are reduced by the geodesic solver. The
72/// returned tuple is `(lat2_deg, lon2_deg, azi2_deg)`.
73pub fn geodesic_direct(
74    lat1_deg: f64,
75    lon1_deg: f64,
76    azi1_deg: f64,
77    s12_m: f64,
78) -> Result<(f64, f64, f64), GeodesicError> {
79    validate_latitude("lat1_deg", lat1_deg)?;
80    validate_finite("lon1_deg", lon1_deg)?;
81    validate_finite("azi1_deg", azi1_deg)?;
82    validate_finite("s12_m", s12_m)?;
83
84    let geodesic = Geodesic::wgs84();
85    let result = geodesic.direct(lat1_deg, lon1_deg, azi1_deg, s12_m);
86    Ok((result.lat2, result.lon2, result.azi2))
87}
sidereon_core/geodesic.rs

sidereon_core/
geodesic.rs
