add some tests from Pixell for scalar SHTs

Project.toml

@@ -26,6 +26,7 @@ julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
 [targets]
-test = ["Test"]
+test = ["Test", "DelimitedFiles"]

src/transforms.jl

@@ -57,7 +57,7 @@ function make_cc_geom_info(shape, wcs₀::AbstractWCSTransform)
         Cint(nsubrings),nph,  offsets,           stride,    phi0s,         theta,        weights,      geom_info_ptr
     )
 
-    GeomInfo(geom_info_ptr[])
+    return GeomInfo(geom_info_ptr[])
 end
 
 
@@ -87,12 +87,12 @@ function map2alm(input_map::Enmap{T,2}; lmax=nothing, mmax=lmax) where T
     nalms = alm_count(alm_info)
     npix = map_size(geom_info)
     alm = Alm(lmax, mmax, zeros(ComplexF64, nalms))
-    map1d = reshape(band.data, npix)
+    flat_map = reshape(band.data, npix)
 
     sharp_execute!(
         SHARP_MAP2ALM, 0, 
         [alm.alm], 
-        [map1d],
+        [flat_map],
         geom_info, alm_info, SHARP_DP
     )
     return alm

test/gen_sht_test_data.py

@@ -0,0 +1,26 @@
+# %%
+from pixell import enmap, utils, curvedsky
+import numpy as np
+
+# add some simple analytic stuff to the map
+def gen_simple(imap, powerlaw_exponent):
+    for i in range(imap.shape[0]):
+        for j in range(imap.shape[1]):
+            imap[i,j] = (i * imap.shape[1] + j + 1)**powerlaw_exponent
+
+shape, wcs = enmap.fullsky_geometry(10.0*utils.degree)
+imap = enmap.zeros(shape, wcs=wcs)
+gen_simple(imap, 2)
+
+alms = curvedsky.map2alm(imap, lmax=18, method="auto")
+np.savetxt("data/simple_analytic_sht.txt", np.column_stack([alms.real, alms.imag]), fmt='%.60g')
+
+alms = curvedsky.map2alm(imap[4:-3, 5:-2], lmax=18, method="auto")
+np.savetxt("data/simple_analytic_sht_sliced.txt", np.column_stack([alms.real, alms.imag]), fmt='%.60g')
+
+box = np.array([[-5,10],[5,-10]]) * utils.degree
+shape, wcs = enmap.geometry(pos=box,res=60 * utils.arcmin,proj='car')
+imap = enmap.zeros(shape, wcs=wcs)
+gen_simple(imap, 2.5)
+alms = curvedsky.map2alm(imap, lmax=100, method="auto")
+np.savetxt("data/simple_box_analytic_sht.txt", np.column_stack([alms.real, alms.imag]), fmt='%.60g')

test/runtests.jl

@@ -2,6 +2,9 @@ using Pixell
 using Test
 import Pixell: degree, arcminute
 
+using DelimitedFiles
+
 include("test_geometry.jl")  # creating geometries and sky ↔ pix
 include("test_enmap.jl")     # enmap features and manipulation
+include("test_transforms.jl")
 include("test_io.jl")

test/test_transforms.jl

@@ -0,0 +1,36 @@
+
+# write some simple stuff to the map
+function gen_simple!(m, powerlaw_exponent)
+    for i in axes(m, 1)
+        for j in axes(m, 2)
+            m[i,j] = ( ((j-1) * size(m, 1)) + i )^powerlaw_exponent
+        end
+    end
+end
+
+@testset "Scalar SHTs" begin
+    shape, wcs = fullsky_geometry(10.0 * Pixell.degree)
+    m = Enmap(zeros(shape), wcs)
+    gen_simple!(m, 2)
+
+    alms = map2alm(m; lmax=18)
+    data = readdlm("data/simple_analytic_sht.txt")
+    ref_alms = data[:,1] + 1im .* data[:,2]
+    @test (alms.alm ≈ ref_alms)
+
+    alms = map2alm(m[6:end-2, 5:end-3]; lmax=18)
+    data = readdlm("data/simple_analytic_sht_sliced.txt")
+    ref_alms = data[:,1] + 1im .* data[:,2]
+    @test (alms.alm ≈ ref_alms)
+
+
+    box = [10   -10;           # RA
+           -5     5] * degree  # DEC
+    shape, wcs = geometry(CarClenshawCurtis, box, 1.0 * degree)
+    m = Enmap(zeros(shape), wcs)
+    gen_simple!(m, 2.5)
+    alms = map2alm(m; lmax=100)
+    data = readdlm("data/simple_box_analytic_sht.txt")
+    ref_alms = data[:,1] + 1im .* data[:,2]
+    @test (alms.alm ≈ ref_alms)
+end