initial commit

src/coordinates.rs

@@ -0,0 +1,192 @@
+use geojson::{Position, Value};
+use std::convert::From;
+use std::f64::consts::{FRAC_PI_2, PI, TAU};
+use crate::ACCURACY_BOUNDARY;
+
+/// Spherical coordinates in radians.
+#[derive(Clone, Copy, Debug, PartialEq, Default)]
+pub struct RadianCoordinate {
+    pub lat: f64,
+    pub lon: f64,
+}
+
+/// Spherical coordinates in degrees.
+#[derive(Clone, Copy, Debug, PartialEq, Default)]
+pub struct DegreeCoordinate {
+    pub lat: f64,
+    pub lon: f64,
+}
+
+impl From<DegreeCoordinate> for RadianCoordinate {
+    fn from(other: DegreeCoordinate) -> Self {
+        RadianCoordinate::from_degrees(other.lat, other.lon)
+    }
+}
+
+impl From<DegreeCoordinate> for geojson::Value {
+    fn from(coordinate: DegreeCoordinate) -> Self {
+        Value::Point(vec![coordinate.lat, coordinate.lon])
+    }
+}
+
+impl From<RadianCoordinate> for Position {
+    fn from(coordinate: RadianCoordinate) -> Self {
+        let coordinate = DegreeCoordinate::from(coordinate);
+        vec![coordinate.lat, coordinate.lon]
+    }
+}
+
+impl From<RadianCoordinate> for geojson::Value {
+    fn from(coordinate: RadianCoordinate) -> Self {
+        let degrees: DegreeCoordinate = coordinate.into();
+        degrees.into()
+    }
+}
+
+impl From<RadianCoordinate> for DegreeCoordinate {
+    fn from(other: RadianCoordinate) -> Self {
+        DegreeCoordinate::from_radians(other.lat, other.lon)
+    }
+}
+
+impl RadianCoordinate {
+    /// Builds a RadianCoordinate from latitude and longitude given in
+    /// degrees.
+    pub fn from_degrees(lat: f64, lon: f64) -> Self {
+        let lat = lat / 90.0 * FRAC_PI_2;
+        let lon = normalize_lon(lon / 180.0 * PI);
+
+        RadianCoordinate { lat, lon }
+    }
+
+    /// gives a normalizes version of the Coordinate
+    pub fn normalize(&self) -> RadianCoordinate {
+        RadianCoordinate {
+            lat: self.lat,
+            lon: normalize_lon(self.lon),
+        }
+    }
+
+    /// returns the longitude of a point when the coordinate system is
+    /// transformed, such that `north_pole` is the north pole of that system.
+    pub fn get_transformed_longitude(&self, north_pole: &RadianCoordinate) -> f64 {
+
+        if self.lat == FRAC_PI_2 {
+            return self.lon
+        };
+
+        let top = (self.lon- north_pole.lon).sin() * self.lat.cos();
+        let bottom = (self.lat.sin() * north_pole.lat.cos())
+            - (self.lat.cos() * north_pole.lat.sin() * (self.lon - north_pole.lon).cos());
+
+        bottom.atan2(top)
+    }
+
+    /// returns `true` when the point is on a great circle with point `a` and `b`
+    /// or numerically very close to that
+    pub fn on_great_circle(&self, a: &RadianCoordinate, b: &RadianCoordinate) -> bool {
+        let lat_a = a.get_transformed_longitude(&self);
+        let lat_b = b.get_transformed_longitude(&self);
+
+        (lat_a - lat_b).abs().rem_euclid(PI) < ACCURACY_BOUNDARY
+    }
+
+    /// calculates whether `other` is antipodal to this point.
+    pub fn antipodal(&self, other: &RadianCoordinate) -> bool {
+        // if the distance between both points is very close to PI, they are
+        // antipodal
+
+        let c = self.distance_to(other);
+
+        let diff = (c - PI).abs();
+
+        diff < ACCURACY_BOUNDARY
+    }
+
+    /// returns the shortest distance to an other RadianCoordinate along
+    /// the surface of the unit-sphere in radians.
+    pub fn distance_to(&self, other: &RadianCoordinate) -> f64 {
+        // using the haversine formula
+        // if the distance between both points is very close to PI, they are
+        // antipodal
+
+        let delta_lat = other.lat - self.lat;
+        let delta_lon = other.lon - self.lon;
+
+        let a = (delta_lat / 2.0).sin().powi(2)
+            + self.lat.cos() * other.lat.cos() * (delta_lon / 2.0).sin().powi(2);
+        2.0 * a.sqrt().atan2((1.0_f64 - a).sqrt())
+    }
+
+    /// checks whether the strike is between the shorter angle between the points
+    /// `a` and `b` while using `pole` as the north pole.
+    pub fn strike_between(
+        &self,
+        pole: &RadianCoordinate,
+        a: &RadianCoordinate,
+        b: &RadianCoordinate,
+    ) -> bool {
+        let mut lon = self.get_transformed_longitude(pole).rem_euclid(2.0 * PI);
+        let a_lon = a.get_transformed_longitude(pole).rem_euclid(2.0 * PI);
+        let b_lon = b.get_transformed_longitude(pole).rem_euclid(2.0 * PI);
+
+        // select the start boundary of the smaller circle segment in the
+        // positive direction
+        let start;
+        let mut end;
+
+        if (b_lon - a_lon) > PI {
+            start = b_lon;
+            end = a_lon;
+        } else {
+            start = a_lon;
+            end = b_lon;
+        }
+
+        // rotate the values such that the start is interpreted as 0 and
+        // handle the wrap around.
+        // The start of the interval is now at 0 and the end should be smaller
+        // than PI.
+        end = (end - start).rem_euclid(2.0 * PI);
+        lon = (lon - start).rem_euclid(2.0 * PI);
+
+        0_f64 <= lon && lon <= end
+    }
+}
+
+impl DegreeCoordinate {
+    /// Builds a DegreeCoordinate from latitude and longitude given in
+    /// radians.
+    pub fn from_radians(lat: f64, lon: f64) -> Self {
+        let lat = lat * 90.0 / FRAC_PI_2;
+        let lon = normalize_lon(lon) * 180.0 / PI;
+
+        DegreeCoordinate { lat, lon }
+    }
+}
+
+/// normalizes longitude values given in radians to the range (-PI, PI]
+pub fn normalize_lon(lon: f64) -> f64 {
+    // restrict values to -/+ TAU
+    let mut lon = lon % (TAU);
+
+    if lon <= -PI {
+        lon = TAU + lon;
+    }
+
+    if lon > PI {
+        lon = -TAU + lon;
+    }
+
+    lon
+}
+
+#[test]
+fn test_normalize_lon() {
+    assert!((normalize_lon(0.0) - 0.0).abs() < 10e-12);
+    assert!((normalize_lon(PI) - PI).abs() < 10e-12);
+    assert!((normalize_lon(-PI) - PI).abs() < 10e-12);
+    assert!((normalize_lon(TAU) - 0.0).abs() < 10e-12);
+    assert!((normalize_lon(PI + FRAC_PI_2) - (-FRAC_PI_2)).abs() < 10e-12);
+    assert!((normalize_lon(-PI - FRAC_PI_2) - (FRAC_PI_2)).abs() < 10e-12);
+}