From bb432f2d064dd976e9d5b7b8200915224bfd5f5e Mon Sep 17 00:00:00 2001
From: xzackli <xzackli@gmail.com>
Date: Mon, 24 Jan 2022 01:56:01 -0500
Subject: [PATCH 1/5] add fullsky_geometry creator

---
 src/Pixell.jl    |  4 +++-
 src/enmap.jl     |  6 +++---
 src/enmap_ops.jl | 27 +++++++++++++++++++++++++++
 3 files changed, 33 insertions(+), 4 deletions(-)
 create mode 100644 src/enmap_ops.jl

diff --git a/src/Pixell.jl b/src/Pixell.jl
index a94fa8f..9263e1f 100644
--- a/src/Pixell.jl
+++ b/src/Pixell.jl
@@ -6,7 +6,9 @@ using FFTW
 using Printf
 
 include("enmap.jl")
+include("enmap_ops.jl")
 
-export Enmap
+export Enmap, CarProjection
+export fullsky_geometry
 
 end
diff --git a/src/enmap.jl b/src/enmap.jl
index a7c15e1..3d1ea10 100644
--- a/src/enmap.jl
+++ b/src/enmap.jl
@@ -1,7 +1,7 @@
 using Base: ViewIndex, @propagate_inbounds, AbstractCartesianIndex
 
-abstract type MapProjection end
-abstract type EquiCylProjection <: MapProjection end  # equidistant cylindrical projection
+abstract type AbstractMapProjection end
+abstract type EquiCylProjection <: AbstractMapProjection end  # equidistant cylindrical projection
 
 struct CarProjection <: EquiCylProjection end  # plate carrée
 struct CeaProjection <: EquiCylProjection end  # cylindrical equal area
@@ -13,7 +13,7 @@ It only implements the subset of Base.Array operations which are common on maps.
 You should work with the data directly using `enmap_instance.data` if you need
 additional Array functions.
 """
-struct Enmap{T,N,AA<:AbstractArray,P<:MapProjection} <: AbstractArray{T,N}
+struct Enmap{T,N,AA<:AbstractArray,P<:AbstractMapProjection} <: AbstractArray{T,N}
     data::AA  # some kind of abstract array
     wcs::WCSTransform  # WCS object from WCS.jl
 end
diff --git a/src/enmap_ops.jl b/src/enmap_ops.jl
new file mode 100644
index 0000000..121e4a2
--- /dev/null
+++ b/src/enmap_ops.jl
@@ -0,0 +1,27 @@
+
+function fullsky_geometry(proj::Type{<:CarProjection}, res; shape = nothing, dims = ())
+    if isnothing(shape)
+        shape = (round.(Int, (2π, π) ./ res .+ (0, 1)))  # CAR has pixels on poles
+    end
+    nx, ny = shape
+    resx, resy = res
+
+    @assert abs(resx * nx - 2π) < 1e-8 "Horizontal resolution does not evenly divide the sky; this is required for SHTs."
+    @assert abs(resy * (ny - 1) - π) < 1e-8 "Vertical resolution does not evenly divide the sky; this is required for SHTs."
+
+    # Note the reference point is shifted by half a pixel to keep
+    # the grid in bounds, from ra=180+cdelt/2 to ra=-180+cdelt/2.
+    wcs = WCSTransform(2;
+        cdelt = [-360.0 / nx, 180.0 / (ny - 1)],
+        ctype = ["RA---CAR", "DEC--CAR"],
+        crpix = [floor(nx / 2) + 0.5, (ny + 1) / 2],
+        crval = [resy * 90 / π, 0])
+
+    return (nx, ny, dims...), wcs
+end
+
+fullsky_geometry(proj::Type{<:CarProjection}, res::Number; shape = nothing, dims = ()) =
+    fullsky_geometry(proj, (res, res); shape = shape, dims = dims)
+
+fullsky_geometry(res; shape = nothing, dims = ()) =
+    fullsky_geometry(CarProjection, res; shape = shape, dims = dims)

From ed3e4d0dd16ecb9a1ca14dfba776f43156a6b84a Mon Sep 17 00:00:00 2001
From: xzackli <xzackli@gmail.com>
Date: Mon, 24 Jan 2022 02:13:35 -0500
Subject: [PATCH 2/5] make CarProjection abstract and have more concrete
 projection

---
 src/Pixell.jl    | 2 +-
 src/enmap.jl     | 9 +++++----
 src/enmap_ops.jl | 6 +++---
 3 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/Pixell.jl b/src/Pixell.jl
index 9263e1f..ec5335e 100644
--- a/src/Pixell.jl
+++ b/src/Pixell.jl
@@ -8,7 +8,7 @@ using Printf
 include("enmap.jl")
 include("enmap_ops.jl")
 
-export Enmap, CarProjection
+export Enmap, CarClenshawCurtis
 export fullsky_geometry
 
 end
diff --git a/src/enmap.jl b/src/enmap.jl
index 3d1ea10..97f84ef 100644
--- a/src/enmap.jl
+++ b/src/enmap.jl
@@ -2,9 +2,10 @@ using Base: ViewIndex, @propagate_inbounds, AbstractCartesianIndex
 
 abstract type AbstractMapProjection end
 abstract type EquiCylProjection <: AbstractMapProjection end  # equidistant cylindrical projection
+abstract type CarProjection <: EquiCylProjection end          # plate carrée
+
+struct CarClenshawCurtis <: CarProjection end      # plate carrée with pixels on poles and equator
 
-struct CarProjection <: EquiCylProjection end  # plate carrée
-struct CeaProjection <: EquiCylProjection end  # cylindrical equal area
 
 """
 Map type, contains an AbstractArray and a WCS object, but behaves like the
