galsci / pysm

PySM 3: Sky emission simulations for Cosmic Microwave Background experiments
https://pysm3.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
36 stars 23 forks source link

[patch] websky radio galaxy units #138

Closed xzackli closed 2 years ago

xzackli commented 2 years ago

This PR removes a single character. The WebSky radio galaxy maps were released as Jy/sr, while the CIB maps were released as MJy/sr. Due to copy-pasting the constructor kwargs, I had both components have input units of MJy/sr.

@zonca @giuspugl sorry about this bug, it evaded by various validation tests. I should have validated on the muK^2 output maps. This should resolve some of your radio source ringing issues, by suppressing them by a factor of a million.

I actually have generated maps smoothed in real space by a 0.9 arcminute Gaussian beam available on perlmutter, but I need to write up an interface somehow. Perhaps just a totally different component?

xzackli commented 2 years ago

I'd like to propose a sanity regression test for this. The middle CMB band is usually chosen to balance the power of the radio and CIB foreground components. Thus, it should have a similar power spectrum for the two components, and indeed this plot is how I found this bug. I think a very sensible unit test would be to simply compute the power spectra as below and then check that they are within i.e. two orders of magnitude.

import numpy as np
import matplotlib.pyplot as plt
import healpy as hp
from pysm3 import WebSkyCIB, WebSkyRadioGalaxies
import pysm3.units as u

model_cib = WebSkyCIB(nside=512)
m_cib = model_cib.get_emission(143*u.GHz)
cl_cib = hp.alm2cl(hp.map2alm(m_cib.value[0,:]))

model_radio = WebSkyRadioGalaxies(nside=512)
m_radio = model_radio.get_emission(143*u.GHz)
cl_radio = hp.alm2cl(hp.map2alm(m_radio.value[0,:]))

plt.plot(cl_cib, label="CIB")
plt.plot(cl_radio, label="Radio")
plt.legend()
plt.xscale("log"); plt.yscale("log")
plt.ylabel(r"$C_{\ell}$ [$\mu$K$^2]$ at 143 GHz")

image

Also here's the shot noise power for reference from Planck. image