laxes and extent_cyl

Project.toml

@@ -1,7 +1,7 @@
 name = "Pixell"
 uuid = "d4cbf8b2-a14f-4080-9f5e-d5a64ada768e"
 authors = ["Zack Li", "Yilun Guan"]
-version = "0.2.8"
+version = "0.2.9"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/Pixell.jl

@@ -40,6 +40,7 @@ export Alm, map2alm, alm2cl, alm2map, alm2map!
 export distance_transform, ApproxSeqSDT, ExactSeqSDT, BruteForceSDT
 export dplanck, RadialFourierTransform, real2harm, harm2real
 export FFTLogPlan
+export extent_cyl, laxes_cyl
 
 # set up some shortcuts for common angles
 const radian = Unitful.rad

src/projections/arbitrary_wcs.jl

@@ -131,6 +131,36 @@ function skyarea_cyl(shape, wcs::AbstractWCSTransform)
     return (sin(δ₂) - sin(δ₁)) * abs(Δα) * N_α
 end
 
+function extent_cyl(shape, wcs::AbstractWCSTransform; signed=false)
+    N_α, N_δ = shape
+    α, δ = pix2sky(shape, wcs, SVector(0., 0.), SVector(0.5, N_δ+0.5); safe=false)
+    δ₁, δ₂ = sort(δ)
+    δ₁, δ₂ = max(-π/2, δ₁), min(π/2, δ₂)
+    δsign = δ₁ ≤ δ₂ ? 1 : -1
+    Δα, Δδ = getcdelt(wcs) .* getunit(wcs)
+    mean_cos   = (sin(δ₂) - sin(δ₁)) / (δ₂ - δ₁)
+    ext = (N_α * Δα * mean_cos, (δ₂ - δ₁)*δsign)
+    if signed
+        return ext
+    else
+        return abs.(ext)
+    end
+end
+
+
+"""
+    laxes_cyl(shape, wcs)
+
+Get multipole axes for a cylindrical pixelization, with a compromise choice for handling the
+spherical geometry that changes the angle between meridians in different declinations.
+"""
+function laxes_cyl(shape, wcs)
+    Δᾱ, Δδ̄ = extent_cyl(shape[1:2], wcs; signed=true) ./ shape
+    ℓ_δ = fftfreq(shape[2], 2π / Δδ̄)
+    ℓ_α = fftfreq(shape[1], 2π / Δᾱ)
+    return ℓ_α, ℓ_δ
+end
+
 function sliced_wcs(wcs::WCSTransform, cdelt′, crpix′)
     new_wcs = deepcopy(wcs)
     new_wcs.cdelt = collect(cdelt′)

test/test_geometry.jl

@@ -314,3 +314,28 @@ end
         @test sum(abs.(pm[1,:] .- pm_py)) < 100eps()
     end
 end
+
+@testset "extent_cyl" begin
+    shape = (3612, 1605)
+    wcs = Pixell.create_car_wcs(CarClenshawCurtis{Float64},  
+        (-0.00833333333333, 0.00833333333333), # cdelt
+        (1806.0, 1358.0),                      # crpix
+        (33.9416666667, 0.0)                   # crval
+    )
+    ref = ([0.5224453478223857, 0.23343778745414823])
+    ext = Pixell.extent_cyl(shape, wcs)
+    @test ref ≈ collect(ext)
+
+    ℓ_α, ℓ_δ = laxes_cyl(shape, wcs)
+    @test ℓ_α[0+1] ≈ -0.0
+    @test ℓ_α[720+1] ≈ -8659.074944442375
+    @test ℓ_α[1440+1] ≈ -17318.14988888475
+    @test ℓ_α[2160+1] ≈ 17462.467804625456
+    @test ℓ_α[2880+1] ≈ 8803.392860183081
+    @test ℓ_α[3600+1] ≈ 144.31791574070627
+
+    @test ℓ_δ[0+1] ≈ 0.0
+    @test ℓ_δ[400+1] ≈ 10766.355140191221
+    @test ℓ_δ[1200+1] ≈ -10900.934579443612
+    @test ℓ_δ[1600+1] ≈ -134.57943925239027
+end