position angle (#32)

* add position_angle * add position_angle to readme * put actual test into inferred test
JuliaAstro · Mar 25, 2020 · 27b4ece · 27b4ece
1 parent 12d19bd
commit 27b4ece
Show file tree

Hide file tree

Showing 5 changed files with 75 additions and 1 deletion.
diff --git a/README.md b/README.md
@@ -37,6 +37,8 @@ FK5Coords{2000,Float64}(1.1102233723050067e-7, 4.411803426976326e-8)
 julia> separation(c1, ICRSCoords(1., 0.)) # radians
 1.0
 
+julia> position_angle(c1, ICRSCoords(1, 0)) |> rad2deg
+90.0
 ```
 
 For more information, visit the [documentation](https://juliaastro.github.io/skycoords.jl/stable)

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -49,4 +49,5 @@ GalCoords{Float64}(1.6814027872278692, -1.0504884034813007)
 ```@docs
 SkyCoords.str2rad
 separation
+position_angle
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -84,6 +84,17 @@ julia> rad2deg(separation(mizar, alcor)) * 60 # Arcminutes
 11.8097230039347
 ```
 
+with an angle
+
+```jldoctest sep
+julia> position_angle(mizar, alcor) # radians
+1.244602401241819
+
+julia> position_angle(mizar, alcor) |> rad2deg # degrees
+71.31046476300408
+
+```
+
 ## Accuracy
 
 All the supported conversions have been compared to the results of

diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
@@ -3,7 +3,12 @@ __precompile__()
 module SkyCoords
 using StaticArrays
 
-export AbstractSkyCoords, ICRSCoords, GalCoords, FK5Coords, separation
+export AbstractSkyCoords, 
+       ICRSCoords,
+       GalCoords,
+       FK5Coords,
+       separation,
+       position_angle
 
 include("types.jl")
 
@@ -224,4 +229,33 @@ separation(c1::T, c2::T) where {T<:AbstractSkyCoords} =
 separation(c1::T1, c2::T2) where {T1<:AbstractSkyCoords,T2<:AbstractSkyCoords} =
     separation(c1, convert(T1, c2))
 
+
+"""
+    position_angle(c1::AbstractSkyCoords, c2::AbstractSkyCoords) -> angle
+
+Return position angle between two sky coordinates, in positive radians.
+
+# Examples
+```jldoctest
+julia> c1 = ICRSCoords(0, 0); c2 = ICRSCoords(deg2rad(1), 0);
+
+julia> position_angle(c1, c2) |> rad2deg
+90.0
+```
+"""
+position_angle(c1::T, c2::T) where {T <: AbstractSkyCoords} = _position_angle(lon(c1), lat(c1), lon(c2), lat(c2))
+position_angle(c1::T1, c2::T2) where {T1 <: AbstractSkyCoords,T2 <: AbstractSkyCoords} = position_angle(c1, convert(T1, c2))
+
+
+function _position_angle(λ1, ϕ1, λ2, ϕ2)
+    sin_Δλ, cos_Δλ = sincos(λ2 - λ1)
+    sin_ϕ1, cos_ϕ1 = sincos(ϕ1)
+    sin_ϕ2, cos_ϕ2 = sincos(ϕ2)
+
+    x = sin_ϕ2 * cos_ϕ1 - cos_ϕ2 * sin_ϕ1 * cos_Δλ
+    y = sin_Δλ * cos_ϕ2
+
+    return mod2pi(atan(y, x))
+end
+
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -137,3 +137,29 @@ end
         @test C("12:0:0", "90:0:0") == C(π, π / 2)
     end
 end
+
+# Test separation between coordinates and conversion with mixed floating types.
+@testset "Position Angles" begin
+    c1 = ICRSCoords(0, 0)
+    c2 = ICRSCoords(deg2rad(1), 0)
+
+    # interface
+    @test @inferred position_angle(c1, c2) ≈ @inferred position_angle(c1, c2 |> GalCoords)
+    @test position_angle(c1, c2) == position_angle(c1, c2 |> GalCoords)
+
+    # accuracy
+    @test position_angle(c1, c2) ≈ π / 2
+
+    c3 = ICRSCoords(deg2rad(1), deg2rad(0.1))
+    @test position_angle(c1, c3) < π / 2
+
+    c4 = ICRSCoords(0, deg2rad(1))
+    @test position_angle(c1, c4) ≈ 0
+
+    # types
+    for T in [ICRSCoords, GalCoords, FK5Coords{2000}]
+        c1 = T(0, 0)
+        c2 = T(deg2rad(1), 0)
+        @test position_angle(c1, c2) ≈ π / 2
+    end
+end
-Original file line number
+Diff line change
@@ Expand Up @@
     julia> separation(c1, ICRSCoords(1., 0.)) # radians
 .0
+    julia> position_angle(c1, ICRSCoords(1, 0)) |> rad2deg
+.0
     ```
     For more information, visit the [documentation](https://juliaastro.github.io/skycoords.jl/stable)
@@ Expand Down @@