The Kinetic Sunyaev–Zeldovich Effect

Let's make some maps of the kinetic Sunyaev-Zeldovich (SZ) effect. First, let's load up the example halo catalog included in this package. These will be automatically downloaded the first time you load them.

using XGPaint, Plots, Unitful, UnitfulAstro

# example , ra and dec in radians, halo mass in M200c (Msun)
ra, dec, redshift, halo_mass = XGPaint.load_example_halos()
print("Number of halos: ", length(halo_mass))

# the current test halos do not include velocities, so we randomly generate some for this
proj_v_over_c = randn(eltype(ra), length(halo_mass)) / 1000
2455621-element Vector{Float32}:
 -0.0003219509
 -0.000118726006
  0.0017187106
  0.00044714267
  0.00064462767
  0.0010979729
  0.00054252613
 -0.0002649792
 -0.0012389072
 -4.4737648f-5
  ⋮
 -0.0006623338
  0.00080121553
  0.00046369727
 -0.0013890561
 -3.6066253f-5
  0.0013778644
 -0.0014373824
  0.00044164242
 -0.0012948895

This small catalog is limited to a relatively small patch of the sky. Before we generate some SZ maps, let's take a look at the halo mass distribution.

b = 10.0 .^ (11:0.25:16)
histogram(halo_mass, bins=b, xaxis=(:log10, (1e11, 1e16)),
    yscale=:log10, label="", xlabel="Halo mass (solar masses)", ylabel="counts")
Example block output

To allow for safe threaded painting on a single map, we'll sort the halo catalog by declinations.

ra, dec, redshift, halo_mass = sort_halo_catalog(ra, dec, redshift, halo_mass);
(Float32[0.03305517, -0.052936707, 0.07140066, -0.06599737, -0.05942857, -0.05192846, 0.022953046, 0.01955623, 0.05080256, 0.00839793  …  0.060685053, -0.061119348, -0.040161338, 0.047465798, -0.041020952, -0.048511826, -0.018889412, 0.058788735, -0.050192643, -0.0027171609], Float32[-0.11871151, -0.118711114, -0.118711084, -0.11871107, -0.1187108, -0.118710734, -0.11871068, -0.118710674, -0.11871066, -0.11871049  …  0.10907136, 0.1090714, 0.10907147, 0.109071486, 0.1090715, 0.109071545, 0.10907169, 0.10907171, 0.10907183, 0.10907187], Float32[2.583177, 0.6798465, 1.4767964, 2.4593022, 1.3147808, 1.7026237, 1.313409, 1.2228533, 2.061896, 1.7019778  …  0.4401552, 3.107578, 1.2217613, 2.3513637, 2.2199023, 1.7342392, 0.6818484, 0.8421148, 1.2661853, 2.0941606], Float32[1.8216356f12, 3.4135917f12, 4.958732f12, 3.3517478f12, 2.6270643f12, 2.1461618f12, 1.5084386f12, 4.1571304f12, 1.3568221f12, 2.2222684f12  …  4.438609f12, 1.3447792f12, 1.5955485f12, 1.7920556f12, 8.5669514f12, 2.6231164f12, 2.2672384f12, 2.409821f13, 5.5713213f12, 2.7740707f12])

Now, we'll set up the map to generate. We will construct a small CAR (Clenshaw-Curtis variant) patch on the sky. You have to construct a workspace for each new sky patch.

using Pixell
box = [4.5   -4.5;           # RA
       -3     3] * Pixell.degree  # DEC
shape, wcs = geometry(CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute)
m = Enmap(zeros(shape), wcs)
1080×720 Pixell.Enmap{Float64, 2, Matrix{Float64}, Pixell.CarClenshawCurtis{Float64}}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮              ⋱            ⋮                   
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

Now let's set up an electron profile.

model = BattagliaTauProfile(Omega_c=0.267, Omega_b=0.0493,  h=0.6712)
BattagliaTauProfile{Float64, Cosmology.FlatLCDM{Float64}}(0.1558646854252292, Cosmology.FlatLCDM{Float64}(0.6712, 0.6836071403864511, 0.3163, 9.285961354891755e-5))

We can now compute the integrated electron density, 

workspace = profileworkspace(shape, wcs)

# this only needs to be done once
# to DIY: interp = build_interpolator(model, cache_file="cached_btau.jld2", overwrite=true)
model_interp = XGPaint.load_precomputed_battaglia_tau()

m = Enmap(zeros(shape), wcs)
paint!(m, workspace, model_interp, halo_mass, redshift, ra, dec, proj_v_over_c)
plot(log10.(abs.(m)), c = :thermal)
Example block output

Healpix

To generate Healpix maps, you'll need the Healpix.jl package.

using Healpix

nside = 2048
m_hp = HealpixMap{Float64,RingOrder}(nside)
max_radius = deg2rad(5.0)  # maximum radius to consider for the profile
w = HealpixProfileWorkspace(nside, max_radius)

@time paint!(m_hp, w, model_interp, halo_mass, redshift, ra, dec, proj_v_over_c)
# Healpix.saveToFITS(m_hp, "!y.fits", typechar="D")  # to save, uncomment

plot(m_hp)
Example block output