Skip to content

Commit

Permalink
Constructors (#26)
Browse files Browse the repository at this point in the history 
* put types into their own file and add constructor from other type

* add tests and write code allowing unit conversion with pipe

* rebase branch and add version

* add conversion docs to documentation
  • Loading branch information
@mileslucas @giordano
mileslucas authored and giordano committed Nov 8, 2019
1 parent ac24a37 commit 8a56a00
Show file tree
Hide file tree
Showing 6 changed files with 148 additions and 108 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "SkyCoords"
uuid = "fc659fc5-75a3-5475-a2ea-3da92c065361"
authors = "Kyle Barbary, Mosé Giordano, and contributors"
version = "0.4.0"
version = "0.5.0-DEV"

[deps]
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Expand Down
29 changes: 29 additions & 0 deletions docs/src/api.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
# API/Reference

```@meta
DocTestSetup = :(using SkyCoords)
```

## Index

```@index
Expand All @@ -15,6 +19,31 @@ GalCoords
FK5Coords
```

## Conversion

To convert between types, there are three (equivalent) methods of doing so.

```jldoctest convsetup
julia> c1 = ICRSCoords(0., 0.)
ICRSCoords{Float64}(0.0, 0.0)
```

- using `convert`
```jldoctest convsetup
julia> convert(GalCoords, c1)
GalCoords{Float64}(1.6814027872278692, -1.0504884034813007)
```
- using constructors
```jldoctest convsetup
julia> GalCoords(c1)
GalCoords{Float64}(1.6814027872278692, -1.0504884034813007)
```
- using `|>`
```jldoctest convsetup
julia> c1 |> GalCoords
GalCoords{Float64}(1.6814027872278692, -1.0504884034813007)
```

## Functions

```@docs
Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ GalCoords{Float64}(1.6814027872278692, -1.0504884034813007)
julia> c2.l # Note that galactic coordinate fields are l, b
1.6814027872278692
julia> convert(FK5Coords{2000}, c1)
julia> c1 |> FK5Coords{2000} # Can use piping syntax for conversion
FK5Coords{2000,Float64}(1.1102233723050067e-7, 4.411803426976326e-8)
```
Expand Down
77 changes: 1 addition & 76 deletions src/SkyCoords.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,82 +11,7 @@ export AbstractSkyCoords,

import Base: convert

# -----------------------------------------------------------------------------
# Types

"""
The supertype for all sky coordinate systems.
"""
abstract type AbstractSkyCoords end

"""
ICRSCoords(ra, dec)
[International Celestial Reference System](https://en.wikipedia.org/wiki/International_Celestial_Reference_System)
This is the current standard adopted by the International Astronomical Union notably due to its high level of accuracy compared to standard equatorial coordinate systems. What sets this apart from [`FK5Coords`](@ref) is that it is completely defined using extragalactic radio sources rather than a geocentric frame, which means the reference frame will not change due to Earth's motion.
# Coordinates
- `ra` - Right ascension in radians (0, 2π)
- `dec` - Declination in radians (-π/2, π/2)
"""
struct ICRSCoords{T <: AbstractFloat} <: AbstractSkyCoords
ra::T
dec::T
ICRSCoords{T}(ra, dec) where {T <: AbstractFloat} = new(mod2pi(ra), dec)
end
ICRSCoords(ra::T, dec::T) where {T <: AbstractFloat} = ICRSCoords{T}(ra, dec)
ICRSCoords(ra::Real, dec::Real) = ICRSCoords(promote(float(ra), float(dec))...)


"""
GalCoords(l, b)
[Galactic Coordinate System](https://en.wikipedia.org/wiki/Galactic_coordinate_system)
This coordinate system is defined based on the projection of the Milky Way galaxy onto our celestial sphere, with (0, 0) being approximately the center of our galaxy.
# Coordinates
- `l` - Galactic longitude in radians (-π, π)
- `b` - Galactic latitude in radians (-π/2, π/2)
"""
struct GalCoords{T <: AbstractFloat} <: AbstractSkyCoords
l::T
b::T
GalCoords{T}(l, b) where {T <: AbstractFloat} = new(mod2pi(l), b)
end
GalCoords(l::T, b::T) where {T <: AbstractFloat} = GalCoords{T}(l, b)
GalCoords(l::Real, b::Real) = GalCoords(promote(float(l), float(b))...)


"""
FK5Coords{equinox}(ra, dec)
[Equatorial Coordinate System](https://en.wikipedia.org/wiki/Equatorial_coordinate_system)
This coordinate system maps the celestial sphere based on a geocentric observer. Historically the oldest, this coordinate system has been shown to be inaccurate due to its definitions based on the Earth, which has long-scale precession causing the reference frame to change. Because of this, an equinox must be provided (typically 2000, commonly known as J2000) which defines the reference frame.
# Coordinates
- `ra` - Right ascension in radians (0, 2π)
- `dec` - Declination in radians (-π/2, π/2)
"""
struct FK5Coords{e,T <: AbstractFloat} <: AbstractSkyCoords
ra::T
dec::T
FK5Coords{e,T}(ra, dec) where {T <: AbstractFloat,e} = new(mod2pi(ra), dec)
end
FK5Coords{e}(ra::T, dec::T) where {e,T <: AbstractFloat} = FK5Coords{e,T}(ra, dec)
FK5Coords{e}(ra::Real, dec::Real) where {e} =
FK5Coords{e}(promote(float(ra), float(dec))...)

# Scalar coordinate conversions
Base.convert(::Type{T}, c::T) where {T <: AbstractSkyCoords} = c

function Base.convert(::Type{T}, c::S) where {T <: AbstractSkyCoords,S <: AbstractSkyCoords}
r = rotmat(T, S) * coords2cart(c)
lon, lat = cart2coords(r)
T(lon, lat)
end
include("types.jl")

# -----------------------------------------------------------------------------
# Helper functions: Create rotation matrix about a given axis (x, y, z)
Expand Down
83 changes: 83 additions & 0 deletions src/types.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
# -----------------------------------------------------------------------------
# Types

"""
The supertype for all sky coordinate systems.
"""
abstract type AbstractSkyCoords end

"""
ICRSCoords(ra, dec)
[International Celestial Reference System](https://en.wikipedia.org/wiki/International_Celestial_Reference_System)
This is the current standard adopted by the International Astronomical Union notably due to its high level of accuracy compared to standard equatorial coordinate systems. What sets this apart from [`FK5Coords`](@ref) is that it is completely defined using extragalactic radio sources rather than a geocentric frame, which means the reference frame will not change due to Earth's motion.
# Coordinates
- `ra` - Right ascension in radians (0, 2π)
- `dec` - Declination in radians (-π/2, π/2)
"""
struct ICRSCoords{T <: AbstractFloat} <: AbstractSkyCoords
ra::T
dec::T
ICRSCoords{T}(ra, dec) where {T <: AbstractFloat} = new(mod2pi(ra), dec)
end
ICRSCoords(ra::T, dec::T) where {T <: AbstractFloat} = ICRSCoords{T}(ra, dec)
ICRSCoords(ra::Real, dec::Real) = ICRSCoords(promote(float(ra), float(dec))...)
ICRSCoords(c::T) where {T <: AbstractSkyCoords} = convert(ICRSCoords, c)
ICRSCoords{F}(c::T) where {F,T <: AbstractSkyCoords} = convert(ICRSCoords{F}, c)

"""
GalCoords(l, b)
[Galactic Coordinate System](https://en.wikipedia.org/wiki/Galactic_coordinate_system)
This coordinate system is defined based on the projection of the Milky Way galaxy onto our celestial sphere, with (0, 0) being approximately the center of our galaxy.
# Coordinates
- `l` - Galactic longitude in radians (-π, π)
- `b` - Galactic latitude in radians (-π/2, π/2)
"""
struct GalCoords{T <: AbstractFloat} <: AbstractSkyCoords
l::T
b::T
GalCoords{T}(l, b) where {T <: AbstractFloat} = new(mod2pi(l), b)
end
GalCoords(l::T, b::T) where {T <: AbstractFloat} = GalCoords{T}(l, b)
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)

