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

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-06-20 12:51:21--  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.201.128
Connecting to portal.nersc.gov (portal.nersc.gov)|128.55.201.128|: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.3’

ST0_UHF1_01_of_20.n 100%[===================>]   3.06M  --.-KB/s    in 0.1s    

2020-06-20 12:51:21 (29.0 MB/s) - ‘ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz.3’ saved [3212515/3212515]

hitmap = hp.read_map("ST0_UHF1_01_of_20.nominal_telescope_all_time_all_hmap.fits.gz")
/home/zonca/zonca/p/software/healpy/healpyvenv/lib/python3.7/site-packages/healpy/fitsfunc.py:352: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
  "If you are not specifying the input dtype and using the default "
NSIDE = 512
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING
hitmap = hitmap > 0
hp.mollview(hitmap)

Uniform full 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

cl = hp.anafast(m) / sky_fraction
cl[100:].mean()
1.5184872877653936e-05
m[hitmap].std()
$1.9568977 \; \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.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$");
white_noise_cl / white_noise_cl_full_sky
$0.38526121 \; \mathrm{}$
sky_fraction
0.38526121775309247