Add function to compute angular separation between coordinates

JuliaAstro · Nov 7, 2016 · 2cc5e67 · 2cc5e67
1 parent 10053c5
commit 2cc5e67
Show file tree

Hide file tree

Showing 2 changed files with 62 additions and 52 deletions.
diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
@@ -4,7 +4,8 @@ using Compat
 export AbstractSkyCoords,
        ICRSCoords,
        GalCoords,
-       FK5Coords
+       FK5Coords,
+       separation
 
 import Base: convert, *, transpose
 
@@ -195,10 +196,13 @@ end
 # -----------------------------------------------------------------------------
 # Type-dependent methods
 
+lon(c::GalCoords) = c.l
+lat(c::GalCoords) = c.b
+lon(c::AbstractSkyCoords) = c.ra
+lat(c::AbstractSkyCoords) = c.dec
+
 # Abstract away specific field names (ra, dec vs l, b)
-coords2cart(c::ICRSCoords) = coords2cart(c.ra, c.dec)
-coords2cart(c::GalCoords) = coords2cart(c.l, c.b)
-coords2cart(c::FK5Coords) = coords2cart(c.ra, c.dec)
+coords2cart(c::AbstractSkyCoords) = coords2cart(lon(c), lat(c))
 
 # Rotation matrix between coordinate systems: `rotmat(to, from)`
 # Note that all of these return Matrix33{Float64}, regardless of
@@ -231,13 +235,10 @@ rotmat{e1, T1, e2, T2}(::Type{FK5Coords{e1,T1}}, ::Type{FK5Coords{e2,T2}}) =
     precess_from_j2000(e1) * precess_from_j2000(e2)'
 
 # get floating point type in coordinates
-_eltype{e,T}(::FK5Coords{e,T}) = T
-_eltype{T}(::GalCoords{T}) = T
-_eltype{T}(::ICRSCoords{T}) = T
 _eltype{e,T}(::Type{FK5Coords{e,T}}) = T
 _eltype{T}(::Type{GalCoords{T}}) = T
 _eltype{T}(::Type{ICRSCoords{T}}) = T
-
+_eltype(c::AbstractSkyCoords) = _eltype(typeof(c))
 
 # Scalar coordinate conversions
 convert{T<:AbstractSkyCoords}(::Type{T}, c::T) = c
