From 7f5297d68718ae7453e6b0cead7ecbb0ee53eae2 Mon Sep 17 00:00:00 2001
From: Alexander Plavin <alexander@plav.in>
Date: Thu, 10 Nov 2022 16:39:56 +0300
Subject: [PATCH 1/3] fix isapprox

---
 src/types.jl     | 8 +++++---
 test/runtests.jl | 3 +++
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/types.jl b/src/types.jl
index 4abdfc5..46e89db 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -82,6 +82,8 @@ function Base.convert(::Type{T}, c::S) where {T<:AbstractSkyCoords,S<:AbstractSk
 end
 
 Base.:(==)(a::T, b::T) where {T<:AbstractSkyCoords} = lon(a) == lon(b) && lat(a) == lat(b)
-Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
+Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(_center_angle(lon(b) - lon(a)) + lon(a), lat(b)); kwargs...)
+Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(_center_angle(lon(b) - lon(a)) + lon(a), lat(b)); kwargs...)
+Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(_center_angle(lon(b) - lon(a)) + lon(a), lat(b)); kwargs...)
+
+_center_angle(x) = mod2pi(x + π) - π
diff --git a/test/runtests.jl b/test/runtests.jl
index 226762a..d4bf943 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -177,6 +177,9 @@ end
         @test !(c1 ≈ c4)
         @test c1 ≈ c2  rtol=1e-3
         @test c1 ≈ c4  rtol=1e-3
+
+        @test ICRSCoords(eps(), 1) ≈ ICRSCoords(0, 1)
+        @test ICRSCoords(eps(), 1) ≈ ICRSCoords(-eps(), 1)
     end
  
     @test_broken (!(ICRSCoords(1, 2) ≈ FK5Coords{2000}(1, 2)); true)

From 17fdff608cac021ea8f30368cb9f338f59a60baa Mon Sep 17 00:00:00 2001
From: Alexander Plavin <alexander@plav.in>
Date: Thu, 10 Nov 2022 16:47:00 +0300
Subject: [PATCH 2/3] add simple projected coordinates

---
 src/SkyCoords.jl |  2 ++
 src/projected.jl | 30 ++++++++++++++++++++++++++++++
 test/runtests.jl | 14 ++++++++++++++
 3 files changed, 46 insertions(+)
 create mode 100644 src/projected.jl

diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
index 9bbf7a8..346eba5 100644
--- a/src/SkyCoords.jl
+++ b/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)
diff --git a/src/projected.jl b/src/projected.jl
new file mode 100644
index 0000000..8376f85
--- /dev/null
+++ b/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
+
diff --git a/test/runtests.jl b/test/runtests.jl
index d4bf943..7cce13b 100644
--- a/test/runtests.jl
+++ b/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)

From 56fbeaeb31d25917af53f3957e0fc4a7316ebd38 Mon Sep 17 00:00:00 2001
From: Alexander Plavin <alexander@plav.in>
Date: Tue, 7 Feb 2023 18:03:00 +0300
Subject: [PATCH 3/3] remove _center_angle

---
 src/projected.jl |  2 +-
 src/types.jl     | 17 ++++++++++++-----
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/src/projected.jl b/src/projected.jl
index 8376f85..474bd46 100644
--- a/src/projected.jl
+++ b/src/projected.jl
@@ -23,7 +23,7 @@ Base.isapprox(a::ProjectedCoords, b::ProjectedCoords; kwargs...) = isapprox(orig
 
 function project(origin::AbstractSkyCoords, c::AbstractSkyCoords)
     cc = convert(typeof(origin), c)
-    Δlon = _center_angle(lon(cc) - lon(origin))
+    Δlon = rem2pi(lon(cc) - lon(origin), RoundNearest)
     offset = SVector(Δlon * cos(lat(origin)), lat(cc) - lat(origin))
     ProjectedCoords(origin, offset)
 end
diff --git a/src/types.jl b/src/types.jl
index 46e89db..f1556d0 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -82,8 +82,15 @@ function Base.convert(::Type{T}, c::S) where {T<:AbstractSkyCoords,S<:AbstractSk
 end
 
 Base.:(==)(a::T, b::T) where {T<:AbstractSkyCoords} = lon(a) == lon(b) && lat(a) == lat(b)
-Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(_center_angle(lon(b) - lon(a)) + lon(a), lat(b)); kwargs...)
-Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(_center_angle(lon(b) - lon(a)) + lon(a), lat(b)); kwargs...)
-Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(_center_angle(lon(b) - lon(a)) + lon(a), lat(b)); kwargs...)
-
-_center_angle(x) = mod2pi(x + π) - π
+Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...)
+Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...)
+Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...)