add a unit test for tsz

Artifacts.toml

@@ -6,3 +6,13 @@ lazy = true
     [[tsz_example.download]]
     url = "https://github.com/WebSky-CITA/XGPaint.jl/releases/download/tsz_example/tsz_example.tar.gz"
     sha256 = "33242449b9a67e9c21e243e80c61bc00123ffc3325cd109aa1f045bf7d798e50"
+
+
+[tsz_test_map]
+git-tree-sha1 = "e3d3d9f9fa5ea25feb5d8dacdd8334867f5db287"
+lazy = true
+
+    [[tsz_test_map.download]]
+    url = "https://github.com/WebSky-CITA/XGPaint.jl/releases/download/tsz_example/test_tsz_map.tar.gz"
+    sha256 = "5142052ac15e9fe0c169ad568439a07e29e274a6b0ed33cf46304413a7117de7"
+

src/util.jl

@@ -2,7 +2,7 @@ using HDF5
 using Healpix
 using Random
 using Random: MersenneTwister
-using Pkg.Artifacts
+using LazyArtifacts
 
 """
 Utility function to read an HDF5 table with x, y, z, M_h as the four rows.
@@ -31,6 +31,12 @@ function load_example_halos()
     return ra, dec, redshift, halo_mass
 end
 
+function load_example_tsz_map()
+    rootpath = artifact"tsz_test_map"
+    ymap = jldopen(joinpath(rootpath, "test_tsz_map.jld2"), "r")["data"]
+    return ymap
+end
+
 
 """
 Reads in a standard Battaglia 2016 model from disk for tSZ.

test/Project.toml

@@ -5,6 +5,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Healpix = "9f4e344d-96bc-545a-84a3-ae6b9e1b672b"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
+Pixell = "d4cbf8b2-a14f-4080-9f5e-d5a64ada768e"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -2,6 +2,7 @@ using XGPaint
 using Cosmology
 using Test
 using HDF5
+using Pixell
 
 println("")  # useful in Atom console
 
@@ -53,7 +54,7 @@ cosmo = XGPaint.get_cosmology(h=0.7f0, OmegaM=0.25f0)
 
     # python: h2fm.fluxmodel.sigma_cen(1e13)
     @test isapprox(XGPaint.sigma_cen(1.0f13, model),
-        3218663770378.3066, rtol=rtol)
+        3218663770378.3066, rtol=rtoXGPaint.Pixelll)
 
     # python: h2fm.fluxmodel.nu2theta(150e9, 0.5)
     @test isapprox(XGPaint.nu2theta(150f9, 0.5f0, model),
@@ -80,9 +81,7 @@ cosmo = XGPaint.get_cosmology(h=0.7f0, OmegaM=0.25f0)
     @test muofn(500.0f0) ≈ 6.14975653e-05
     @test muofn(1.0f0) ≈ 0.11765558
 
-    @test XGPaint.z_evo(0.0f0, model) ≈ 1.0f0
-end
-
+    @test XGPaint.z_evo(0.0f0, modelXGPaint.Pixell
 
 radio_model = Radio_Sehgal2009{Float32}()
 @testset "radio" begin
@@ -114,8 +113,25 @@ end
     @test all(XGPaint.chunk(10, 3) == [(1,3), (4,6), (7,9), (10,10)])
     @test all(XGPaint.chunk(10, 4) == [(1,4), (5,8), (9,10)])
     @test all(XGPaint.chunk(10, 10) == [(1,10)])
+end
 
+##
+@testset "tsz" begin
+    ra, dec, redshift, halo_mass = XGPaint.load_example_halos()
+    ra, dec, redshift, halo_mass = sort_halo_catalog(ra, dec, redshift, halo_mass);
+    model, interp = XGPaint.load_precomputed_battaglia()
+    box = [4.5   -4.5;           # RA
+       -3     3] * Pixell.degree  # DEC
+    shape, wcs = geometry(Pixell.CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute)
+    m = Enmap(zeros(shape), wcs)
+    workspace = profileworkspace(shape, wcs)
+    @time paint!(m, model, workspace, interp, halo_mass, redshift, ra, dec)
+
+    ymap_ref = XGPaint.load_example_tsz_map()
+    @test sum(abs, m.data .- ymap_ref) ≈ 0 atol=1e-7
 end
 
+
+
 ##
 include("test_query.jl")