Handle floating point errors in Haversine distance

Author: juliohm

src/haversine.jl

@@ -12,7 +12,7 @@ end
 
 const VecOrLengthTwoTuple{T} = Union{AbstractVector{T}, NTuple{2, T}}
 
-function evaluate(dist::Haversine, x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple) 
+function evaluate(dist::Haversine, x::VecOrLengthTwoTuple{T}, y::VecOrLengthTwoTuple{T}) where {T}
     length(x) == length(y) == 2 || haversine_error()
 
     @inbounds begin
@@ -28,12 +28,14 @@ function evaluate(dist::Haversine, x::VecOrLengthTwoTuple, y::VecOrLengthTwoTupl
 
     # haversine formula
     a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2
-    c = 2atan2(√a, √(1-a))
+
+    # take care of floating point errors
+    a = min(a, one(T))
 
     # distance on the sphere
-    c*dist.radius
+    2*dist.radius*atan2(√a, √(1-a))
 end
 
 haversine(x::VecOrLengthTwoTuple, y::VecOrLengthTwoTuple, radius::Real) = evaluate(Haversine(radius), x, y)
 
-@noinline haversine_error() = throw(ArgumentError("expected both inputs to have length 2 in Haversine distance"))
+@noinline haversine_error() = throw(ArgumentError("expected both inputs to have length 2 in Haversine distance"))

test/test_dists.jl

@@ -284,6 +284,7 @@ end #testset
         @test haversine([0.,-90.],  [0.,90.],  1.) ≈ π atol=1e-10
         @test haversine((-180.,0.), (180.,0.), 1.) ≈ 0 atol=1e-10
         @test haversine((0.,-90.),  (0.,90.),  1.) ≈ π atol=1e-10
+        @test haversine((1.,-15.625), (-179.,15.625), 6371.) ≈ 20015. atol=1e0
         @test_throws ArgumentError haversine([0.,-90., 0.25], [0.,90.], 1.)
     end
 end