Handle white noise with healpy 2 partial sky coverage

Simulate white noise maps and use hitmaps
cosmology
python
healpy
Published

June 20, 2020

import healpy as hp
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import astropy.units as u
hp.disable_warnings()

In this second notebook, we will handle a case of partial sky coverage.

# Number based on Simons Observatory SAT UHF1 array of detectors

net = 10. * u.Unit("uK * sqrt(s)")

5 years with a efficiency of 20%:

integration_time_total = 5 * u.year * .2

Download a hitmap

We can download a simulated hitmap for a Simons Observatory band, for now however, we assume a uniform coverage.

hitmap_url = "https://portal.nersc.gov/project/sobs/so_mapsims_data/v0.2/healpix/ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz"
!wget $hitmap_url
--2020-12-11 10:32:00--  https://portal.nersc.gov/project/sobs/so_mapsims_data/v0.2/healpix/ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz
Resolving portal.nersc.gov (portal.nersc.gov)... 128.55.206.24, 128.55.206.26
Connecting to portal.nersc.gov (portal.nersc.gov)|128.55.206.24|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 3212515 (3.1M) [application/x-gzip]
Saving to: ‘ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz.8’

ST0_UHF1_01_of_20.n 100%[===================>]   3.06M  3.12MB/s    in 1.0s    

2020-12-11 10:32:06 (3.12 MB/s) - ‘ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz.8’ saved [3212515/3212515]

hitmap = hp.read_map("ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz")
hitmap = hitmap > 0
hp.mollview(hitmap)

Uniform partial sky survey

As a reference, let’s first start with the trivial case of uniform full sky coverage, i.e. we spend the same amount of observation time in each pixel.

nside = 512
npix = hp.nside2npix(nside)
npix_hit = hitmap.sum()
standard_deviation_per_pixel = (net / np.sqrt(integration_time_total/npix_hit)).decompose()
standard_deviation_per_pixel

\(1.9596849 \times 10^{-6} \; \mathrm{K}\)

m = np.nan * np.ones(npix, dtype=np.double) * standard_deviation_per_pixel.unit
m[hitmap] = np.random.normal(scale = standard_deviation_per_pixel.value, size=npix_hit) * standard_deviation_per_pixel.unit
m = m.to(u.uK)
m.value[np.isnan(m)] = hp.UNSEEN
hp.mollview(m, unit=m.unit, title="White noise map")

Power spectrum

sky_fraction = np.mean(hitmap)
cl = hp.anafast(m) / sky_fraction
cl[100:].mean()
1.5205726153293572e-05
m[hitmap].std()

\(1.9595938 \; \mathrm{\mu K}\)

pixel_area = hp.nside2pixarea(nside)
white_noise_cl = (standard_deviation_per_pixel**2 * pixel_area).to(u.uK**2)
white_noise_cl

\(1.5341266 \times 10^{-5} \; \mathrm{\mu K^{2}}\)

white_noise_cl_full_sky = 3.9820426e-5 * u.uK**2
plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.hlines(white_noise_cl.value, 0, len(cl), color="blue",
           label="White noise level")
plt.hlines(white_noise_cl_full_sky.value, 0, len(cl),
           label="White noise level full sky")
plt.title("Fullsky white noise spectrum")
plt.legend()
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$");

white_noise_cl / white_noise_cl_full_sky

\(0.38526121 \; \mathrm{}\)

sky_fraction
0.38526121775309247