Inputs are latitudes and longitudes in radians. Outputs are meters.
HaversinesPackedMatrix and VincentysPackedMatrix pack a set of points once each.
Their methods score every point of one set, the rows, against every point of another, the columns, into a caller’s row-major slice of rows × columns entries.
use numkong::{HaversinesPackedMatrix, VincentysPackedMatrix};
// Statue of Liberty (40.6892°N, 74.0445°W) → Big Ben (51.5007°N, 0.1246°W)
let liberty_latitudes = [0.7101605100_f64];
let liberty_longitudes = [-1.2923203180_f64];
let big_ben_latitudes = [0.8988567821_f64];
let big_ben_longitudes = [-0.0021746802_f64];
let mut distance = [0.0_f64; 1];
let liberty = VincentysPackedMatrix::new(&liberty_latitudes, &liberty_longitudes).unwrap();
let big_ben = VincentysPackedMatrix::new(&big_ben_latitudes, &big_ben_longitudes).unwrap();
liberty.vincentys(&big_ben, &mut distance).unwrap(); // ≈ 5,589,857 m (ellipsoidal, baseline)
let liberty = HaversinesPackedMatrix::new(&liberty_latitudes, &liberty_longitudes).unwrap();
let big_ben = HaversinesPackedMatrix::new(&big_ben_latitudes, &big_ben_longitudes).unwrap();
liberty.haversines(&big_ben, &mut distance).unwrap(); // ≈ 5,574,848 m (spherical, ~15 km less)Haversine uses the IUGG mean radius of 6,371,008.8 m, and Vincenty the WGS-84 ellipsoid. SqEuclideans3D scores [x, y, z] points the same way, with no packing. i32 coordinates, like Earth-fixed centimeters, give exact u64 squares, and f64 ones give f64.
use numkong::SqEuclideans3D;
let satellites = [[637_813_700_i32, 0, 0], [0, 637_813_700, 0]];
let stations = [[-637_813_700_i32, 0, 0], [0, 0, 637_813_700], [637_813_600, 0, 0]];
let mut squares = [0_u64; 6];
i32::sqeuclideans3d(&satellites, &stations, &mut squares).unwrap();
assert_eq!(squares[2], 10_000); // a meter apartHaversinesPackedMatrixPoints packed once as unit vectors on the mean sphere, scored with great-circle distances as rows against the columns of another pack.6pub struct HaversinesPackedMatrix<Alloc: Allocator = Global>The distance is 2R × atan2(‖a − b‖, ‖a + b‖) over unit vectors a and b, which stays
conditioned at every separation, from a millimeter to the antipode, with R the mean radius.
Debug, Send, Syncnew_inPacks the points at latitudes and longitudes, in radians, using a custom allocator.pub fn new_in(latitudes: &[f64], longitudes: &[f64], alloc: Alloc) -> Result<Self, Error>Returns Err if the two slices differ in length or allocation fails.
*× 2Returns a reference to the allocator.pub fn allocator(&self) -> &Alloc*× 2Number of packed points.pub fn points(&self) -> usizepack_*_size× 2Bytes one pack of points points takes, to size a buffer without packing.pub fn pack__size(points: usize) -> Result<usize, Error>haversinesGreat-circle distances in meters from every point of self to every point of columns, filling results with [rows,columns] entries.pub fn haversines<OtherAlloc: Allocator>(
&self,
columns: &HaversinesPackedMatrix<OtherAlloc>,
results: &mut [f64],
) -> Result<(), Error>Returns Err if results does not hold rows × columns entries.
newPacks the points at latitudes and longitudes, in radians, using the global allocator.pub fn new(latitudes: &[f64], longitudes: &[f64]) -> Result<Self, Error>SqEuclideans3DComputes squared Euclidean distances between every row point and every column point in 3D.2pub trait SqEuclideans3D: Sizedi32 coordinates yield exact u64 squares whenever every axis difference fits 32 bits, which
Earth-fixed centimeters below 4,359 km of altitude always satisfy. f64 coordinates sum their
three squares in Dot2.
impl SqEuclideans3D for i32impl SqEuclideans3D for f64DistanceSquared distance type: u64 for i32 points, f64 for f64 ones.type Distance: Copy + Default;sqeuclideans3dFills results with the [rows,columns] squared distances from every point of a to every point of b.fn sqeuclideans3d(
a: &[[Self; 3]],
b: &[[Self; 3]],
results: &mut [Self::Distance],
) -> Result<(), Error>Returns Err if results does not hold rows × columns entries.
VincentysPackedMatrixPoints packed once as longitudes and the sines and cosines of their reduced latitudes, scored with Vincenty geodesic distances on the WGS-84 ellipsoid as rows against the columns of another pack.3pub struct VincentysPackedMatrix<Alloc: Allocator = Global>Uses Vincenty's iterative formula for oblate spheroid geodesics:
tan(U) = (1−f) × tan(φ)λ → L + (1−C) × f × sin(α) × [σ + C × sin(σ) × ...]u² = cos²(α) × (a² − b²)/b²s = b × A × (σ − Δσ)Where a = equatorial radius, b = polar radius, f = flattening. ~20× more accurate than Haversine for long distances. Inputs are in radians, outputs in meters. A pair still unsettled after the iteration limit, which happens only near the antipode, reports its last iterate.
Debug, Send, Syncnew_inPacks the points at latitudes and longitudes, in radians, using a custom allocator.pub fn new_in(latitudes: &[f64], longitudes: &[f64], alloc: Alloc) -> Result<Self, Error>Returns Err if the two slices differ in length or allocation fails.
vincentysGeodesic distances in meters from every point of self to every point of columns, filling results with [rows,columns] entries.pub fn vincentys<OtherAlloc: Allocator>(
&self,
columns: &VincentysPackedMatrix<OtherAlloc>,
results: &mut [f64],
) -> Result<(), Error>Returns Err if results does not hold rows × columns entries.
newPacks the points at latitudes and longitudes, in radians, using the global allocator.pub fn new(latitudes: &[f64], longitudes: &[f64]) -> Result<Self, Error>