@@ -22,9 +23,9 @@ end
 function Enmap(data::A, wcs, ::Type{PT}) where {A<:AbstractArray,PT}
     Enmap{eltype(A),ndims(A),A,PT}(data, wcs)
 end
-# create CAR maps by default
+# create CAR (Clenshaw-Curtis variant) maps by default
 function Enmap(data::A, wcs) where {A<:AbstractArray}
-    Enmap{eltype(A),ndims(A),A,CarProjection}(data, wcs)
+    Enmap{eltype(A),ndims(A),A,CarClenshawCurtis}(data, wcs)
 end
 
 Base.parent(x::Enmap) = x.data
diff --git a/src/enmap_ops.jl b/src/enmap_ops.jl
index 121e4a2..c10e19c 100644
--- a/src/enmap_ops.jl
+++ b/src/enmap_ops.jl
@@ -1,5 +1,5 @@
 
-function fullsky_geometry(proj::Type{<:CarProjection}, res; shape = nothing, dims = ())
+function fullsky_geometry(proj::Type{<:CarClenshawCurtis}, res; shape = nothing, dims = ())
     if isnothing(shape)
         shape = (round.(Int, (2π, π) ./ res .+ (0, 1)))  # CAR has pixels on poles
     end
@@ -20,8 +20,8 @@ function fullsky_geometry(proj::Type{<:CarProjection}, res; shape = nothing, dim
     return (nx, ny, dims...), wcs
 end
 
-fullsky_geometry(proj::Type{<:CarProjection}, res::Number; shape = nothing, dims = ()) =
+fullsky_geometry(proj::Type{<:CarClenshawCurtis}, res::Number; shape = nothing, dims = ()) =
     fullsky_geometry(proj, (res, res); shape = shape, dims = dims)
 
 fullsky_geometry(res; shape = nothing, dims = ()) =
-    fullsky_geometry(CarProjection, res; shape = shape, dims = dims)
+    fullsky_geometry(CarClenshawCurtis, res; shape = shape, dims = dims)

From 46f85d0f663ee7f7740a23e06ce6c3caf079ffde Mon Sep 17 00:00:00 2001
From: xzackli <xzackli@gmail.com>
Date: Mon, 24 Jan 2022 02:33:08 -0500
Subject: [PATCH 3/5] add full sky geometry tests

---
 test/runtests.jl | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index fd1d8d2..eccca28 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,18 @@
 using Pixell
 using Test
 
-@testset "Pixell.jl" begin
-    # Write your tests here.
+
+@testset "Enmap geometry" begin
+    shape, wcs = fullsky_geometry(deg2rad(1 / 60))
+    @test wcs.cdelt ≈ [-0.016666666666666666, 0.016666666666666666]
+    @test wcs.crpix ≈ [10800.5, 5401.0]
+    @test wcs.crval ≈ [0.008333333333333333, 0.0]
+
+    shape, wcs = fullsky_geometry(deg2rad(1 / 61))
+    @test wcs.cdelt ≈ [-0.01639344262295082, 0.01639344262295082]
+    @test wcs.crpix ≈ [10980.5, 5491.0]
+    @test wcs.crval ≈ [0.00819672131147541, 0.0]
+
+    shape, wcs = fullsky_geometry(deg2rad(5); dims=(3,))
+    @test shape == (72, 37, 3)
 end

From 3487d0519fa73b7a2850adcaf11cbf518f4927ed Mon Sep 17 00:00:00 2001
From: xzackli <xzackli@gmail.com>
Date: Mon, 24 Jan 2022 02:47:04 -0500
Subject: [PATCH 4/5] add broadcasting testset and also fix bug with
 broadcasting

---
 src/enmap.jl     |  2 +-
 test/runtests.jl | 11 +++++++++++
 2 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/src/enmap.jl b/src/enmap.jl
index 97f84ef..f8156ba 100644
--- a/src/enmap.jl
+++ b/src/enmap.jl
@@ -88,7 +88,7 @@ struct NoWCS end
 combine(x::NoWCS, y) = y
 combine(x, y::NoWCS) = x
 combine(x::NoWCS, ::NoWCS) = x
-combine(x::Enmap, y::Enmap) = x  # TODO: check compatibility
+combine(x::WCSTransform, y::WCSTransform) = x  # TODO: check compatibility
 
 
 ####
diff --git a/test/runtests.jl b/test/runtests.jl
index eccca28..9c485e5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,3 +16,14 @@ using Test
     shape, wcs = fullsky_geometry(deg2rad(5); dims=(3,))
     @test shape == (72, 37, 3)
 end
+
+
+@testset "Enmap broadcasting" begin
+    shape, wcs = fullsky_geometry(deg2rad(1); dims=(3,))
+    A, B = rand(shape...), rand(shape...)
+    ma = Enmap(A, wcs)
+    mb = Enmap(B, wcs)
+    @test A .+ B == ma .+ mb
+    @test A .+ B == ma .+ B
+    @test A .+ B .* sin.(A.^2) == (ma .+ mb .* sin.(ma.^2))
+end

From e11badd5cf9efbecfe321ebfcc9d0e77f83bd884 Mon Sep 17 00:00:00 2001
From: xzackli <xzackli@gmail.com>
Date: Mon, 24 Jan 2022 03:06:32 -0500
Subject: [PATCH 5/5] add docstring to fullsky_geometry

---
 src/enmap_ops.jl | 24 +++++++++++++++++++++++-
 test/runtests.jl |  1 +
 2 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/src/enmap_ops.jl b/src/enmap_ops.jl
index c10e19c..65199ee 100644
--- a/src/enmap_ops.jl
+++ b/src/enmap_ops.jl
@@ -1,5 +1,27 @@
 
-function fullsky_geometry(proj::Type{<:CarClenshawCurtis}, res; shape = nothing, dims = ())
+"""
+    fullsky_geometry([P=CarClenshawCurtis], res; shape = nothing, dims = ())
+
+# Arguments:
+- `proj=CarClenshawCurtis`: [optional] projection
+- `res`: resolution in radians. Passing a Number produces a square pixel.
+    Passing a tuple with (ΔRA, ΔDEC) produces a rectangular pixel.
+
+# Keywords
+- `shape::NTuple=nothing`: shape of the map. If not specified, will be computed.
+- `dims::NTuple=()`: additional dimensions to append to the shape, such as (3,) for IQU
+    to generate a map with `(nx, ny, 3)`.
+
+# Returns: 
+- `shape::Tuple, wcs::WCSTransform`: a tuple containing the shape of the map and the WCS
+
+# Examples
+```julia-repl
+julia> shape, wcs = fullsky_geometry(deg2rad(30/60))  # 30 arcmin pixel
+((720, 361), WCSTransform(naxis=2,cdelt=[-0.5, 0.5],crval=[0.25, 0.0],crpix=[360.5, 181.0]))
+```
+"""
+function fullsky_geometry(P::Type{<:CarClenshawCurtis}, res; shape = nothing, dims = ())
     if isnothing(shape)
         shape = (round.(Int, (2π, π) ./ res .+ (0, 1)))  # CAR has pixels on poles
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9c485e5..35bbfad 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,5 +25,6 @@ end
     mb = Enmap(B, wcs)
     @test A .+ B == ma .+ mb
     @test A .+ B == ma .+ B
+    @test A .+ B == A .+ mb
     @test A .+ B .* sin.(A.^2) == (ma .+ mb .* sin.(ma.^2))
 end