Add offset function (#34)

.travis.yml

@@ -1,30 +1,36 @@
 language: julia
+
 os:
-  - osx
   - linux
+  - osx
+
 julia:
   - 1.0
-  - 1.1
-  - 1.2
-  - 1.3
+  - 1.4
   - nightly
 
 notifications:
   email: false
 
-matrix:
+jobs:
   allow_failures:
     - julia: nightly
-
-codecov: true
-
-jobs:
   include:
-    - stage: Documentation
-      julia: 1.2
-      script: julia --project=docs -e '
-        using Pkg;
-        Pkg.develop(PackageSpec(path=pwd()));
-        Pkg.instantiate();
-        include("docs/make.jl");'
+    - stage: "Documentation"
+      julia: 1.3
+      os: linux
+      script: |
+        julia --project=docs -e '
+          using Pkg
+          Pkg.develop(PackageSpec(path=pwd()))
+          Pkg.instantiate()
+          include("docs/make.jl")'
       after_success: skip
+
+after_success:
+  - |
+    julia -e '
+      using Pkg
+      Pkg.add("Coverage")
+      using Coverage
+      Codecov.submit(process_folder())'

Project.toml

@@ -1,7 +1,7 @@
 name = "SkyCoords"
 uuid = "fc659fc5-75a3-5475-a2ea-3da92c065361"
 authors = "Kyle Barbary, Mosé Giordano, and contributors"
-version = "0.5.0"
+version = "0.6.0"
 
 [deps]
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

appveyor.yml

@@ -1,9 +1,7 @@
 environment:
   matrix:
   - julia_version: 1
-  - julia_version: 1.1
-  - julia_version: 1.2
-  - julia_version: 1.3
+  - julia_version: 1.4
   - julia_version: nightly
 
 platform:

docs/make.jl

@@ -11,4 +11,4 @@ makedocs(
     strict = true,
 )
 
-deploydocs(repo = "github.com/JuliaAstro/SkyCoords.jl.git")
+deploydocs(repo = "github.com/JuliaAstro/SkyCoords.jl.git", push_preview=true)

docs/src/api.md

@@ -50,4 +50,5 @@ GalCoords{Float64}(1.6814027872278692, -1.0504884034813007)
 SkyCoords.str2rad
 separation
 position_angle
+offset
 ```

docs/src/index.md

@@ -78,10 +78,10 @@ Their angular separation is given by
 
 ```jldoctest sep
 julia> separation(mizar, alcor) # Radians
-0.0034353091694529292
+0.0034353091694529297
 
 julia> rad2deg(separation(mizar, alcor)) * 60 # Arcminutes
-11.8097230039347
+11.809723003934701
 ```
 
 with an angle

src/SkyCoords.jl

@@ -8,7 +8,8 @@ export AbstractSkyCoords,
        GalCoords,
        FK5Coords,
        separation,
-       position_angle
+       position_angle,
+       offset
 
 include("types.jl")
 
@@ -258,4 +259,68 @@ function _position_angle(λ1, ϕ1, λ2, ϕ2)
     return mod2pi(atan(y, x))
 end
 
+"""
+    offset(::AbstractSkyCoords, separation, pa) -> coordinate
+
+Offset a coordinate by a given angular separation, `separation`, in radians and position angle, `pa`, in radians.
+
+Uses the sine and cosine rules in spherical coordinates with corrections for the antipodes. Returns a sky coordinate of the same type as input.
+
+# Examples
+```jldoctest
+julia> c1 = ICRSCoords(0, 0);
+
+julia> c2 = offset(c1, deg2rad(1), deg2rad(90))
+ICRSCoords{Float64}(0.017453292519943295, 1.0686516840418957e-18)
+
+julia> offset(c1, c2) .|> rad2deg
+(1.0, 90.0)
+```
+
+# See Also
+* [`separation`](@ref), [`position_angle`](@ref)
+"""
+offset(c::T, sep, pa) where T <: AbstractSkyCoords = T(_offset(lon(c), lat(c), sep, pa)...)
+
+"""
+    offset(::AbstractSkyCoords, AbstractSkyCoords) -> angle, angle
+
+Return the separation and position angle in radians between two sky coordinates.
+
+# Examples
+```jldoctest
+julia> c1 = ICRSCoords(0, 0); c2 = ICRSCoords(deg2rad(1), 0);
+
+julia> offset(c1, c2) .|> rad2deg
+(1.0, 90.0)
+```
+
+# See Also
+* [`separation`](@ref), [`position_angle`](@ref)
+"""
+offset(c1::AbstractSkyCoords, c2::AbstractSkyCoords) = separation(c1, c2), position_angle(c1, c2)
+
+#= use the cosine rule in spherical geometry with three points, the north pole, the starting point, 
+and the final point.
+angles: (change in lon), (position angle), (-1/position angle)
+sides: (separation), (final co-latitude), (starting co-latitude)
+=#
+function _offset(λ, ϕ, separation, pa)
+    sin_a, cos_a = sincos(separation)
+    cos_c, sin_c = sincos(ϕ)
+    sin_B, cos_B = sincos(pa)
+
+    # solving cosine rule
+    cos_b = cos_c * cos_a + sin_c * sin_a * cos_B
+
+    # solving sine rule
+    xsin_A = sin_a * sin_B * sin_c
+    xcos_A = cos_a - cos_b * cos_c
+
+    # correction for antipodes, otherwise atan2
+    ang = sin_c < 1e-12 ? π/2 + cos_c * (π/2 - pa) : atan(xsin_A, xcos_A)
+
+    return mod2pi(λ + ang), asin(cos_b)
+end
+
 end # module

test/runtests.jl

@@ -163,3 +163,55 @@ end
         @test position_angle(c1, c2) ≈ π / 2
     end
 end
+
+
+
+@testset "Offset ($T1, $T2)" for T1 in [ICRSCoords, GalCoords, FK5Coords{2000}], T2 in [ICRSCoords, GalCoords, FK5Coords{2000}]
+    # simple integration tests, depend that separation and position_angle are accurate
+    c1s = [
+        T1(0, -π/2), # south pole
+        T1(0, π/2), # north pole
+        T1(deg2rad(1), deg2rad(2))
+    ]
+    c2 = T2(deg2rad(5), deg2rad(10))
+
+    for c1 in c1s
+        sep, pa = @inferred offset(c1, c2)
+        test_c2 = @inferred offset(c1, sep, pa)
+        @test test_c2 isa T1
+        test_c2 = T2(test_c2) 
+        @test lon(test_c2) ≈ lon(c2)
+        @test lat(test_c2) ≈ lat(c2)
+    end
+
+    # specific cases to cover special cases.
+    c1 = T1(0, deg2rad(89))
+    for (pa, sep) in [(0, 2), (180, 358)]
+        sep = deg2rad(sep)
+        pa = deg2rad(pa)
+        c2 = offset(c1, sep, pa)
+        @test lon(c2) |> rad2deg ≈ 180
+        @test lat(c2) |> rad2deg ≈ 89
+
+        c2 = offset(c1, 2sep, pa)
+        @test lon(c2) |> rad2deg ≈ 180
+        @test lat(c2) |> rad2deg ≈ 87
+    end
+
+    # verify antipode
+    c1 = T1(deg2rad(10), deg2rad(47))
+    for pa in range(0, stop=377, length=10)
+        c2 = offset(c1, deg2rad(180), deg2rad(pa))
+        @test lon(c2) |> rad2deg ≈ 190
+        @test lat(c2) |> rad2deg ≈ -47
+
+        c2 = offset(c1, deg2rad(360), deg2rad(pa))
+        @test lon(c2) |> rad2deg ≈ 10
+        @test lat(c2) |> rad2deg ≈ 47
+    end
+
+    c1 = T1(deg2rad(10), deg2rad(60))
+    c2 = offset(c1, deg2rad(1), deg2rad(90))
+    @test 11.9 < lon(c2) |> rad2deg < 12.0
+    @test 59.9 < lat(c2) |> rad2deg < 60.0
+end