From 4ddeba87a4f0684a056aeb296814120931a455ee Mon Sep 17 00:00:00 2001 From: Cedric Liegeois Date: Mon, 27 Nov 2023 23:19:03 +1100 Subject: [PATCH] spherical caps --- Cargo.toml | 2 +- ChangeLog.md | 4 + README.md | 1 + src/spherical/cap.rs | 379 +++++++++++++++++++++++++++++++++++++ src/spherical/mod.rs | 3 + src/spherical/rectangle.rs | 95 +++++++++- src/spherical/sphere.rs | 76 +++++--- 7 files changed, 525 insertions(+), 35 deletions(-) create mode 100644 src/spherical/cap.rs diff --git a/Cargo.toml b/Cargo.toml index ad3c1f8..0d8a1d2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "jord" -version = "0.9.0" +version = "0.10.0" edition = "2021" rust-version = "1.65" authors = ["Cedric Liegeois ", diff --git a/ChangeLog.md b/ChangeLog.md index 0af0786..0dfe3db 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,7 @@ +### 0.10.0 + +- Spherical caps + ### 0.9.0 - Tests diff --git a/README.md b/README.md index c30b2ff..d20cc72 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,7 @@ The `jord` crate implements various geographical position calculations, featurin - [Great circle](https://en.wikipedia.org/wiki/Great_circle) ([spherical](crate::spherical::Sphere)) navigation: surface distance, initial & final bearing, interpolated position, [minor arc](crate::spherical::MinorArc) intersection, cross track distance, angle turned, side of point..., - Kinematics ([spherical](crate::spherical::Sphere)): closest point of approach between tracks, minimum speed for intercept and time to intercept, - [Spherical Loop](crate::spherical::Loop)s ('simple polygons'): convex/concave, clockwise/anti-clockwise, contains point, [minimum bounding rectangle](crate::spherical::Rectangle), triangulation, spherical excess..., +- [Spherical Cap](crate::spherical::Cap)s and [Rectangular Region](crate::spherical::Rectangle)s - Location-dependent radii of [ellispoid](crate::ellipsoidal::Ellipsoid)s. ## Literature diff --git a/src/spherical/cap.rs b/src/spherical/cap.rs new file mode 100644 index 0000000..e68684c --- /dev/null +++ b/src/spherical/cap.rs @@ -0,0 +1,379 @@ +use std::f64::consts::PI; + +use crate::{Angle, NVector, Vec3}; + +use super::Sphere; + +/// A [spherical cap](https://en.wikipedia.org/wiki/Spherical_cap): a portion of a sphere cut off by a plane. +/// This struct and implementation is very much based on [S2Cap](https://github.com/google/s2geometry/blob/master/src/s2/s2cap.h). +#[derive(PartialEq, Clone, Copy, Debug, Default)] +pub struct Cap { + centre: NVector, + chord_radius2: f64, +} + +impl Cap { + const MAX_CHORD_RADIUS_2: f64 = 4.0; + + /// Empty spherical cap: contains no point. + pub const EMPTY: Cap = Self { + centre: NVector::new(Vec3::UNIT_Z), + chord_radius2: -1.0, + }; + + /// Full spherical cap: contains all points. + pub const FULL: Cap = Self { + centre: NVector::new(Vec3::UNIT_Z), + chord_radius2: Self::MAX_CHORD_RADIUS_2, + }; + + /// Constructs a new cap from the given centre and given radius expressed as the angle between the + /// centre and all points on the boundary of the cap. + pub fn from_centre_and_radius(centre: NVector, radius: Angle) -> Self { + let chord_radius2 = if radius < Angle::ZERO { + -1.0 + } else { + Self::radius_to_chord_radius2(radius) + }; + Self { + centre, + chord_radius2, + } + } + + /// Constructs a new cap from the given centre and a given point on the boundary of the cap. + pub fn from_centre_and_boundary_point(centre: NVector, boundary_point: NVector) -> Self { + Self { + centre, + chord_radius2: Self::chord_radius2(centre, boundary_point), + } + } + + /// Constructs a new cap whose boundary passes by the 3 given points: the returned cap is the circumcircle of the + /// triangle defined by the 3 given points. + pub fn from_points(a: NVector, b: NVector, c: NVector) -> Self { + // see STRIPACK: http://orion.math.iastate.edu/burkardt/f_src/stripack/stripack.f90 + // 3 points must be in anti-clockwise order + let clockwise = Sphere::side(a, b, c) < 0; + let v1 = a.as_vec3(); + let v2 = if clockwise { c.as_vec3() } else { b.as_vec3() }; + let v3 = if clockwise { b.as_vec3() } else { c.as_vec3() }; + let e1 = v2 - v1; + let e2 = v3 - v1; + let centre = NVector::new(e1.orthogonal_to(e2)); + // all chord radius should be equal, still take maximum to account for floating point errors. + let chord_radius2 = Self::chord_radius2(a, centre) + .max(Self::chord_radius2(b, centre).max(Self::chord_radius2(c, centre))); + Self { + centre, + chord_radius2, + } + } + + /// Determines whether this cap is [full](crate::spherical::Cap::FULL). + pub fn is_full(&self) -> bool { + self.chord_radius2 >= Self::MAX_CHORD_RADIUS_2 + } + + /// Determines whether this cap is [empty](crate::spherical::Cap::EMPTY). + pub fn is_empty(&self) -> bool { + self.chord_radius2 < 0.0 + } + + /// Returns the complement of this cap. Both caps have the same boundary but + /// disjoint interiors (the union of both caps is [full](crate::spherical::Cap::FULL)). + pub fn complement(&self) -> Self { + if self.is_empty() { + Self::FULL + } else if self.is_full() { + Self::EMPTY + } else { + Self { + centre: self.centre.antipode(), + chord_radius2: (Self::MAX_CHORD_RADIUS_2 - self.chord_radius2), + } + } + } + + /// Determines whether this cap contains the given point (including the boundary). + /// + /// # Examples + /// + /// ``` + /// use jord::NVector; + /// use jord::spherical::Cap; + /// + /// let cap = Cap::from_centre_and_boundary_point( + /// NVector::from_lat_long_degrees(90.0, 0.0), + /// NVector::from_lat_long_degrees(0.0, 0.0) + /// ); + /// + /// assert!(cap.contains_point(NVector::from_lat_long_degrees(0.0, 0.0))); + /// assert!(cap.contains_point(NVector::from_lat_long_degrees(45.0, 45.0))); + /// ``` + pub fn contains_point(&self, p: NVector) -> bool { + Self::chord_radius2(self.centre, p) <= self.chord_radius2 + } + + /// Determines whether the interior of this cap contains the given point. + /// + /// # Examples + /// + /// ``` + /// use jord::NVector; + /// use jord::spherical::Cap; + /// + /// let cap = Cap::from_centre_and_boundary_point( + /// NVector::from_lat_long_degrees(90.0, 0.0), + /// NVector::from_lat_long_degrees(0.0, 0.0) + /// ); + /// + /// assert!(!cap.interior_contains_point(NVector::from_lat_long_degrees(0.0, 0.0))); + /// assert!(cap.interior_contains_point(NVector::from_lat_long_degrees(45.0, 45.0))); + /// ``` + pub fn interior_contains_point(&self, p: NVector) -> bool { + Self::chord_radius2(self.centre, p) < self.chord_radius2 + } + + /// Determines whether this cap contains the given cap. The full cap contains all caps + /// and the empty cap is contained by all caps. + /// + /// # Examples + /// + /// ``` + /// use jord::NVector; + /// use jord::spherical::Cap; + /// + /// let cap1 = Cap::from_centre_and_boundary_point( + /// NVector::from_lat_long_degrees(90.0, 0.0), + /// NVector::from_lat_long_degrees(0.0, 0.0) + /// ); + /// + /// let cap2 = Cap::from_centre_and_boundary_point( + /// NVector::from_lat_long_degrees(90.0, 0.0), + /// NVector::from_lat_long_degrees(45.0, 0.0) + /// ); + /// + /// assert!(cap1.contains_cap(cap2)); + /// ``` + pub fn contains_cap(&self, other: Self) -> bool { + if self.is_full() || other.is_empty() { + true + } else { + self.chord_radius2 + >= Self::chord_radius2(self.centre, other.centre) + other.chord_radius2 + } + } + + ///Returns the smallest cap which encloses this cap and the other given cap. + pub fn union(&self, other: Self) -> Self { + if self.chord_radius2 < other.chord_radius2 { + return other.union(*self); + } + if self.is_full() || other.is_empty() { + return *self; + } + + let self_radius = self.radius(); + let other_radius = other.radius(); + let distance = Sphere::angle(self.centre, other.centre); + if self_radius >= distance + other_radius { + return *self; + } + let union_radius = 0.5 * (distance + self_radius + other_radius); + let ang = 0.5 * (distance - self_radius + other_radius); + let centre = Sphere::position_on_great_circle(self.centre, other.centre, ang); + Self { + centre, + chord_radius2: Self::radius_to_chord_radius2(union_radius), + } + } + + /// Returns the centre of this cap. + pub fn centre(&self) -> NVector { + self.centre + } + + /// Returns the radius of this cap: central angle between the centre of this cap and + /// any point on the boundary (negative for [empty](crate::spherical::Cap::EMPTY) caps). + /// The returned value may not exactly equal the value passed + /// to [from_centre_and_boundary_point](crate::spherical::Cap::from_centre_and_boundary_point). + /// + /// # Examples + /// + /// ``` + /// use std::f64::consts::PI; + /// + /// use jord::{Angle, NVector}; + /// use jord::spherical::Cap; + /// + /// let cap = Cap::from_centre_and_boundary_point( + /// NVector::from_lat_long_degrees(90.0, 0.0), + /// NVector::from_lat_long_degrees(45.0, 45.0) + /// ); + /// + /// assert_eq!(Angle::from_radians(PI / 4.0), cap.radius().round_d7()); + /// ``` + pub fn radius(&self) -> Angle { + if self.is_empty() { + Angle::from_radians(-1.0) + } else { + Angle::from_radians(2.0 * (self.chord_radius2.sqrt() * 0.5).asin()) + } + } + + fn chord_radius2(a: NVector, b: NVector) -> f64 { + (a.as_vec3() - b.as_vec3()) + .squared_norm() + .min(Self::MAX_CHORD_RADIUS_2) + } + + fn radius_to_chord_radius2(radius: Angle) -> f64 { + // max angle is PI + let chord_radius = 2.0 * ((radius.as_radians().min(PI)) * 0.5).sin(); + (chord_radius * chord_radius).min(Self::MAX_CHORD_RADIUS_2) + } +} + +#[cfg(test)] +mod tests { + use crate::{positions::assert_nv_eq_d7, spherical::Cap, Angle, NVector}; + use std::f64::consts::PI; + + #[test] + fn full() { + assert!(Cap::FULL.contains_point(NVector::from_lat_long_degrees(90.0, 0.0))); + assert!(Cap::FULL.contains_point(NVector::from_lat_long_degrees(-90.0, 0.0))); + assert_eq!(Angle::from_radians(PI), Cap::FULL.radius()); + assert_eq!(Cap::EMPTY, Cap::FULL.complement()); + } + + #[test] + fn empty() { + assert!(!Cap::EMPTY.contains_point(NVector::from_lat_long_degrees(90.0, 0.0))); + assert!(!Cap::EMPTY.contains_point(NVector::from_lat_long_degrees(-90.0, 0.0))); + assert_eq!(Angle::from_radians(-1.0), Cap::EMPTY.radius()); + assert_eq!(Cap::FULL, Cap::EMPTY.complement()); + } + + #[test] + fn from_points() { + let a = NVector::from_lat_long_degrees(0.0, 0.0); + let b = NVector::from_lat_long_degrees(20.0, 0.0); + let c = NVector::from_lat_long_degrees(10.0, 10.0); + let cap = Cap::from_points(a, b, c); + assert!(cap.contains_point(a)); + assert!(cap.contains_point(b)); + assert!(cap.contains_point(c)); + + let o = Cap::from_points(c, b, a); + assert_nv_eq_d7(o.centre, cap.centre); + assert!((o.chord_radius2 - cap.chord_radius2).abs() < 1e-16); + } + + #[test] + fn complement() { + let np = NVector::from_lat_long_degrees(90.0, 0.0); + let sp = NVector::from_lat_long_degrees(-90.0, 0.0); + let northern = Cap::from_centre_and_radius(sp, Angle::QUARTER_CIRCLE); + let southern = Cap::from_centre_and_radius(np, Angle::QUARTER_CIRCLE); + + let northern_complement = northern.complement(); + assert_eq!(southern.centre, northern_complement.centre); + assert!((southern.chord_radius2 - northern_complement.chord_radius2).abs() < 1e15); + + let southern_complement = southern.complement(); + assert_eq!(northern.centre, southern_complement.centre); + assert!((northern.chord_radius2 - southern_complement.chord_radius2).abs() < 1e15); + } + + #[test] + fn contains_point() { + let cap = Cap::from_centre_and_boundary_point( + NVector::from_lat_long_degrees(90.0, 0.0), + NVector::from_lat_long_degrees(0.0, 0.0), + ); + assert!(cap.contains_point(NVector::from_lat_long_degrees(0.0, 0.0))); + assert!(cap.contains_point(NVector::from_lat_long_degrees(45.0, 45.0))); + } + + #[test] + fn interior_contains_point() { + let cap = Cap::from_centre_and_boundary_point( + NVector::from_lat_long_degrees(90.0, 0.0), + NVector::from_lat_long_degrees(0.0, 0.0), + ); + assert!(!cap.interior_contains_point(NVector::from_lat_long_degrees(0.0, 0.0))); + assert!(cap.interior_contains_point(NVector::from_lat_long_degrees(45.0, 45.0))); + } + + #[test] + fn contains_cap() { + let c = Cap::from_centre_and_radius( + NVector::from_lat_long_degrees(30.0, 30.0), + Angle::from_degrees(10.0), + ); + assert!(Cap::FULL.contains_cap(c)); + assert!(c.contains_cap(Cap::EMPTY)); + + let o = Cap::from_centre_and_radius( + NVector::from_lat_long_degrees(30.0, 30.0), + Angle::from_degrees(20.0), + ); + assert!(!c.contains_cap(o)); + assert!(o.contains_cap(c)); + } + + #[test] + fn radius() { + assert_eq!( + Angle::QUARTER_CIRCLE, + Cap::from_centre_and_boundary_point( + NVector::from_lat_long_degrees(90.0, 0.0), + NVector::from_lat_long_degrees(0.0, 0.0) + ) + .radius() + .round_d7() + ); + assert_eq!( + Angle::from_radians(PI / 4.0), + Cap::from_centre_and_boundary_point( + NVector::from_lat_long_degrees(90.0, 0.0), + NVector::from_lat_long_degrees(45.0, 45.0) + ) + .radius() + .round_d7() + ); + } + + #[test] + fn union() { + assert!(Cap::FULL.union(Cap::EMPTY).is_full()); + assert!(Cap::EMPTY.union(Cap::FULL).is_full()); + + // a and b have same centre, but radius of a < radius of b. + let a = Cap::from_centre_and_radius( + NVector::from_lat_long_degrees(50.0, 10.0), + Angle::from_radians(0.2), + ); + let b = Cap::from_centre_and_radius( + NVector::from_lat_long_degrees(50.0, 10.0), + Angle::from_radians(0.3), + ); + + assert!(b.contains_cap(a)); + assert_eq!(b, a.union(b)); + + assert_eq!(Cap::FULL, a.union(Cap::FULL)); + assert_eq!(a, a.union(Cap::EMPTY)); + + // a and c have different centers, one entirely encompasses the other. + let c = Cap::from_centre_and_radius( + NVector::from_lat_long_degrees(51.0, 11.0), + Angle::from_radians(1.5), + ); + assert!(c.contains_cap(a)); + assert_eq!(a.union(c).centre(), c.centre()); + assert_eq!(a.union(c).radius(), c.radius()); + } +} diff --git a/src/spherical/mod.rs b/src/spherical/mod.rs index 0b48387..5bad0dc 100644 --- a/src/spherical/mod.rs +++ b/src/spherical/mod.rs @@ -2,6 +2,9 @@ mod base; +mod cap; +pub use cap::Cap; + mod great_circle; pub use great_circle::GreatCircle; diff --git a/src/spherical/rectangle.rs b/src/spherical/rectangle.rs index e527411..2872a56 100644 --- a/src/spherical/rectangle.rs +++ b/src/spherical/rectangle.rs @@ -16,15 +16,15 @@ pub struct Rectangle { lng: LongitudeInterval, } -// TODO(CL): Exmaples +// TODO(CL): Exmaples + (Interior)Intersection impl Rectangle { - /// Empty rectangle. + /// Empty rectangle: contains no point. pub const EMPTY: Rectangle = Self { lat: LatitudeInterval::EMPTY, lng: LongitudeInterval::EMPTY, }; - /// Creates a full rectangle. + /// Full rectangle: contains all points. pub const FULL: Rectangle = Self { lat: LatitudeInterval::FULL, lng: LongitudeInterval::FULL, @@ -97,7 +97,7 @@ impl Rectangle { } } - /// Determines whether this rectangle contains the given point. + /// Determines whether this rectangle contains the given point (boundaries included). /// /// # Examples /// @@ -114,6 +114,9 @@ impl Rectangle { /// /// assert!(a.contains_point(LatLong::from_degrees(10.0, 10.0))); /// + /// /// point on boundary + /// assert!(a.contains_point(LatLong::from_degrees(30.0, 00.0))); + /// /// // latitude above north. /// assert!(!a.contains_point(LatLong::from_degrees(40.0, 10.0))); /// @@ -130,11 +133,52 @@ impl Rectangle { self.lat.contains_lat(p.latitude()) && self.lng.contains_lng(p.longitude()) } + /// Determines whether the interior of this rectangle contains the given point (i.e. boundaries excluded). + /// # Examples + /// + /// ``` + /// use jord::{Angle, LatLong}; + /// use jord::spherical::Rectangle; + /// + /// let a = Rectangle::from_nesw( + /// Angle::from_degrees(30.0), + /// Angle::from_degrees(30.0), + /// Angle::ZERO, + /// Angle::ZERO + /// ); + /// + /// assert!(a.interior_contains_point(LatLong::from_degrees(10.0, 10.0))); + /// + /// /// point on boundary + /// assert!(!a.interior_contains_point(LatLong::from_degrees(30.0, 00.0))); + /// + /// // latitude above north. + /// assert!(!a.interior_contains_point(LatLong::from_degrees(40.0, 10.0))); + /// + /// // latitude below south. + /// assert!(!a.interior_contains_point(LatLong::from_degrees(-1.0, 10.0))); + /// + /// // longitude after east. + /// assert!(!a.interior_contains_point(LatLong::from_degrees(10.0, 40.0))); + /// + /// // longitude after west. + /// assert!(!a.interior_contains_point(LatLong::from_degrees(10.0, -10.0))); + /// ``` + pub fn interior_contains_point(&self, p: LatLong) -> bool { + self.lat.interior_contains_lat(p.latitude()) + && self.lng.interior_contains_lng(p.longitude()) + } + /// Determines whether this rectangle contains the given rectangle. pub fn contains_rectangle(&self, r: Rectangle) -> bool { self.lat.contains_int(r.lat) && self.lng.contains_int(r.lng) } + /// Determines whether the interior of this rectangle contains the given rectangle (including its boundary). + pub fn interior_contains_rectangle(&self, r: Rectangle) -> bool { + self.lat.interior_contains_int(r.lat) && self.lng.interior_contains_int(r.lng) + } + /// Determines whether this rectangle is [full](crate::spherical::Rectangle::FULL). pub fn is_full(&self) -> bool { self.is_latitude_full() && self.is_longitude_full() @@ -321,11 +365,16 @@ impl LatitudeInterval { Self::new(lo, hi) } - /// Returns true if and only if this latitude interval contains the given latitude. + /// Returns true if and only if this latitude interval (closed) contains the given latitude. fn contains_lat(&self, latitude: Angle) -> bool { latitude >= self.lo && latitude <= self.hi } + /// Returns true if and only if this latitude interval (opened) contains the given latitude. + fn interior_contains_lat(&self, latitude: Angle) -> bool { + latitude > self.lo && latitude < self.hi + } + /// Returns true if and only if this latitude interval contains the given latitude interval. fn contains_int(&self, o: Self) -> bool { if o.is_empty() { @@ -335,6 +384,15 @@ impl LatitudeInterval { } } + /// Returns true if and only if the interior of this latitude interval contains the given latitude interval. + fn interior_contains_int(&self, o: Self) -> bool { + if o.is_empty() { + true + } else { + o.lo > self.lo && o.hi < self.hi + } + } + /// Returns an interval that has been expanded/shrinked on each side by the given amount. fn expand(&self, amount: Angle) -> Self { if self.is_empty() { @@ -474,8 +532,9 @@ impl LongitudeInterval { } } + /// Returns true if and only if this longitude interval (closed) contains the given longitude. fn contains_lng(&self, longitude: Angle) -> bool { - let lng = Self::normalised_longitude(longitude); + let lng: Angle = Self::normalised_longitude(longitude); if self.is_inverted() { (lng >= self.lo || lng <= self.hi) && !self.is_empty() } else { @@ -483,6 +542,16 @@ impl LongitudeInterval { } } + /// Returns true if and only if this longitude interval (opened) contains the given longitude. + fn interior_contains_lng(&self, longitude: Angle) -> bool { + let lng: Angle = Self::normalised_longitude(longitude); + if self.is_inverted() { + lng > self.lo || lng < self.hi + } else { + (lng > self.lo && lng < self.hi) || self.is_full() + } + } + /// Returns true if and only if this longitude interval contains the given longitude interval. fn contains_int(&self, o: Self) -> bool { if self.is_inverted() { @@ -497,6 +566,20 @@ impl LongitudeInterval { o.lo >= self.lo && o.hi <= self.hi } + /// Returns true if and only if the interior of this longitude interval contains the given longitude interval. + fn interior_contains_int(&self, o: Self) -> bool { + if self.is_inverted() { + if !o.is_inverted() { + return o.lo > self.lo || o.hi < self.hi; + } + return (o.lo > self.lo && o.hi < self.hi) || o.is_empty(); + } + if o.is_inverted() { + return self.is_full() || o.is_empty(); + } + (o.lo > self.lo && o.hi < self.hi) || self.is_full() + } + /// Returns true if this longitude interval is full. fn is_full(&self) -> bool { self.lo == Angle::NEG_HALF_CIRCLE && self.hi == Angle::HALF_CIRCLE diff --git a/src/spherical/sphere.rs b/src/spherical/sphere.rs index 5f0c150..d51a762 100644 --- a/src/spherical/sphere.rs +++ b/src/spherical/sphere.rs @@ -251,24 +251,41 @@ impl Sphere { } } + /// Returns the position at the given distance `a` from given point `p1` along the great circle `p1, p2`. + /// + /// ``` + /// use std::f64::consts::PI; + /// use jord::{Angle, NVector}; + /// use jord::spherical::Sphere; + /// + /// let r = Sphere::position_on_great_circle( + /// NVector::from_lat_long_degrees(0.0, 0.0), + /// NVector::from_lat_long_degrees(0.0, 1.0), + /// Angle::from_radians(PI / 4.0) + /// ); + /// + /// assert_eq!(NVector::from_lat_long_degrees(0.0, 45.0), r); + /// ``` + pub fn position_on_great_circle(p1: NVector, p2: NVector, a: Angle) -> NVector { + let rads = a.as_radians(); + // direction from v0 to v1. + let dir = (p1.as_vec3().orthogonal_to(p2.as_vec3())).cross_prod_unit(p1.as_vec3()); + let v = (p1.as_vec3() * rads.cos() + dir * rads.sin()).unit(); + NVector::new(v) + } + /// Computes the position at given fraction between this position and the given position. - /// Returns `None` if: - /// - the fraction is `< 0` or `> 1`, - /// - this position and the given position are the antipodes of one another. + /// Returns `None` if the given fraction is `< 0` or `> 1`.` pub fn interpolated_pos(p1: NVector, p2: NVector, f: f64) -> Option { - if !(0.0..=1.0).contains(&f) || p1.is_antipode_of(p2) { + if !(0.0..=1.0).contains(&f) { None } else if f == 0.0 { Some(p1) } else if f == 1.0 { Some(p2) } else { - // angular distance in radians multiplied by the fraction: how far from v0. - let distance_radians = f * angle_radians_between(p1.as_vec3(), p2.as_vec3(), None); - // a vector representing the direction from v0 to v1. - let dir = (p1.as_vec3().stable_cross_prod(p2.as_vec3())).cross_prod_unit(p1.as_vec3()); - let v = (p1.as_vec3() * distance_radians.cos() + dir * distance_radians.sin()).unit(); - Some(NVector::new(v)) + let a = f * Self::angle(p1, p2); + Some(Self::position_on_great_circle(p1, p2, a)) } } @@ -866,7 +883,7 @@ mod tests { // destination. #[test] - fn destination_across_date_line() { + fn destination_pos_across_date_line() { let p = NVector::from_lat_long_degrees(0.0, 154.0); let actual = Sphere::EARTH.destination_pos( p, @@ -878,7 +895,7 @@ mod tests { } #[test] - fn destination_from_north_pole() { + fn destination_pos_from_north_pole() { let expected = NVector::from_lat_long_degrees(45.0, 0.0); let distance = Sphere::EARTH.radius() * (PI / 4.0); let actual = Sphere::EARTH.destination_pos( @@ -890,7 +907,7 @@ mod tests { } #[test] - fn destination_from_south_pole() { + fn destination_pos_from_south_pole() { let expected = NVector::from_lat_long_degrees(-45.0, 0.0); let distance = Sphere::EARTH.radius() * (PI / 4.0); let actual = Sphere::EARTH.destination_pos( @@ -902,7 +919,7 @@ mod tests { } #[test] - fn destination_negative_distance() { + fn destination_pos_negative_distance() { let p = NVector::from_lat_long_degrees(0.0, 0.0); // equivalent of -10 degree of longitude. let d = Sphere::EARTH.radius() * (-2.0 * PI / 36.0); @@ -912,7 +929,7 @@ mod tests { } #[test] - fn destination_travelled_longitude_greater_than_90() { + fn destination_pos_travelled_longitude_greater_than_90() { let p = NVector::from_lat_long_degrees(60.2, 11.1); let d = Sphere::EARTH.destination_pos( p, @@ -924,7 +941,7 @@ mod tests { } #[test] - fn destination_zero_distance() { + fn destination_pos_zero_distance() { let p = NVector::from_lat_long_degrees(55.6050, 13.0038); assert_eq!( p, @@ -1153,27 +1170,30 @@ mod tests { // interpolated #[test] - fn interpolated_antipodal() { + fn interpolated_pos_antipodal() { let p = NVector::from_lat_long_degrees(90.0, 0.0); - assert!(Sphere::interpolated_pos(p, p.antipode(), 0.5).is_none()); + assert_opt_nv_eq_d7( + NVector::from_lat_long_degrees(0.0, 90.0), + Sphere::interpolated_pos(p, p.antipode(), 0.5), + ); } #[test] - fn interpolated_f0() { + fn interpolated_pos_f0() { let p1 = NVector::from_lat_long_degrees(90.0, 0.0); let p2 = NVector::from_lat_long_degrees(0.0, 0.0); assert_eq!(Some(p1), Sphere::interpolated_pos(p1, p2, 0.0)); } #[test] - fn interpolated_f1() { + fn interpolated_pos_f1() { let p1 = NVector::from_lat_long_degrees(90.0, 0.0); let p2 = NVector::from_lat_long_degrees(0.0, 0.0); assert_eq!(Some(p2), Sphere::interpolated_pos(p1, p2, 1.0)); } #[test] - fn interpolated_invalid_f() { + fn interpolated_pos_invalid_f() { let p1 = NVector::from_lat_long_degrees(0.0, 0.0); let p2 = NVector::from_lat_long_degrees(1.0, 0.0); assert!(Sphere::interpolated_pos(p1, p2, -0.1).is_none()); @@ -1181,19 +1201,19 @@ mod tests { } #[test] - fn interpolated_half() { - assert_eq!( - Some(NVector::from_lat_long_degrees(0.0, 0.0)), + fn interpolated_pos_half() { + assert_opt_nv_eq_d7( + NVector::from_lat_long_degrees(0.0, 0.0), Sphere::interpolated_pos( NVector::from_lat_long_degrees(10.0, 0.0), NVector::from_lat_long_degrees(-10.0, 0.0), - 0.5 - ) + 0.5, + ), ); } #[test] - fn interpolated_side() { + fn interpolated_pos_side() { let p0 = NVector::from_lat_long_degrees(154.0, 54.0); let p1 = NVector::from_lat_long_degrees(155.0, 55.0); let i = Sphere::interpolated_pos(p0, p1, 0.25).unwrap(); @@ -1201,7 +1221,7 @@ mod tests { } #[test] - fn interpolated_transitivity() { + fn interpolated_pos_transitivity() { let p0 = NVector::from_lat_long_degrees(10.0, 0.0); let p1 = NVector::from_lat_long_degrees(-10.0, 0.0); let expected = Sphere::interpolated_pos(p0, p1, 0.5).unwrap();