Merge pull request #46 from aplavin/ecliptic

JuliaAstro · Feb 11, 2023 · 07653d1 · 07653d1
2 parents 6c8849b + 04d5038
commit 07653d1
Show file tree

Hide file tree

Showing 3 changed files with 39 additions and 7 deletions.
diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
@@ -9,6 +9,7 @@ export AbstractSkyCoords,
        ICRSCoords,
        GalCoords,
        FK5Coords,
+       EclipticCoords,
        CartesianCoords,
        separation,
        position_angle,
@@ -101,11 +102,20 @@ function precess_from_j2000(equinox)
     zrotmat(-z) * yrotmat(theta) * zrotmat(-zeta)
 end
 
+function ecliptic_obliquity(y)
+    # https://github.com/JuliaAstro/AstroBase.jl/blob/master/src/EarthAttitude/obliquity.jl
+    T = (y - 2000.0) / 100
+    obl = @evalpoly(T, 84381.406, -46.836769, -0.0001831, 0.00200340, -0.000000576, -0.0000000434)
+    return deg2rad(obl / 3600)
+end
+
 # -----------------------------------------------------------------------------
 # Type-dependent methods
 
 lon(c::GalCoords) = c.l
 lat(c::GalCoords) = c.b
+lat(c::EclipticCoords) = c.lat
+lon(c::EclipticCoords) = c.lon
 lon(c::AbstractSkyCoords) = c.ra
 lat(c::AbstractSkyCoords) = c.dec
 
@@ -121,6 +131,13 @@ rotmat(::Type{<:ICRSCoords}, ::Type{<:ICRSCoords}) = I
 rotmat(::Type{<:GalCoords}, ::Type{<:GalCoords}) = I
 rotmat(::Type{<:FK5Coords{e}}, ::Type{<:FK5Coords{e}}) where {e} = I
 
+@generated rotmat(::Type{<:EclipticCoords{e}}, ::Type{<:FK5Coords{e}}) where {e} = xrotmat(ecliptic_obliquity(e))
+@generated rotmat(::Type{<:FK5Coords{e}}, ::Type{<:EclipticCoords{e}}) where {e} = xrotmat(-ecliptic_obliquity(e))
+@generated rotmat(::Type{<:EclipticCoords{e}}, ::Type{T}) where {e,T<:AbstractSkyCoords} = rotmat(EclipticCoords{e}, FK5Coords{e}) * rotmat(FK5Coords{e}, T)
+@generated rotmat(::Type{T}, ::Type{<:EclipticCoords{e}}) where {e,T<:AbstractSkyCoords} = rotmat(T, FK5Coords{e}) * rotmat(FK5Coords{e}, EclipticCoords{e})
+# disambiguation:
+@generated rotmat(::Type{<:EclipticCoords{e_to}}, ::Type{<:EclipticCoords{e_from}}) where {e_to,e_from} = rotmat(EclipticCoords{e_to}, FK5Coords{e_to}) * rotmat(FK5Coords{e_to}, EclipticCoords{e_from})
+
 @generated rotmat(::Type{<:FK5Coords{e1}}, ::Type{<:ICRSCoords}) where {e1} =
     precess_from_j2000(e1) * ICRS_TO_FK5J2000
 @generated rotmat(::Type{<:FK5Coords{e1}}, ::Type{<:GalCoords}) where {e1} =

diff --git a/src/types.jl b/src/types.jl
@@ -72,6 +72,18 @@ FK5Coords{e,F}(c::T) where {e,F,T<:AbstractSkyCoords} = convert(FK5Coords{e,F},
 constructorof(::Type{<:FK5Coords{e}}) where {e} = FK5Coords{e}
 
 
+struct EclipticCoords{e,T<:AbstractFloat} <: AbstractSkyCoords
+    lon::T
+    lat::T
+    EclipticCoords{e,T}(lon, lat) where {e,T<:AbstractFloat} = new(mod2pi(lon), lat)
+end
+EclipticCoords{e}(lon::T, lat::T) where {e,T<:AbstractFloat} = EclipticCoords{e,T}(lon, lat)
+EclipticCoords{e}(lon::Real, lat::Real) where {e} = EclipticCoords(promote(float(lon), float(lat))...)
+EclipticCoords{e}(c::AbstractSkyCoords) where {e} = convert(EclipticCoords{e}, c)
+EclipticCoords{e,F}(c::AbstractSkyCoords) where {e,F} = convert(EclipticCoords{e,F}, c)
+constructorof(::Type{<:EclipticCoords{e}}) where {e} = EclipticCoords{e}
+
+
 # Scalar coordinate conversions
 Base.convert(::Type{T}, c::T) where {T<:AbstractSkyCoords} = c
 
@@ -85,3 +97,4 @@ Base.:(==)(a::T, b::T) where {T<:AbstractSkyCoords} = lon(a) == lon(b) && lat(a)
 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::EclipticCoords{e}, b::EclipticCoords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
diff --git a/test/astropy.jl b/test/astropy.jl
@@ -12,8 +12,9 @@ const lats = pi .* (rand(rng, N) .- 0.5) # (-π, π)
 astropy_conversion(::Type{<:ICRSCoords}) = apc.ICRS
 astropy_conversion(::Type{<:FK5Coords{F}}) where {F} = apc.FK5(equinox="J$F")
 astropy_conversion(::Type{<:GalCoords}) = apc.Galactic
+astropy_conversion(::Type{<:EclipticCoords{F}}) where {F} = apc.GeocentricMeanEcliptic(equinox="J$F")
 
-function test_against_astropy(intype, outtype)
+function test_against_astropy(intype, outtype; atol=0)
     ## get julia values
     output_coords = map((lon, lat) -> convert(outtype, intype(lon, lat)), lons, lats)
     output_lons = map(lon, output_coords)
@@ -28,8 +29,8 @@ function test_against_astropy(intype, outtype)
     ap_output_lons = pyconvert(Vector, output_coord_list.spherical.lon.rad)
     ap_output_lats = pyconvert(Vector, output_coord_list.spherical.lat.rad)
 
-    @test output_lons ≈ ap_output_lons
-    @test output_lats ≈ ap_output_lats
+    @test output_lons ≈ ap_output_lons  atol=atol*√N
+    @test output_lats ≈ ap_output_lats  atol=atol*√N
 
 end
 
@@ -47,10 +48,11 @@ end
         FK5Coords{2000,F},
         FK5Coords{1975,F},
         GalCoords{F},
+        EclipticCoords{2000,F},
+        EclipticCoords{1975,F},
     )
-
-    @testset "$IN_SYS --> $OUT_SYS"for IN_SYS in systems, OUT_SYS in systems
-        test_against_astropy(IN_SYS, OUT_SYS)
-
+    @testset "$IN_SYS --> $OUT_SYS" for IN_SYS in systems, OUT_SYS in systems
+        atol = IN_SYS <: EclipticCoords || OUT_SYS <: EclipticCoords ? 1e-3 : 0
+        test_against_astropy(IN_SYS, OUT_SYS; atol)
     end
 end