"""
FK5Coords{equinox}(ra, dec)
[Equatorial Coordinate System](https://en.wikipedia.org/wiki/Equatorial_coordinate_system)
This coordinate system maps the celestial sphere based on a geocentric observer. Historically the oldest, this coordinate system has been shown to be inaccurate due to its definitions based on the Earth, which has long-scale precession causing the reference frame to change. Because of this, an equinox must be provided (typically 2000, commonly known as J2000) which defines the reference frame.
# Coordinates
- `ra` - Right ascension in radians (0, 2π)
- `dec` - Declination in radians (-π/2, π/2)
"""
struct FK5Coords{e,T <: AbstractFloat} <: AbstractSkyCoords
ra::T
dec::T
FK5Coords{e,T}(ra, dec) where {T <: AbstractFloat,e} = new(mod2pi(ra), dec)
end
FK5Coords{e}(ra::T, dec::T) where {e,T <: AbstractFloat} = FK5Coords{e,T}(ra, dec)
FK5Coords{e}(ra::Real, dec::Real) where {e} =
FK5Coords{e}(promote(float(ra), float(dec))...)
FK5Coords{e}(c::T) where {e,T <: AbstractSkyCoords} = convert(FK5Coords{e}, c)
FK5Coords{e,F}(c::T) where {e,F,T <: AbstractSkyCoords} = convert(FK5Coords{e,F}, c)


# Scalar coordinate conversions
Base.convert(::Type{T}, c::T) where {T <: AbstractSkyCoords} = c

function Base.convert(::Type{T}, c::S) where {T <: AbstractSkyCoords,S <: AbstractSkyCoords}
r = rotmat(T, S) * coords2cart(c)
lon, lat = cart2coords(r)
T(lon, lat)
end

