Merge pull request #49 from LudwigBoess/master

Supergalactic Coordinates
JuliaAstro · Oct 16, 2024 · 6d1db5e · 6d1db5e
2 parents a900b63 + f2df76f
commit 6d1db5e
Show file tree

Hide file tree

Showing 10 changed files with 119 additions and 44 deletions.
diff --git a/README.md b/README.md
@@ -18,11 +18,12 @@ julia> Pkg.add("SkyCoords")
 ## Usage
 
 
-There are currently three supported coordinate systems. The following
+There are currently five supported coordinate systems. The following
 immutable types are used to represent coordinates in each system:
 
 - `ICRSCoords`: ICRS coordinates system
 - `GalCoords`: Galactic coordinates system
+- `SuperGalCoords`: Super-Galactic coordinate system
 - `FK5Coords`: FK5 coordinates system (with arbitrary equinox)
 - `EclipticCoords`: Ecliptic coordinates system
 

diff --git a/bench/bench.jl b/bench/bench.jl
@@ -18,6 +18,8 @@ for n in [1, 10, 100, 1000, 10_000, 100_000, 1_000_000]
         myprintln(io, "$n,fk5j2000,$t2")
         t3 = minimum(@benchmark(convert(FK5Coords{1975}, $c)).times) / 1e9
         myprintln(io, "$n,fk5j1975,$t3")
+        t4 = minimum(@benchmark(convert(SuperGalCoords, $c)).times) / 1e9
+        myprintln(io, "$n,supergalactic,$t4")
     else
         c = [ICRSCoords(2pi * rand(), pi * (rand() - 0.5)) for i = 1:n]
         t1 = minimum(@benchmark(convert(Vector{GalCoords{Float64}}, $c)).times) / 1e9
@@ -26,6 +28,8 @@ for n in [1, 10, 100, 1000, 10_000, 100_000, 1_000_000]
         myprintln(io, "$n,fk5j2000,$t2")
         t3 = minimum(@benchmark(convert(Vector{FK5Coords{1975,Float64}}, $c)).times) / 1e9
         myprintln(io, "$n,fk5j1975,$t3")
+        t4 = minimum(@benchmark(convert(Vector{SuperGalCoords{Float64}}, $c)).times) / 1e9
+        myprintln(io, "$n,supergalactic,$t4")
     end
     println()
 end

diff --git a/bench/bench.py b/bench/bench.py
@@ -22,9 +22,11 @@
     t1 = timeit.repeat("c.galactic", setup, repeat=3, number=nloops)
     t2 = timeit.repeat("c.transform_to('fk5')", setup, repeat=3, number=nloops)
     t3 = timeit.repeat("c.transform_to(FK5(equinox='J1975'))", setup, repeat=3, number=nloops)
+    t4 = timeit.repeat("c.supergalactic", setup, repeat=3, number=nloops)
     line = """{n:d},galactic,{t1:f}
 {n:d},fk5j2000,{t2:f}
 {n:d},fk5j1975,{t3:f}
+{n:d},supergalactic,{t3:f}
 """.format(n=n, t1=min(t1)/nloops, t2=min(t2)/nloops, t3=min(t3)/nloops)
     io.write(line)
     print(line,)

diff --git a/bench/julia_times.csv b/bench/julia_times.csv
@@ -1,22 +1,29 @@
 n,system,time
