add simple projected coordinates

Author: aplavin

src/SkyCoords.jl

@@ -10,6 +10,7 @@ export AbstractSkyCoords,
        GalCoords,
        FK5Coords,
        CartesianCoords,
+       ProjectedCoords,
        separation,
        position_angle,
        offset,
@@ -18,6 +19,7 @@ export AbstractSkyCoords,
 
 include("types.jl")
 include("cartesian.jl")
+include("projected.jl")
 
 # -----------------------------------------------------------------------------
 # Helper functions: Create rotation matrix about a given axis (x, y, z)

src/projected.jl

@@ -0,0 +1,30 @@
+abstract type AbstractProjectedCoords <: AbstractSkyCoords end
+
+""" Projected coordinates with origin and offset.
+Approximation for small FoV. """
+struct ProjectedCoords{TC,T} <: AbstractProjectedCoords
+    origin::TC
+    offset::SVector{2,T}
+end
+
+origin(c::ProjectedCoords) = c.origin
+lon(c::AbstractProjectedCoords) = lon(origin(c)) + c.offset[1] / cos(lat(origin(c)))
+lat(c::AbstractProjectedCoords) = lat(origin(c)) + c.offset[2]
+
+Base.convert(::Type{T}, c::T) where {T<:AbstractProjectedCoords} = c
+Base.convert(::Type{TCto}, c::AbstractProjectedCoords) where {TCto<:AbstractSkyCoords} =
+    convert(TCto,
+        constructorof(typeof(origin(c)))(
+            lon(c), lat(c)
+        )
+    )
+
+Base.isapprox(a::ProjectedCoords, b::ProjectedCoords; kwargs...) = isapprox(origin(a), origin(b); kwargs...) && isapprox(a.offset, b.offset; kwargs...)
+
+function project(origin::AbstractSkyCoords, c::AbstractSkyCoords)
+    cc = convert(typeof(origin), c)
+    Δlon = _center_angle(lon(cc) - lon(origin))
+    offset = SVector(Δlon * cos(lat(origin)), lat(cc) - lat(origin))
+    ProjectedCoords(origin, offset)
+end
+

test/runtests.jl

@@ -3,6 +3,7 @@ using ConstructionBase: setproperties
 using DelimitedFiles
 using LinearAlgebra: normalize
 using SkyCoords
+using SkyCoords: project, origin
 using StableRNGs
 using Statistics
 using Test
@@ -15,6 +16,19 @@ rad2arcsec(r) = 3600 * rad2deg(r)
 # tests against astropy.coordinates
 include("astropy.jl")
 
+@testset "projected coords" begin
+    c0 = ICRSCoords(0.1, -0.2)
+    c1 = ICRSCoords(0.1 + 1e-5, -0.2 + 3e-5)
+    cp = project(c0, c1)::ProjectedCoords
+    @test origin(cp) == c0
+    @test cp.offset[1] ≈ 0.98 * 1e-5  rtol=1e-4
+    @test cp.offset[2] ≈ 3e-5
+    @test convert(ICRSCoords, cp) ≈ c1
+    @test convert(GalCoords, cp) ≈ convert(GalCoords, c1)
+    @test cp == cp
+    @test cp ≈ cp
+end
+
 # Test separation between coordinates and conversion with mixed floating types.
 @testset "Separation" begin
     c1 = ICRSCoords(ℯ, pi / 2)