Matter distribution

This example simulates only the matter field in nested shells up to redshift 1.

Setup

Set up a matter-only GLASS simulation, which requires angular matter power spectra and the sampling itself (here: lognormal).

The setup for angular matter power spectra matches the definition from the Matter shell definition example.

import healpy as hp
import matplotlib.pyplot as plt
import numpy as np

# uses the CAMB cosmology which produced the cls
import camb
from cosmology import Cosmology

# almost all GLASS functionality is available from the `glass` namespace
import glass
import glass.ext.camb

# creating a numpy random number generator for sampling
rng = np.random.default_rng(seed=42)

# cosmology for the simulation
h = 0.7
Oc = 0.25
Ob = 0.05

# basic parameters of the simulation
nside = 1024
lmax = 1000

# set up CAMB parameters for matter angular power spectrum
pars = camb.set_params(
    H0=100 * h,
    omch2=Oc * h**2,
    ombh2=Ob * h**2,
    NonLinear=camb.model.NonLinear_both,
)

# get the cosmology from CAMB
cosmo = Cosmology.from_camb(pars)

# shells of 200 Mpc in comoving distance spacing
zb = glass.distance_grid(cosmo, 0.0, 1.0, dx=200.0)

# triangular radial windows, equivalent to linear interpolation of n(z)
shells = glass.linear_windows(zb)

# compute angular matter power spectra with CAMB
cls = glass.ext.camb.matter_cls(pars, lmax, shells)

# apply discretisation to the full set of spectra:
# - HEALPix pixel window function (`nside=nside`)
# - maximum angular mode number (`lmax=lmax`)
# - number of correlated shells (`ncorr=3`)
cls = glass.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)

# compute Gaussian spectra for lognormal fields from discretised spectra
gls = glass.lognormal_gls(cls)

# this generator will yield the matter fields in each shell
matter = glass.generate_lognormal(gls, nside, ncorr=3, rng=rng)

Simulation

Run the simulation. For each shell, plot an orthographic annulus of the matter distribution.

# make a 2d grid in redshift
n = 2000
zend = 1.05 * zb[-1]
x, y = np.mgrid[-zend : zend : 1j * n, -zend : zend : 1j * n]
z = np.hypot(x, y)
grid = np.full(z.shape, np.nan)

# set up the plot
ax = plt.subplot(111)
ax.axis("off")

# simulate and project an annulus of each matter shell onto the grid
for i, delta_i in enumerate(matter):
    zmin, zmax = zb[i], zb[i + 1]
    g = (zmin <= z) & (z < zmax)
    zg = np.sqrt(1 - (z[g] / zmax) ** 2)
    theta, phi = hp.vec2ang(np.transpose([x[g] / zmax, y[g] / zmax, zg]))
    grid[g] = hp.get_interp_val(delta_i, theta, phi)
    ax.add_patch(
        plt.Circle((0, 0), zmax / zend, fc="none", ec="k", lw=0.5, alpha=0.5, zorder=1),
    )

# show the grid of shells
ax.imshow(grid, extent=[-1, 1, -1, 1], zorder=0, cmap="bwr", vmin=-2, vmax=2)

# show the resulting plot
plt.show()