-1,galactic,6.690992835209827e-8
-1,fk5j2000,6.678118609406953e-8
-1,fk5j1975,6.694779938587513e-8
-10,galactic,7.389423076923077e-7
-10,fk5j2000,7.478246753246754e-7
-10,fk5j1975,7.499803921568628e-7
-100,galactic,6.81075e-6
-100,fk5j2000,6.79725e-6
-100,fk5j1975,6.77775e-6
-1000,galactic,9.0615e-5
-1000,fk5j2000,9.2829e-5
-1000,fk5j1975,9.2004e-5
-10000,galactic,0.00100854
-10000,fk5j2000,0.001019878
-10000,fk5j1975,0.001015691
-100000,galactic,0.010132403
-100000,fk5j2000,0.01024148
-100000,fk5j1975,0.010239616
-1000000,galactic,0.105910954
-1000000,fk5j2000,0.110236743
-1000000,fk5j1975,0.106532673
+1,galactic,9.070188284518828e-8
+1,fk5j2000,8.07080745341615e-8
+1,fk5j1975,8.1971875e-8
+1,supergalactic,9.603601694915254e-8
+10,galactic,8.433142857142857e-7
+10,fk5j2000,8.527272727272727e-7
+10,fk5j1975,8.526060606060606e-7
+10,supergalactic,8.368108108108109e-7
+100,galactic,8.163333333333334e-6
+100,fk5j2000,8.174666666666666e-6
+100,fk5j1975,8.192333333333335e-6
+100,supergalactic,8.367e-6
+1000,galactic,0.000101975
+1000,fk5j2000,0.000104378
+1000,fk5j1975,0.000104571
+1000,supergalactic,0.000103604
+10000,galactic,0.001081324
+10000,fk5j2000,0.001096162
+10000,fk5j1975,0.001100317
+10000,supergalactic,0.001094572
+100000,galactic,0.010836767
+100000,fk5j2000,0.010945083
+100000,fk5j1975,0.011005974
+100000,supergalactic,0.010936259
+1000000,galactic,0.108670411
+1000000,fk5j2000,0.110063537
+1000000,fk5j1975,0.110422317
+1000000,supergalactic,0.10979276
diff --git a/bench/plot_times.jl b/bench/plot_times.jl
@@ -1,4 +1,4 @@
-#!/usr/bin/env julia
+##!/usr/bin/env julia
 using CSV
 using StatsPlots
 

diff --git a/bench/python_times.csv b/bench/python_times.csv
@@ -1,22 +1,29 @@
 n,system,time
-1,galactic,0.002752
-1,fk5j2000,0.001413
-1,fk5j1975,0.001870
-10,galactic,0.002823
-10,fk5j2000,0.001398
-10,fk5j1975,0.001842
-100,galactic,0.003179
-100,fk5j2000,0.001619
-100,fk5j1975,0.001818
-1000,galactic,0.003131
-1000,fk5j2000,0.001613
-1000,fk5j1975,0.002028
-10000,galactic,0.005890
-10000,fk5j2000,0.003130
-10000,fk5j1975,0.003510
-100000,galactic,0.028635
-100000,fk5j2000,0.015479
-100000,fk5j1975,0.016745
-1000000,galactic,0.301621
-1000000,fk5j2000,0.161274
-1000000,fk5j1975,0.164234
+1,galactic,0.001644
+1,fk5j2000,0.001053
+1,fk5j1975,0.001588
+1,supergalactic,0.001588
+10,galactic,0.001670
+10,fk5j2000,0.001013
+10,fk5j1975,0.001603
+10,supergalactic,0.001603
+100,galactic,0.001711
+100,fk5j2000,0.001020
+100,fk5j1975,0.001641
+100,supergalactic,0.001641
+1000,galactic,0.002037
+1000,fk5j2000,0.001189
+1000,fk5j1975,0.001815
+1000,supergalactic,0.001815
+10000,galactic,0.005010
+10000,fk5j2000,0.002675
+10000,fk5j1975,0.003280
+10000,supergalactic,0.003280
+100000,galactic,0.033002
+100000,fk5j2000,0.016734
+100000,fk5j1975,0.017235
+100000,supergalactic,0.017235
+1000000,galactic,0.328802
+1000000,fk5j2000,0.165695
+1000000,fk5j1975,0.165792
+1000000,supergalactic,0.165792
diff --git a/docs/src/api.md b/docs/src/api.md
@@ -16,6 +16,7 @@ Pages = ["api.md"]
 AbstractSkyCoords
 ICRSCoords
 GalCoords
