The Relativistic Sunyaev–Zeldovich Effect
An example from Kuhn et al. (in prep). Begin with the tSZ example but on a smaller tile,
using Unitful, UnitfulAstro, XGPaint, Plots, Pixell
import PhysicalConstants.CODATA2018 as constants
M_200 = 1e15 * XGPaint.M_sun
z = 0.4
θ = 60*4.84814e-6
box = [-30 30; # RA
-30 30] * Pixell.arcminute # DEC
shape, wcs = geometry(Pixell.CarClenshawCurtis, box, 1 * Pixell.arcminute)
workspace = profileworkspace(shape, wcs)
model_tsz, tsz_interp = XGPaint.load_precomputed_battaglia();
tsz_snap = Enmap(zeros(shape), wcs)
mass_in_Msun = ustrip(M_200 / XGPaint.M_sun)
paint!(tsz_snap, model_tsz, workspace, tsz_interp, [mass_in_Msun], [z], [0.0], [0.0])
plot(log10.(tsz_snap))
tsz_snap = Enmap(zeros(shape), wcs)
mass_in_Msun = ustrip(M_200 / XGPaint.M_sun)
paint!(tsz_snap, model_tsz, workspace, tsz_interp, [mass_in_Msun], [z], [0.0], [0.0])
plot(log10.(tsz_snap))
i_ctr, j_ctr = round.(Int, sky2pix(tsz_snap, 0.0, 0.0))
ras = [pix2sky(tsz_snap, i, j_ctr)[1] for i in 1:size(tsz_snap,1)]
plot(ras, log10.(tsz_snap[:,j_ctr]), label="tSZ", xlabel="RA [rad]", ylabel="log10 y")
# generate a snapshot of a cluster at a specific frequency
function cluster_snapshot(nu, shape, wcs, model_tsz, tsz_interp, workspace)
X = nu_to_X(uconvert(u"s^-1",nu))
m = Enmap(zeros(shape), wcs);
# construct a new SZpack object for each frequency; this loads a table from disk.
# keep these around if you're looping over multiple maps at fixed frequency
model_szp = Battaglia16SZPackProfile(model_tsz, tsz_interp, X)
paint!(m, model_szp, workspace, [mass_in_Msun], [z], [0.0], [0.0])
return m
end
snap = cluster_snapshot(500u"GHz", shape, wcs, model_tsz, tsz_interp, workspace)
plot(log10.(snap))
freqs = LinRange(30.0, 800.0, 100) .* u"GHz"
snaps = [cluster_snapshot(nu, shape, wcs, model_tsz, tsz_interp, workspace) for nu in freqs]
I_rsz = [snap[i_ctr, j_ctr] * u"MJy/sr" for snap in snaps]
I_tsz = similar(I_rsz)
y0 = tsz_snap[i_ctr, j_ctr]
for (i, nu) in enumerate(freqs)
X = nu_to_X(nu)
I_tsz[i] = uconvert(u"MJy/sr",(X*(ℯ^X + 1)/(ℯ^X - 1) - 4) * y0 *
abs((2 * constants.h^2 * nu^4 * ℯ^X) / (
constants.k_B * constants.c_0^2 * XGPaint.T_cmb * (ℯ^X - 1)^2)))
end
plot(freqs, I_tsz, label="tSZ")
plot!(freqs, I_rsz, label="rSZ")