Struct Ellipsoid

Source

pub struct Ellipsoid { /* private fields */ }

Expand description

A reference ellipsoid of revolution, oblate or spherical.

Implementations§

Source §

impl Ellipsoid

Source

pub fn geodesic_inverse( &self, from: Geodetic, to: Geodetic, ) -> Result<GeodesicInverse, CoreError>

The geodesic from from to to (Karney 2013 §4, §5): its length and the azimuths at both ends. Heights are ignored.

Two coincident points give a zero distance. Where the points are nearly antipodal, or one is near a pole, the azimuths are ill-conditioned: a tiny move of a point turns them through large angles while the distance hardly changes. And the shortest path is not always unique: when φ₂ = −φ₁ exactly, two geodesics of the same length join the points (unless α₁ = α₂), the second with α₁ and α₂ swapped (GeographicLib’s GeodSolve manual, Multiple solutions); either may be returned.

use hpr_core::geodesy::{Ellipsoid, Geodetic};

// Two points a degree of longitude apart on the equator.
let a = Geodetic::from_degrees(0.0, 0.0, 0.0)?;
let b = Geodetic::from_degrees(0.0, 1.0, 0.0)?;
let g = Ellipsoid::WGS84.geodesic_inverse(a, b)?;
// An arc of the equator: a·(π/180).
assert!((g.distance_m - 111_319.490_793_273_57).abs() < 1e-8);
assert!((g.initial_azimuth_rad.to_degrees() - 90.0).abs() < 1e-12);

§Errors

CoreError::Domain if the ellipsoid’s flattening is past GEODESIC_MAX_FLATTENING, or either point fails Geodetic::validated’s checks.

Source

pub fn geodesic_direct( &self, from: Geodetic, azimuth_rad: f64, distance_m: f64, ) -> Result<GeodesicDirect, CoreError>

The point distance_m along the geodesic leaving from at azimuth azimuth_rad (clockwise from north), with the azimuth there (Karney 2013 §2). from’s height is ignored. A negative distance runs backwards along the same geodesic. Karney’s accuracy is shown up to half a meridian (20,004 km); a distance of many circuits also carries its own rounding, one ulp of distance_m (at least 15 nm past 2²⁶ m, 67,109 km).

use hpr_core::geodesy::{Ellipsoid, Geodetic};

// A degree of longitude east along the equator, as above.
let a = Geodetic::from_degrees(0.0, 0.0, 0.0)?;
let d = Ellipsoid::WGS84.geodesic_direct(a, 90f64.to_radians(), 111_319.490_793_273_57)?;
assert!((d.longitude_rad.to_degrees() - 1.0).abs() < 1e-12);
assert!(d.latitude_rad.abs() < 1e-15);

§Errors

CoreError::Domain if the ellipsoid’s flattening is past GEODESIC_MAX_FLATTENING, from fails Geodetic::validated’s checks, or the azimuth or the distance is not finite.

Source §

impl Ellipsoid

Source

pub const WGS84: Self

The WGS 84 ellipsoid: a = 6378137.0 m, 1/f = 298.257223563 (NGA.STND.0036_1.0.0_WGS84, Table 3.1).

Source

pub fn new( semi_major_axis_m: f64, inverse_flattening: f64, ) -> Result<Self, CoreError>

An ellipsoid from its semi-major axis a and inverse flattening 1/f. An infinite 1/f is a sphere.

§Errors

CoreError::Domain unless a is finite and positive and 1/f is greater than one (or +∞).

Source

pub fn from_flattening( semi_major_axis_m: f64, flattening: f64, ) -> Result<Self, CoreError>

An ellipsoid from its semi-major axis a and flattening f (zero for a sphere).

§Errors

CoreError::Domain unless a is finite and positive and 0 ≤ f < 1.

Source

pub fn semi_major_axis_m(&self) -> f64

Semi-major (equatorial) axis a, m.

Source

pub fn flattening(&self) -> f64

Flattening f = (a − b)/a.

Source

pub fn inverse_flattening(&self) -> f64

Inverse flattening 1/f (+∞ for a sphere).

Source

pub fn semi_minor_axis_m(&self) -> f64

Semi-minor (polar) axis b = a(1 − f), m.

Source

pub fn eccentricity_squared(&self) -> f64

First eccentricity squared e² = f(2 − f) = (a² − b²)/a².

Source

pub fn linear_eccentricity_m(&self) -> f64

Linear eccentricity E = √(a² − b²) = a e, m.

Source

pub fn prime_vertical_radius_m(&self, latitude_rad: f64) -> f64

Radius of curvature in the prime vertical, N(φ) = a / √(1 − e² sin²φ), m (eq. 4-15).

Source

pub fn ecef_from_geodetic(&self, point: Geodetic) -> DVec3

ECEF position of a geodetic point (NGA.STND.0036 eq. 4-14):

X = (N + h) cos φ cos λ,   Y = (N + h) cos φ sin λ,   Z = ((b²/a²) N + h) sin φ

Source

pub fn geodetic_from_ecef( &self, position_ecef_m: DVec3, ) -> Result<Geodetic, CoreError>

Geodetic position of an ECEF point (Karney 2011, appendix B). With R = √(X² + Y²), x = R/a and y = √(1 − e²) Z/a, the largest real root κ of

κ⁴ + 2e²κ³ − (x² + y² − e⁴)κ² − 2e²y²κ − e⁴y² = 0                        (B1)

gives φ = ph(R/(κ + e²) + iZ/κ) (B2) and h = (1 − (1 − e²)/κ) √(D² + Z²) with D = κR/(κ + e²) (B3). κ comes from the resolvent cubic and the factorization (B4)–(B5) in the round-off-avoiding form the paper gives. λ = ph(X + iY), which is 0 on the axis.

§Errors

CoreError::Domain if the position is not finite, or lies in the degenerate set where the closed form needs limiting forms: the equatorial plane within a e² of the center (about 42.7 km for WGS 84), deep inside the Earth.