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

Noise-Equivalent-Temperature, it is a measure of sensitivity of a detector, in cosmology, it is often quoted in $\mu K \sqrt(s)$, i.e. it is the sensitivity per unit time and can be divided by the integration time to get the actual standard deviation of the white noise of the instrument.

For example let's consider a white noise NET of $200 \mu K \sqrt(s)$

it means that if you integrate for 100 seconds for each pixel, the standard deviation will be $20 \mu K$.

net = 200 * u.Unit("uK * sqrt(s)")
net
$200 \; \mathrm{\mu K\,s^{1/2}}$
integration_time_per_pixel = 100 * u.s
standard_deviation = net / np.sqrt(integration_time_per_pixel)
standard_deviation = 2 * u.uK

## Create a white noise map

Now that we have an estimate of the standard deviation per pixel, we can use numpy to create a map of gaussian white noise.

nside = 128
npix = hp.nside2npix(nside)
m = np.random.normal(scale = standard_deviation.value, size=npix) * standard_deviation.unit
hp.mollview(m, unit=m.unit, title="White noise map")

## Power spectrum

Finally we can compute the angular power spectrum with anafast, i.e. the power as a function of the angular scales, from low $\ell$ values for large angular scales, to high $\ell$ values for small angular scales.

At low $\ell$ there is not much statistics and the power spectrum is biased, but if we exclude lower ells, we can have an estimate of the white noise $C_\ell$ coefficients. We can then compare with the theoretical power computed as:

$$C_\ell = \Omega_{pix}\sigma$$

Where: $\Omega_{pix}$ is the pixel are in square-radians and $\sigma$ is the white noise standard deviation.

cl = hp.anafast(m)
cl[100:].mean()
0.00025151565299227235
pixel_area = hp.nside2pixarea(nside)
white_noise_cl = standard_deviation.value**2 * pixel_area
white_noise_cl
0.00025566346464760684
plt.figure(figsize=(6,4))
plt.loglog(cl, label="Map power spectrum", alpha=.7)
plt.hlines(white_noise_cl, 0, len(cl), label="White noise level")
plt.xlabel("$\ell$")
plt.ylabel("$C_\ell [\mu K ^ 2]$");