@@ -265,4 +266,47 @@ function convert{T<:AbstractSkyCoords, n, S<:AbstractSkyCoords}(
     result
 end
 
+# ------------------------------------------------------------------------------
+# Distance between coordinates
+
+function _separation(λ_1, ϕ_1, λ_2, ϕ_2)
+    Δλ = λ_2 - λ_1
+    sin_Δλ = sin(Δλ)
+    cos_Δλ = cos(Δλ)
+    sin_ϕ1 = sin(ϕ_1)
+    sin_ϕ2 = sin(ϕ_2)
+    cos_ϕ1 = cos(ϕ_1)
+    cos_ϕ2 = cos(ϕ_2)
+    return atan2(hypot(cos_ϕ2 * sin_Δλ,
+                       cos_ϕ1 * sin_ϕ2 - sin_ϕ1 * cos_ϕ2 * cos_Δλ),
+                 sin_ϕ1 * sin_ϕ2 + cos_ϕ1 * cos_ϕ2 * cos_Δλ)
+end
+
+"""
+    separation(c1::AbstractSkyCoords, c2::AbstractSkyCoords) -> distance
+
+Return angular separation between two sky coordinates, in radians.
+
+The angular separation is calculated using the Vincenty formula
+(http://en.wikipedia.org/wiki/Great-circle_distance), which is slightly more
+complex and computationally expensive than some alternatives, but is stable at
+at all distances, including the poles and antipodes.
+"""
+separation{T<:AbstractSkyCoords}(c1::T, c2::T) =
+    _separation(lon(c1), lat(c1), lon(c2), lat(c2))
+
+separation{T1<:AbstractSkyCoords,T2<:AbstractSkyCoords}(c1::T1, c2::T2) =
+    separation(c1, convert(T1, c2))
+
+function separation{T<:AbstractSkyCoords}(c1::AbstractArray{T},
+                                          c2::AbstractArray{T})
+    @assert size(c1) == size(c2) "Size mismatch"
+    result = similar(c1, _eltype(first(c1)))
+    for i in eachindex(c1)
+        result[i] = separation(c1[i], c2[i])
+    end
+    return result
+end
+
+
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -5,49 +5,9 @@
 using SkyCoords
 using Base.Test
 
-datapath = joinpath(dirname(@__FILE__), "data")
-
-# Angular separation between two points (angles in radians)
-#
-#   The angular separation is calculated using the Vincenty formula [1]_,
-#   which is slighly more complex and computationally expensive than
-#   some alternatives, but is stable at at all distances, including the
-#   poles and antipodes.
-#
-#   [1] http://en.wikipedia.org/wiki/Great-circle_distance
-function angsep(lon1, lat1, lon2, lat2)
-
-    sdlon = sin(lon2 - lon1)
-    cdlon = cos(lon2 - lon1)
-    slat1 = sin(lat1)
-    slat2 = sin(lat2)
-    clat1 = cos(lat1)
-    clat2 = cos(lat2)
-
-    num1 = clat2 * sdlon
-    num2 = clat1 * slat2 - slat1 * clat2 * cdlon
-    denom = slat1 * slat2 + clat1 * clat2 * cdlon
-
-    atan2(sqrt(num1*num1 + num2*num2), denom)
-end
+import SkyCoords: lat, lon
 
-lon(c::GalCoords) = c.l
-lat(c::GalCoords) = c.b
-lon(c::AbstractSkyCoords) = c.ra
-lat(c::AbstractSkyCoords) = c.dec
-
-function angsep{T<:AbstractSkyCoords}(c1::T, c2::T)
-    angsep(lon(c1), lat(c1), lon(c2), lat(c2))
-end
-
-function angsep{T<:AbstractSkyCoords}(c1::Array{T}, c2::Array{T})
-    size(c1) == size(c2) || error("size mismatch")
-    result = similar(c1, SkyCoords._eltype(first(c1)))
-    for i=1:length(c1)
-        result[i] = angsep(c1[i], c2[i])
-    end
-    result
-end
+datapath = joinpath(dirname(@__FILE__), "data")
 
 rad2arcsec(r) = 3600.*rad2deg(r)
 
@@ -77,7 +37,7 @@ for (F, TOL) in ((Float32, 0.2), (Float64, 0.0001), (BigFloat, 0.0001))
             c_ref = S[S(refdata[i, 1], refdata[i, 2]) for i=1:size(refdata,1)]
 
             # compare
-            sep = angsep(c_out, c_ref)
+            sep = separation(c_out, c_ref)
             maxsep = rad2arcsec(maximum(sep))
             meansep = rad2arcsec(mean(sep))
             minsep = rad2arcsec(minimum(sep))
@@ -87,8 +47,12 @@ for (F, TOL) in ((Float32, 0.2), (Float64, 0.0001), (BigFloat, 0.0001))
     end
 end
 
-# Test conversion with mixed floating types.
+# Test separation between coordinates and conversion with mixed floating types.
 c1 = ICRSCoords(e, pi/2)
+c5 = ICRSCoords(e, 1 + pi/2)
+@test separation(c1, c5) ≈ separation(c5, c1) ≈
+    separation(c1, convert(GalCoords, c5)) ≈
+    separation(convert(FK5Coords{1980}, c5), c1) ≈ 1
 for T in (GalCoords, FK5Coords{2000})
     c2 = convert(T{Float32}, c1)
     c3 = convert(T{Float64}, c1)
@@ -100,6 +64,8 @@ for T in (GalCoords, FK5Coords{2000})
     @test isapprox(lat(c3), lat(c4), rtol=sqrt(eps(Float64)))
     @test isapprox(lon(c2), lon(c3), rtol=sqrt(eps(Float32)))
     @test isapprox(lon(c3), lon(c4), rtol=sqrt(eps(Float64)))
+    c6 = convert(T, c5)
+    @test separation(c3, c6) ≈ separation(c6, c3) ≈ 1
 end
 
 println()