<jupyter_start><jupyter_text>Tutorial This notebook demonstrates the usage of `redshifted_gaussian_fields` to generate realizations of an intensity field on the sky generated by redshifted 21cm emission with a chosen power spectrum. Table of Contents1&nbsp;&nbsp;Define Power Spectrum2&nbsp;&nbsp;Compute Cross-frequency Covariance3&nbsp;&nbsp;Generate Random Realizations4&nbsp;&nbsp;Save and Restore Completed Calculations<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from astropy import cosmology
from redshifted_gaussian_fields import generator
plt.rcParams['figure.figsize'] = [10,8]

plt.rcParams['xtick.labelsize'] = 'large'
plt.rcParams['ytick.labelsize'] = 'large'
def plot_power_spectrum(k_axis, Pspec):
    plt.plot(k_axis, Pspec(k_axis), label=r'$P(k)$')
    plt.plot(k_axis, (k_axis/0.2)**(-2.7), label=r'$ \propto k^{-2.7}$')

    plt.legend(loc='upper right', fontsize=16)
    plt.yscale('log')
    plt.xscale('log')
    plt.ylim(1e-6,1e4)
    plt.xlabel(r'$k$ (Mpc$^{-1}$)',fontsize=20)
    plt.ylabel(r'$P(k)$ (K$^2$)',fontsize=20)
    
    plt.show()<jupyter_output><empty_output><jupyter_text>Define Power Spectrum Begin by defining the power spectrum. In this example we will choose a power spectrum that is effectively a power law over the given range.<jupyter_code>k0 = np.logspace(-2.,1.,11)
a = k0**(-2.7)

normalization_point = 0.2
normalization_amplitude = 1.

Pspec = generator.ParameterizedGaussianPowerSpectrum(a, k0, renormalization=(normalization_point, normalization_amplitude), term_type='flat_gauss')
k_axis = np.logspace(-2.5,1.5,1001)
plot_power_spectrum(k_axis, Pspec)<jupyter_output><empty_output><jupyter_text>The components used to define the power spectrum do not span all the possible power spectra, but they are flexible enough to produce a variety of forms. For example, we might add a run on the spectral index:<jupyter_code># Using the same locations parameters,but with a different set
# of amplitude parameters
a_2 = k0**(-2.7 - 0.1*np.log(k0))

Pspec_2 = generator.ParameterizedGaussianPowerSpectrum(a_2, k0, renormalization=(normalization_point, normalization_amplitude), term_type='flat_gauss')
plot_power_spectrum(k_axis, Pspec_2)<jupyter_output><empty_output><jupyter_text>Compute Cross-frequency Covariance Next we choose the parameters over which to evaluate the cross-covariance function which is used to produce sky model realizations.Frequency $\nu$ points in MHz, 140MHz to 150MHz with a channel width of 0.1MHz:<jupyter_code>nu_axis = np.linspace(140., 150., 101, endpoint=True)
del_nu = 0.1<jupyter_output><empty_output><jupyter_text>Angular harmonic order $\ell$, from 0 to 100<jupyter_code>ell_axis = np.arange(101)<jupyter_output><empty_output><jupyter_text>A choice of cosmological parameters is needed to specify the comoving distance $r$ as a function of $\nu$. This package uses an `astropy` cosmology object for this purpose:<jupyter_code>cosmo = cosmology.Planck15<jupyter_output><empty_output><jupyter_text>Choose the number of samples used to compute the frequency channel integrals<jupyter_code>Np = 15<jupyter_output><empty_output><jupyter_text>Now we can initialize a GaussianCosmologicalFieldGenerator object, which is used to manage the various computations and well as saving to disk and restoring completed calculations.<jupyter_code>gcfg = generator.GaussianCosmologicalFieldGenerator(cosmo, Pspec, nu_axis, del_nu, ell_axis, Np=Np)<jupyter_output><empty_output><jupyter_text>Then we can compute the cross-frequency angular power spectrum over the specified parameter domain for the chosen 3D power spectrum `Pspec`. This is most computationally intensive step and for useful ranges of the parameters may take minutes to hours depending on the number of frequency samples, angular modes, number of terms in the power spectrum function, number of points used in the channelization integral, and whether a matrix element cutoff value $\epsilon$ can be effectively utilized.<jupyter_code>%%time
gcfg.compute_cross_frequency_angular_power_spectrum()<jupyter_output>CPU times: user 2min 42s, sys: 3.09 s, total: 2min 45s
Wall time: 26 s<jupyter_text>The computational scaling in terms of* $N_\nu$, the number of frequency samples* $N_\ell$, the number of angular modes* $N_i$, the number of power spectrum terms* $N_p$, the number of samples per frequency channelis approximately $\mathcal{O}(N_\nu^2 N_\ell N_i N_p^2)$ if $\epsilon = 0$. However if $\epsilon > 0$, it becomes $\sim \mathcal{O}(N_\nu N_\ell N_i N_p^2)$, with an overall nearly-constant factor that depends on $\epsilon$ and the small-$k_0$ extent of the power spectrum. The computed cross-frequency angular power spectrum is stored in a $N_\ell \times N_\nu \times N_\nu$ array<jupyter_code>print(gcfg.barC.shape)
plt.imshow(gcfg.barC[50]);
plt.plot(gcfg.nu_axis, gcfg.barC[50,0,:], '.')
plt.plot(gcfg.nu_axis, gcfg.barC[50,50,:], '.')
plt.plot(gcfg.nu_axis, gcfg.barC[50,75,:], '.');<jupyter_output><empty_output><jupyter_text>Generate Random Realizations Generate realizations of the harmonic coefficients of the specific intensity of the sky. Each realization is generated from the input random seed so that the same realization can be re-generated on demand.<jupyter_code>seed = 2983
a_lm = gcfg.generate_realization(seed)
a_lm2 = gcfg.generate_realization(seed)
np.allclose(a_lm, a_lm2)<jupyter_output><empty_output><jupyter_text>Generate a healpix map sampling of the specific intensity with resolution specified by `nside`:<jupyter_code>nside = 128
hmap = gcfg.generate_healpix_map_realization(seed, nside)<jupyter_output>/lustre/aoc/projects/hera/zmartino/anaconda3/envs/main/lib/python3.7/site-packages/healpy/sphtfunc.py:824: UserWarning: Sigma is 0.000000 arcmin (0.000000 rad) 
  sigma * 60 * 180 / np.pi, sigma
/lustre/aoc/projects/hera/zmartino/anaconda3/envs/main/lib/python3.7/site-packages/healpy/sphtfunc.py:829: UserWarning: -> fwhm is 0.000000 arcmin
  sigma * 60 * 180 / np.pi * (2.0 * np.sqrt(2.0 * np.log(2.0)))<jupyter_text>The first axis of the array is frequnency, the second axis is healpix pixels.<jupyter_code>print(hmap.shape, hmap.dtype)<jupyter_output>(101, 196608) float64<jupyter_text>The frequency spectrum of at a point in the map:<jupyter_code>plt.plot(gcfg.nu_axis, hmap[:,0], '.-')
plt.ylabel('Jy / sr')
plt.xlabel('MHz');<jupyter_output><empty_output>