switch to using Rotations.jl

JuliaAstro · Dec 29, 2022 · df31123 · df31123
1 parent 4c2f0fe
commit df31123
Show file tree

Hide file tree

Showing 3 changed files with 28 additions and 41 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,16 +1,18 @@
 name = "SkyCoords"
 uuid = "fc659fc5-75a3-5475-a2ea-3da92c065361"
 authors = "Kyle Barbary, Mosé Giordano, and contributors"
-version = "1.1.0"
+version = "1.1.1"
 
 [deps]
 AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 AstroAngles = "0.1"
 ConstructionBase = "1"
+Rotations = "1"
 StaticArrays = "0.8, 0.9, 1"
 julia = "1.3"
diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
@@ -1,9 +1,9 @@
-__precompile__()
-
 module SkyCoords
-using StaticArrays
-using LinearAlgebra: I, norm
+
 import ConstructionBase: constructorof
+using LinearAlgebra: I, norm
+using Rotations
+using StaticArrays
 
 export AbstractSkyCoords, 
        ICRSCoords,
@@ -22,21 +22,6 @@ include("cartesian.jl")
 # -----------------------------------------------------------------------------
 # Helper functions: Create rotation matrix about a given axis (x, y, z)
 
-function xrotmat(angle)
-    s, c = sincos(angle)
-    SMatrix{3,3}(1, 0, 0, 0, c, -s, 0, s, c)
-end
-
-function yrotmat(angle)
-    s, c = sincos(angle)
-    SMatrix{3,3}(c, 0, s, 0, 1, 0, -s, 0, c)
-end
-
-function zrotmat(angle)
-    s, c = sincos(angle)
-    SMatrix{3,3}(c, -s, 0, s, c, 0, 0, 0, 1)
-end
-
 # (lon, lat) -> [x, y, z] unit vector
 function coords2cart(lon, lat)
     sinlon, coslon = sincos(lon)
@@ -45,16 +30,21 @@ function coords2cart(lon, lat)
 end
 
 # [x, y, z] unit vector -> (lon, lat)
-cart2coords(v) = atan(v[2], v[1]), atan(v[3], sqrt(v[1] * v[1] + v[2] * v[2]))
+function cart2coords(v)
+    lon = atan(v[begin + 1], v[begin])
+    xy_norm = sqrt(v[begin]^2 + v[begin + 1]^2)
+    lat = atan(v[begin + 2], xy_norm)
+    return lon, lat
+end
 
 # -----------------------------------------------------------------------------
 # Constant rotation matricies and precession matrix function
 
 # ICRS --> FK5 at J2000 (See USNO Circular 179, section 3.5)
-eta0 = deg2rad(-19.9 / 3600000.0)
-xi0 = deg2rad(9.1 / 3600000.0)
-da0 = deg2rad(-22.9 / 3600000.0)
-const ICRS_TO_FK5J2000 = xrotmat(-eta0) * yrotmat(xi0) * zrotmat(da0)
+eta0 = deg2rad(-19.9 / 3600000)
+xi0 = deg2rad(9.1 / 3600000)
+da0 = deg2rad(-22.9 / 3600000)
+const ICRS_TO_FK5J2000 = RotXYZ(eta0, -xi0, -da0)
 const FK5J2000_TO_ICRS = ICRS_TO_FK5J2000'
 
 # FK5J2000 --> Gal
@@ -68,8 +58,7 @@ const FK5J2000_TO_ICRS = ICRS_TO_FK5J2000'
 ngp_fk5j2000_ra = deg2rad(192.8594812065348)
 ngp_fk5j2000_dec = deg2rad(27.12825118085622)
 lon0_fk5j2000 = deg2rad(122.9319185680026)
-const FK5J2000_TO_GAL = (zrotmat(pi - lon0_fk5j2000) *
-                         yrotmat(pi / 2.0 - ngp_fk5j2000_dec) * zrotmat(ngp_fk5j2000_ra))
+const FK5J2000_TO_GAL = RotZYZ(lon0_fk5j2000 - π, ngp_fk5j2000_dec - π/2, -ngp_fk5j2000_ra)
 const GAL_TO_FK5J2000 = FK5J2000_TO_GAL'
 
 # Gal --> ICRS: simply chain through FK5J2000
@@ -80,25 +69,21 @@ const ICRS_TO_GAL = GAL_TO_ICRS'
 # Computes the precession matrix from J2000 to the given Julian equinox.
 # Expression from from Capitaine et al. 2003 as expressed in the USNO
 # Circular 179.  This should match the IAU 2006 standard from SOFA.
-const pzeta = [2.650545, 2306.083227, 0.2988499, 0.01801828, -0.000005971, -0.0000003173]
-const pz = [-2.650545, 2306.077181, 1.0927348, 0.01826837, -0.000028596, -0.0000002904]
-const ptheta = [0.0, 2004.191903, -0.4294934, -0.04182264, -0.000007089, -0.0000001274]
+const pzeta = SA[2.650545, 2306.083227, 0.2988499, 0.01801828, -0.000005971, -0.0000003173]
+const pz = SA[-2.650545, 2306.077181, 1.0927348, 0.01826837, -0.000028596, -0.0000002904]
+const ptheta = SA[0.0, 2004.191903, -0.4294934, -0.04182264, -0.000007089, -0.0000001274]
+
 function precess_from_j2000(equinox)
-    t = (equinox - 2000.0) / 100.0
+    t = (equinox - 2000) / 100
     tn = 1.0
-    zeta = pzeta[1]
-    z = pz[1]
-    theta = ptheta[1]
-    for i = 2:6
-        tn *= t
+    zeta = z = theta = 0.0
+    for i in eachindex(pzeta, pz, ptheta)
         zeta += pzeta[i] * tn
         z += pz[i] * tn
         theta += ptheta[i] * tn
+        tn *= t
     end
-    zeta = deg2rad(zeta / 3600.0)
-    z = deg2rad(z / 3600.0)
-    theta = deg2rad(theta / 3600.0)
-    zrotmat(-z) * yrotmat(theta) * zrotmat(-zeta)
+    return RotZYZ(deg2rad(z / 3600), -deg2rad(theta / 3600), deg2rad(zeta / 3600))
 end
 
 # -----------------------------------------------------------------------------

diff --git a/src/cartesian.jl b/src/cartesian.jl
@@ -4,7 +4,7 @@ end
 
 CartesianCoords{TC}(args::Real...) where {TC} = CartesianCoords{TC}(SVector(float.(args)...))
 CartesianCoords{TC}(vec::AbstractVector{TF}) where {TC, TF} = CartesianCoords{TC, TF}(vec)
-CartesianCoords(c::AbstractSkyCoords) where {TC} = convert(CartesianCoords, c)
+CartesianCoords(c::AbstractSkyCoords) = convert(CartesianCoords, c)
 CartesianCoords{TC}(c::AbstractSkyCoords) where {TC} = convert(CartesianCoords{TC}, c)
 CartesianCoords{TC,TF}(c::AbstractSkyCoords) where {TC<:AbstractSkyCoords,TF} = convert(CartesianCoords{TC,TF}, c)
 constructorof(::Type{<:CartesianCoords{TC}}) where {TC} = CartesianCoords{TC}