+SuperGalCoords
 FK5Coords
 EclipticCoords
 ```

diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
@@ -8,6 +8,7 @@ using StaticArrays
 export AbstractSkyCoords, 
        ICRSCoords,
        GalCoords,
+       SuperGalCoords,
        FK5Coords,
        EclipticCoords,
        CartesianCoords,
@@ -66,6 +67,20 @@ const GAL_TO_FK5J2000 = FK5J2000_TO_GAL'
 const GAL_TO_ICRS = FK5J2000_TO_ICRS * GAL_TO_FK5J2000
 const ICRS_TO_GAL = GAL_TO_ICRS'
 
+
+# Gal --> SuperGal
+# we use the same parameters as in astropy
+sgp_l = deg2rad(47.37)
+sgp_b = deg2rad(6.32)
+# rotation matrix compared to astropy
+const GAL_TO_SUPERGAL = RotZYZ(π/2, π/2 - sgp_b, π - sgp_l)
+const SUPERGAL_TO_GAL = GAL_TO_SUPERGAL'
+# SuperGal --> ICRS: chain through GAL
+const SUPERGAL_TO_ICRS = GAL_TO_ICRS * SUPERGAL_TO_GAL
+const ICRS_TO_SUPERGAL = SUPERGAL_TO_ICRS'
+# SuperGal --> FK5J2000: chain through GAL
+const SUPERGAL_TO_FK5J2000 = GAL_TO_FK5J2000 * SUPERGAL_TO_GAL
+const FK5J2000_TO_SUPERGAL = SUPERGAL_TO_FK5J2000'
 # FK5J2000 --> FK5{epoch}
 # Computes the precession matrix from J2000 to the given Julian equinox.
 # Expression from from Capitaine et al. 2003 as expressed in the USNO
@@ -99,6 +114,8 @@ end
 
 lon(c::GalCoords) = c.l
 lat(c::GalCoords) = c.b
+lon(c::SuperGalCoords) = c.l
+lat(c::SuperGalCoords) = c.b
 lat(c::EclipticCoords) = c.lat
 lon(c::EclipticCoords) = c.lon
 lon(c::AbstractSkyCoords) = c.ra
@@ -115,6 +132,12 @@ rotmat(::Type{<:ICRSCoords}, ::Type{<:GalCoords}) = GAL_TO_ICRS
 rotmat(::Type{<:ICRSCoords}, ::Type{<:ICRSCoords}) = I
 rotmat(::Type{<:GalCoords}, ::Type{<:GalCoords}) = I
 rotmat(::Type{<:FK5Coords{e}}, ::Type{<:FK5Coords{e}}) where {e} = I
+rotmat(::Type{<:SuperGalCoords}, ::Type{<:SuperGalCoords}) = I
+rotmat(::Type{<:GalCoords}, ::Type{<:SuperGalCoords}) = SUPERGAL_TO_GAL
+rotmat(::Type{<:SuperGalCoords}, ::Type{<:GalCoords}) = GAL_TO_SUPERGAL
+rotmat(::Type{<:SuperGalCoords}, ::Type{<:ICRSCoords}) = ICRS_TO_SUPERGAL
+rotmat(::Type{<:ICRSCoords}, ::Type{<:SuperGalCoords}) = SUPERGAL_TO_ICRS
+
 
 @generated rotmat(::Type{<:EclipticCoords{e}}, ::Type{<:FK5Coords{e}}) where {e} = RotX(-ecliptic_obliquity(e))
 @generated rotmat(::Type{<:FK5Coords{e}}, ::Type{<:EclipticCoords{e}}) where {e} = RotX(ecliptic_obliquity(e))
@@ -127,10 +150,14 @@ rotmat(::Type{<:FK5Coords{e}}, ::Type{<:FK5Coords{e}}) where {e} = I
     precess_from_j2000(e1) * ICRS_TO_FK5J2000
 @generated rotmat(::Type{<:FK5Coords{e1}}, ::Type{<:GalCoords}) where {e1} =
     precess_from_j2000(e1) * GAL_TO_FK5J2000
+@generated rotmat(::Type{<:FK5Coords{e1}}, ::Type{<:SuperGalCoords}) where {e1} =
+    precess_from_j2000(e1) * SUPERGAL_TO_FK5J2000
 @generated rotmat(::Type{<:ICRSCoords}, ::Type{<:FK5Coords{e2}}) where {e2} =
     FK5J2000_TO_ICRS * precess_from_j2000(e2)'
 @generated rotmat(::Type{<:GalCoords}, ::Type{<:FK5Coords{e2}}) where {e2} =
     FK5J2000_TO_GAL * precess_from_j2000(e2)'
+@generated rotmat(::Type{<:SuperGalCoords}, ::Type{<:FK5Coords{e2}}) where {e2} =
+    FK5J2000_TO_SUPERGAL * precess_from_j2000(e2)'
 @generated rotmat(::Type{<:FK5Coords{e1}}, ::Type{<:FK5Coords{e2}}) where {e1,e2} =
     precess_from_j2000(e1) * precess_from_j2000(e2)'
 

diff --git a/src/types.jl b/src/types.jl
@@ -48,6 +48,29 @@ GalCoords(l::Real, b::Real) = GalCoords(promote(float(l), float(b))...)
 GalCoords(c::T) where {T<:AbstractSkyCoords} = convert(GalCoords, c)
 GalCoords{F}(c::T) where {F,T<:AbstractSkyCoords} = convert(GalCoords{F}, c)
 
+"""
+    SuperGalCoords(l, b)
+
+[Supergalactic Coordinate System](https://en.wikipedia.org/wiki/Supergalactic_coordinate_system)
+
+The supergalactic plane is part of a reference frame for the supercluster of galaxies that contains the Milky Way galaxy.
+The supergalactic plane as so-far observed is more or less perpendicular to the plane of the Milky Way, the angle is 84.5 degrees. Viewed from the Earth, the plane traces a great circle across the sky through the constellations 
+
+# Coordinates
+- `l` - SuperGalCoords longitude in radians (-π, π)
+- `b` - SuperGalCoords latitude in radians (-π/2, π/2)
+"""
+struct SuperGalCoords{T<:AbstractFloat} <: AbstractSkyCoords
+    l::T
+    b::T
+    SuperGalCoords{T}(l, b) where {T<:AbstractFloat} = new(mod2pi(l), b)
+end
+SuperGalCoords(l::T, b::T) where {T<:AbstractFloat} = SuperGalCoords{T}(l, b)
+SuperGalCoords(l::Real, b::Real) = SuperGalCoords(promote(float(l), float(b))...)
+SuperGalCoords(c::T) where {T<:AbstractSkyCoords} = convert(SuperGalCoords, c)
+SuperGalCoords{F}(c::T) where {F,T<:AbstractSkyCoords} = convert(SuperGalCoords{F}, c)
+
+
 """
     FK5Coords{equinox}(ra, dec)
 
@@ -107,5 +130,6 @@ 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::SuperGalCoords, b::SuperGalCoords; 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,6 +12,7 @@ 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{<:SuperGalCoords}) = apc.Supergalactic
 astropy_conversion(::Type{<:EclipticCoords{F}}) where {F} = apc.GeocentricMeanEcliptic(equinox="J$F")
 
 function test_against_astropy(intype, outtype; atol=0)
@@ -48,6 +49,7 @@ end
         FK5Coords{2000,F},
         FK5Coords{1975,F},
         GalCoords{F},
+        SuperGalCoords{F},
         EclipticCoords{2000,F},
         EclipticCoords{1975,F},
     )