Base.:(==)(a::T, b::T) where {T <: AbstractSkyCoords} = lon(a) == lon(b) && lat(a) == lat(b)
63 changes: 33 additions & 30 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,28 +14,35 @@ rad2arcsec(r) = 3600 * rad2deg(r)

# input coordinates
fname = joinpath(datapath, "input_coords.csv")
indata, inhdr = readdlm(fname, ','; header=true)
indata, inhdr = readdlm(fname, ','; header = true)

# Float32 has a large tolerance compared to Float64 and BigFloat, but here we
# are more interested in making sure that the infrastructure works for different
# floating types.
for (F, TOL) in ((Float32, 0.2), (Float64, 0.0001), (BigFloat, 0.0001))
println("Testing type ", F)
@testset "Testing $F" for (F, TOL) in ((Float32, 0.2), (Float64, 0.0001), (BigFloat, 0.0001))

for (insys, T) in (("icrs", ICRSCoords{F}), ("fk5j2000", FK5Coords{2000,F}),
("fk5j1975", FK5Coords{1975,F}), ("gal", GalCoords{F}))

c_in = T[T(indata[i, 1], indata[i, 2]) for i=1:size(indata,1)]
c_in = T[T(indata[i, 1], indata[i, 2]) for i = 1:size(indata, 1)]
for c in c_in
@test convert(T, c) == c
end

for (outsys, S) in (("icrs", ICRSCoords{F}), ("fk5j2000", FK5Coords{2000,F}),
("fk5j1975", FK5Coords{1975,F}), ("gal", GalCoords{F}))
(outsys == insys) && continue

c_out = S[convert(S, c) for c in c_in]

# Test pipe and constructor conversion
@test c_out == S[S(c) for c in c_in]
@test c_out == S[c |> S for c in c_in]

# Read in reference answers.
ref_fname = joinpath(datapath, "$(insys)_to_$(outsys).csv")
refdata, hdr = readdlm(ref_fname, ','; header=true)
c_ref = S[S(refdata[i, 1], refdata[i, 2]) for i=1:size(refdata,1)]
refdata, hdr = readdlm(ref_fname, ','; header = true)
c_ref = S[S(refdata[i, 1], refdata[i, 2]) for i = 1:size(refdata, 1)]

# compare
sep = separation.(c_out, c_ref)
Expand All @@ -49,28 +56,24 @@ for (F, TOL) in ((Float32, 0.2), (Float64, 0.0001), (BigFloat, 0.0001))
end

# Test separation between coordinates and conversion with mixed floating types.
c1 = ICRSCoords(ℯ, pi/2)
c5 = ICRSCoords(ℯ, 1 + pi/2)
@test separation(c1, c5) ≈ separation(c5, c1) ≈
separation(c1, convert(GalCoords, c5)) ≈
separation(convert(FK5Coords{1980}, c5), c1) ≈ 1
for T in (GalCoords, FK5Coords{2000})
c2 = convert(T{Float32}, c1)
c3 = convert(T{Float64}, c1)
c4 = convert(T{BigFloat}, c1)
@test typeof(c2) === T{Float32}
@test typeof(c3) === T{Float64}
@test typeof(c4) === T{BigFloat}
@test isapprox(lat(c2), lat(c3), rtol=sqrt(eps(Float32)))
@test isapprox(lat(c3), lat(c4), rtol=sqrt(eps(Float64)))
@test isapprox(lon(c2), lon(c3), rtol=sqrt(eps(Float32)))
@test isapprox(lon(c3), lon(c4), rtol=sqrt(eps(Float64)))
c6 = convert(T, c5)
@test separation(c3, c6) ≈ separation(c6, c3) ≈ 1
@testset "Separation" begin
c1 = ICRSCoords(ℯ, pi / 2)
c5 = ICRSCoords(ℯ, 1 + pi / 2)
@test separation(c1, c5) ≈ separation(c5, c1) ≈
separation(c1, convert(GalCoords, c5)) ≈
separation(convert(FK5Coords{1980}, c5), c1) ≈ 1
for T in (GalCoords, FK5Coords{2000})
c2 = convert(T{Float32}, c1)
c3 = convert(T{Float64}, c1)
c4 = convert(T{BigFloat}, c1)
@test typeof(c2) === T{Float32}
@test typeof(c3) === T{Float64}
@test typeof(c4) === T{BigFloat}
@test isapprox(lat(c2), lat(c3), rtol = sqrt(eps(Float32)))
@test isapprox(lat(c3), lat(c4), rtol = sqrt(eps(Float64)))
@test isapprox(lon(c2), lon(c3), rtol = sqrt(eps(Float32)))
@test isapprox(lon(c3), lon(c4), rtol = sqrt(eps(Float64)))
c6 = convert(T, c5)
@test separation(c3, c6) ≈ separation(c6, c3) ≈ 1
end
end

# Constructor of FK5Coords{1950} with non-float arguments
@test typeof(FK5Coords{1950}(1, big(2))) == FK5Coords{1950,BigFloat}

println()
println("All tests passed.")

0 comments on commit 8a56a00

Please sign in to comment.