Running and Plotting LightCones

This tutorial follows on from the coeval cube tutorial, and provides an introduction to creating lightcones with 21cmFAST. If you are new to 21cmFAST you should go through the coeval cube tutorial first.

The main entry point into creating lightcones in 21cmFAST is the run_lightcone function. The function takes a few different arguments, most of which will be familiar to you if you’ve gone through the coeval tutorial. All simulation parameters can be passed (i.e. user_params, cosmo_params, flag_options and astro_params). As an alternative to the first two, an InitialConditions and/or PerturbField box can be passed.

Furthermore, the evolution can be managed with the zprime_step_factor and z_heat_max arguments.

There are in principle multiple ways to set up a lightcone: there is a choice of which line-of-sight slices are desired (eg. regular in frequency, comoving distance or redshift?), and also whether the box itself should be regular in transverse comoving distance (like the coeval boxes) or exist on an angular lattice. In 21cmFAST < 3.3.0, all lightcones were output in regular, comoving, rectilinear coordinates such that each “cell” was cubic and had the size of the coeval cubes used to make it. In 21cmFAST v3.3+, there is much more flexibility. You use a Lightconer subclass to tell the code what your lightcone should look like, and you can specify your own sub-class if the builtin ones don’t do what you require. The default lightconer is the same regular rectilinear grid as was used in 21cmFAST<3.3.

Note

Until v4, you will be able to call run_lightcone by directly passing a redshift and max_redshift (which was the only way to do it pre-3.3). However, since this API is now deprecated, we will exclusively work directly with the new interface, using Lightconer subclasses, in this tutorial.

import py21cmfast as p21c
from py21cmfast import plotting
import os
import matplotlib.pyplot as plt
from matplotlib import colors
from matplotlib import cm
import numpy as np
from scipy.spatial.transform import Rotation
from astropy import units as un
%matplotlib inline

print(f"21cmFAST version is {p21c.__version__}")

21cmFAST version is 3.2.1.dev196+g12e3ef7

Let’s set up a standard rectilinear lightcone subclass.

The final minimum redshift of the lightcone is set by the min_redshift argument, and the maximum redshift of the lightcone is defined by the max_redshift argument (note that this is not the maximum redshift coeval box computed, which is instead controlled by z_heat_max, merely the maximum saved into the returned lightcone).

You can specify which 3D quantities are interpolated as lightcones as well (these must be fields that exist on any standard output box).

user_params = p21c.UserParams(
    HII_DIM=64, BOX_LEN=256, KEEP_3D_VELOCITIES=True
)

lcn = p21c.RectilinearLightconer.with_equal_cdist_slices(
    min_redshift=7.0,
    max_redshift=12.0,
    quantities=('brightness_temp', 'density', 'velocity_z'),
    resolution=user_params.cell_size,
    # index_offset=0,
)

Notice the index_offset argument above. When creating the lightcone, the slice of the coeval box that makes the first slice of the lightcone is arbitrary. By setting index_offset to zero, we set the very back of the lightcone to be equal to the last slice of the coeval. The default is to have the first slice of the lightcone correspond to the first slice of the coeval.

lightcone = p21c.run_lightcone(
    lightconer=lcn,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=user_params
)

/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/inputs.py:510: UserWarning: The USE_INTERPOLATION_TABLES setting has changed in v3.1.2 to be default True. You can likely ignore this warning, but if you relied onhaving USE_INTERPOLATION_TABLES=False by *default*, please set it explicitly. To silence this warning, set it explicitly to True. Thiswarning will be removed in v4.
  warnings.warn(
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vx
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vy
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vz
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vx_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vy_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vz_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")

fig, ax = plt.subplots(1, 1, figsize=(12, 6), constrained_layout=True)
plotting.lightcone_sliceplot(lightcone, ax=ax, fig=fig);

fig, ax = plt.subplots(1, 1, figsize=(12, 6), constrained_layout=True)
plotting.lightcone_sliceplot(lightcone, 'density', ax=ax, fig=fig);

Simple!

Using an Angular Lightcone

Say we would like to save our lightcone in angular coordinates. We can do this with the AngularLightconer subclass. This class requires you to pass a set of latitude (from \(-\pi/2\) to \(\pi/2\)) and longitude (0 to 2\(\pi\)) points. It does not require these to be on a grid of any form (though often that is what you’ll want). A common use-case would be to use a regular grid in (lat, lon) that roughly covers the coeval box size (at the low-z end). Another common use-case would be to use some subset of a HEALPix map (or maybe the full map, if you’d like to go on and use these simulations to produce visibilities!).

This lightconer is able to deal with full sky or any subset.

In this case, for ease of visualization, we will set up a regular (lon, lat) grid that is similar in size and resolution to our coeval box. Since we’d like to see the “transverse face” of the lightcone and see that it looks like the face of the rectilinear lightcone, we need to set a few options. First, set the angular region:

cosmo = lightcone.cosmo_params.cosmo

box_size_radians = user_params.BOX_LEN / cosmo.comoving_distance(7.0).value

lon = np.linspace(0, box_size_radians, user_params.HII_DIM)
lat = np.linspace(0, box_size_radians, user_params.HII_DIM)[::-1]  # This makes the X-values increasing from 0.

LON, LAT = np.meshgrid(lon, lat)
LON = LON.flatten()
LAT = LAT.flatten()

So far, we have just specified a small grid that should roughly line up with the coeval box, but it is close to the \(x\)-axis. Natively, the rectilinear lightcone uses the last axis (the \(z\) axis) as the line-of-sight. We can use the options origin and rotation to align our specified co-ordinates with the rectilinear lightcone.

The first operation applied is the rotation. In our case, we want to rotate our angles from being close to the \(x\)-axis, to instead being close to the \(z\)-axis (imagine a small squarish region of pixels perpendicular to and close to the \(z\)-axis). This means we need to rotate the co-ordinates by 90 degrees around the \(y\)-axis. Secondly, we apply an offset. The rectilinear lightcone is constructed in such a way that the first slice at the lowest redshift corresponds to the first slice along the \(z\)-axis. We achieve this by displacing the origin by the comoving distance to the minimum redshift. Note also that we constructed our initial latitude grid to go from high elevation down to the horizon. When the grid is rotated up to be close to the \(z\)-axis, this gets the \(x\)-coordinates going from left-to-right instead of the other way around.

offset = cosmo.comoving_distance(lightcone.lightcone_redshifts.min()).to(un.pixel, un.pixel_scale(user_params.cell_size/un.pixel))
origin = np.array([0, 0, offset.value])*offset.unit
rot = Rotation.from_euler('Y', -np.pi/2)

Now we can construct the lightconer, using the AngularLightconer class. Note that we specify that we want to get the line-of-sight velocity lightcone as well as our other quantities. This will allow us to apply RSDs to the temperature lightcone at the end.

ang_lcn = p21c.AngularLightconer.with_equal_cdist_slices(
    min_redshift=7.0,
    max_redshift=12.0,
    quantities=('brightness_temp', 'density','velocity_x', 'velocity_y', 'velocity_z'),
    resolution=user_params.cell_size,
    latitude=LAT,
    longitude=LON,
    origin=-origin,
    rotation=rot,
    get_los_velocity=True
)

Notice that we set get_los_velocity to True. This ensures that we take the 3D (coeval) velocities, and make a lightcone out of them where the field is the line-of-sight component (really line-of-sight, not just the third dimension of the coeval field). Setting this to True switches on application of RSDs in the run_lightcone function below, once all the lightcones have been computed.

Since it is a common thing to want to get an angular lightcone that “looks like” the coeval box at some redshift, we provide a convenience method for creating it:

ang_lcn2 = p21c.AngularLightconer.like_rectilinear(
    match_at_z=7.0,
    max_redshift=12.0,
    user_params=user_params,
    get_los_velocity=True,
    quantities=('brightness_temp', 'density','velocity_x', 'velocity_y', 'velocity_z'),
)

This convenience method gives the same lightcone we manually got above:

ang_lcn == ang_lcn2

True

Once we have that, constructing the lightcone itself looks exactly the same as before, except that we need to turn off the RSDs. We apply RSDs to the angular lightcone itself, rather than to the coeval boxes that make it up. This is because the RSDs are applied in the transverse plane, and so we need to know the transverse co-ordinates of each pixel. This is not possible for a rectilinear lightcone, but is possible for an angular lightcone.

ang_lightcone = p21c.run_lightcone(
    lightconer=ang_lcn,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=user_params,
    flag_options={"APPLY_RSDS": False}
)

In this case, the output lightcone has a different type than the default rectilinear Lightcone:

type(ang_lightcone)

py21cmfast.outputs.AngularLightcone

The boxes in ang_lightcone have their transverse coordinates packed up into one dimension (because, again, there is no assumption that the angular coordinates were in a regular grid). In our specific case here, they are in a grid, so let’s reshape the array to make it look like our rectilinear lightcone:

bt = ang_lightcone.brightness_temp.reshape((user_params.HII_DIM , user_params.HII_DIM, ang_lightcone.brightness_temp.shape[-1]))

Now we can plot the rectilinear and angular lightcones beside eachother to check that they are similar. Here we simply plot the face of the lightcone at the lowest redshift:

fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,16))
ax[0].imshow(lightcone.brightness_temp[:, :, 0], origin='lower', extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[0].set_title("Rectilinear Lightcone Face")
ax[1].imshow(bt[:, :, 0], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[1].set_title("Angular Lightcone Face")

btrsd = ang_lightcone.lightcones['brightness_temp_with_rsds'].reshape((user_params.HII_DIM , user_params.HII_DIM, ang_lightcone.brightness_temp.shape[-1]))
ax[2].imshow(btrsd[:, :, 0], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[2].set_title("Angular Lightcone Face (RSDs)")

ax[0].set_xlabel("Transverse (x) dimension [cMpc]")
ax[1].set_xlabel("Transverse (x) dimension [cMpc]")
ax[0].set_ylabel("Transverse (y) dimension [cMpc]")

Text(0, 0.5, 'Transverse (y) dimension [cMpc]')

Here, we see that the upper left corner of the plots is very similar – this is where the two grids overlap almost exactly. Moving to the bottom right corner, the angular lightcone “curves away” from the rectilinear, and so the plots become less similar. Plotting the same thing but at the highest redshift:

fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,16))
ax[0].imshow(lightcone.brightness_temp[:, :, -1], origin='lower', extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[0].set_title("Rectilinear Lightcone Face")
ax[1].imshow(bt[:, :, -1], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[1].set_title("Angular Lightcone Face")

ax[2].imshow(btrsd[:, :, -1], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[2].set_title("Angular Lightcone Face (RSDs)")

ax[0].set_xlabel("Transverse (x) dimension [cMpc]")
ax[1].set_xlabel("Transverse (x) dimension [cMpc]")
ax[0].set_ylabel("Transverse (y) dimension [cMpc]")

Text(0, 0.5, 'Transverse (y) dimension [cMpc]')

We see that the upper left is still very similar, but changes much more radidly moving away from that point, as we would expect, since we have tried to ensure that the angular coordinates align at the lowest redshift, and the angular lightcone “fans out” towards higher redshift.

Redshift-Space Distortions in the Lightcone

To obtain RSDs in the angular lightcone requires setting the get_los_velocity=True on the lightconer instance. Contrary to the rectilinear lightcone, the angular lightcone only applies the RSDs to the lightcone itself, rather than the coeval boxes that make it up. It does this after the full lightcone has been created, and the result is saved in the “brightness_temp_with_rsds” element of the .lightcones dict.

Three “variants” of RSDs can be obtained: (i) No RSDs, (ii) linear-only RSDs, and (iii) linear + subcell RSDs. We compare these now:

lightcone_norsd = p21c.run_lightcone(
    lightconer=lcn,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=user_params,
    flag_options={"APPLY_RSDS": False}
)

lightcone_withsubcell = p21c.run_lightcone(
    lightconer=lcn,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=user_params,
    flag_options={"SUBCELL_RSD": True}
)

/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/inputs.py:510: UserWarning: The USE_INTERPOLATION_TABLES setting has changed in v3.1.2 to be default True. You can likely ignore this warning, but if you relied onhaving USE_INTERPOLATION_TABLES=False by *default*, please set it explicitly. To silence this warning, set it explicitly to True. Thiswarning will be removed in v4.
  warnings.warn(
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vx
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vy
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vz
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vx_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vy_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/emu/lib/python3.10/site-packages/py21cmfast/_utils.py:821: UserWarning: Trying to remove array that isn't yet created: hires_vz_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")

First, let’s look at the line-of-sight “apparent displacement” of the first lightcone slice incurred by the rectilinear and angular lightcones. For the rectilinear, this is simply the z-component of the velocity, while for the angular lightcone, it is a combination of the full 3D velocity, pointing towards the observer:

def topix(losv):
    if losv.ndim ==2:
        n = user_params.HII_DIM
        losv = losv.reshape((n,n, -1))

    los_displacement = (losv*un.Mpc/un.s /  ang_lightcone.cosmo_params.cosmo.H(7.0)).to(un.Mpc)
    los_displacement_pix = los_displacement.to(un.pixel, un.pixel_scale(un.pixel / user_params.cell_size))
    return los_displacement_pix.value

fig, ax = plt.subplots(1, 2, figsize=(12, 6), constrained_layout=True)
los_displacement_pix = topix(ang_lightcone.lightcones['los_velocity'])
cax =ax[0].imshow(los_displacement_pix[...,0], origin='lower', vmin=-1, vmax=1)

los_displacement_rec_pix = topix(lightcone.velocity_z)
cax = ax[1].imshow(-los_displacement_rec_pix[..., 0], origin='lower', vmin=-1, vmax=1)
plt.colorbar(cax, ax=ax[1])

<matplotlib.colorbar.Colorbar at 0x7f548e1d1c30>

Again, we see the same pattern here – the upper left is very similar, but the lower right differs substantially.

We can compute the difference between the angular and rectilinear lightcones (on the closest slice to the observer) for each combination of RSDs:

dvdx = np.gradient(-los_displacement_pix, ang_lightcone.lightcone_distances / (user_params.cell_size / un.pixel), axis=-1)

fig, ax = plt.subplots(3, 3, figsize=(12, 12), constrained_layout=True, sharex=True, sharey=True)

extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN)

btlin = bt[:, :, 0]/(1 + dvdx.value[...,0])

ang_norsd_lin = bt[:, :, 0] - btlin
ax[0, 0].imshow(ang_norsd_lin, origin='lower', extent=extent, vmin=-3, vmax=3)
ax[0, 0].set_title(r"$T_{\rm norsd} - T_{\rm linear}$")
ax[0, 0].set_ylabel("ANGULAR")

ang_lin_full = btlin - btrsd[:, :, 0]
cax = ax[0, 1].imshow(ang_lin_full, origin='lower',extent=extent, vmin=-3, vmax=3)
ax[0, 1].set_title(r"$T_{\rm linear} - T_{\rm fullrsd}$")

ang_norsd_full = bt[:, :, 0] - btrsd[:, :, 0]
cax = ax[0, 2].imshow(ang_norsd_full, origin='lower',extent=extent, vmin=-3, vmax=3)
ax[0, 2].set_title(r"$T_{\rm norsd} - T_{\rm fullrsd}$")


rct_norsd_lin = lightcone_norsd.brightness_temp[:, :, 0] - lightcone.brightness_temp[:, :, 0]
cax= ax[1, 0].imshow(rct_norsd_lin, origin='lower',extent=extent, vmin=-3, vmax=3)
ax[1, 0].set_ylabel("RECT")

rct_norsd_full = (lightcone_norsd.brightness_temp[:, :, 0] - lightcone_withsubcell.brightness_temp[:, :, 0])
cax= ax[1, 2].imshow(rct_norsd_full, origin='lower',extent=extent, vmin=-3, vmax=3)

rct_lin_full = (lightcone.brightness_temp[:, :, 0] - lightcone_withsubcell.brightness_temp[:, :, 0])
cax= ax[1, 1].imshow(rct_lin_full, origin='lower',extent=extent, vmin=-3, vmax=3)

ax[2, 0].set_ylabel("DIFF")
ax[2, 0].imshow(rct_norsd_lin - ang_norsd_lin, origin='lower',extent=extent, vmin=-3, vmax=3)
ax[2, 2].imshow(rct_norsd_full - ang_norsd_full, origin='lower',extent=extent, vmin=-3, vmax=3)
ax[2, 1].imshow(rct_lin_full - ang_lin_full, origin='lower',extent=extent, vmin=-3, vmax=3)

<matplotlib.image.AxesImage at 0x7f548db27af0>

The lower row shows the difference between angular and rectilinear, again showing that same pattern from top-left to lower-right. The left panels show the effect of just the linear RSD approximation, the middle panels show just the effect of the subcell RSDs, and the right panel shows the effect of both being applied.

Plotting evolution of the angular lightcone

We can also plot the evolution of the lightcone, choosing one transverse slice:

fig, ax = plt.subplots(2, 1, sharex=True)
ax[0].imshow(lightcone.brightness_temp[:, 0])
ax[1].imshow(bt[:, 0])

<matplotlib.image.AxesImage at 0x7f548416da50>

Here we see a very similar evolution, since we have chosen a slice that corresponds best to the rectilinear lightcone.

We can get a more accurate representation of the angular lightcone by plotting it in polar coordinates. Here, we can just plot a slice corresponding to the horizon:

fig, ax = plt.subplots(1,1,figsize=[14, 5], constrained_layout=True,subplot_kw={'projection': 'polar'})

ax.set_thetamax(lat.max() * 180/np.pi)
ax.set_ylim(ang_lcn.lc_distances.min().value, ang_lcn.lc_distances.max().value)
ax.set_rorigin(0)
H, D = np.meshgrid(lat, ang_lcn.lc_distances.value)
plt.grid(False)
img = ax.pcolormesh(H, D, bt[:, 0].T,edgecolors='face', cmap='EoR', vmin=-150, vmax=30)
plt.colorbar(img)

<matplotlib.colorbar.Colorbar at 0x7f5483e9c850>

Here, the radial axis corresponds to comoving distance, and we have used the built-in “EoR” colormap provided by 21cmFAST.

Saving the Lightcones

You can easily save a HDF5 file containing the interpolated lightcone:

filename = lightcone.save(direc='_cache', clobber=True)

print(os.path.basename(filename))

And for the angular lightcone:

ang_lightcone.save(fname='big_angular_lightcone.h5')

'/data4/smurray/Projects/radio/EOR/OthersCodes/21cmFAST/21cmFAST-git/docs/tutorials/big_angular_lightcone.h5'

ang_filename = ang_lightcone.save(direc='_cache', clobber=True)

print(os.path.basename(ang_filename))

Overcoming blurry edges

Since the angular lightcone has to be interpolated (in the angular direction) from a rectilinear coeval box, it will tend to “blur” edges of ionized regions that would otherwise be sharp. This will reduce the power on small scales compared to the truth. To overcome this, you can try setting the resolution higher, or the interpolation order, or both:

lon = np.linspace(0, box_size_radians, user_params.HII_DIM)
lat = np.linspace(0, box_size_radians, user_params.HII_DIM)[::-1]  # This makes the X-values increasing from 0.

lonhires = np.linspace(0, box_size_radians, user_params.HII_DIM*4)
lathires = np.linspace(0, box_size_radians, user_params.HII_DIM*4)[::-1]  # This makes the X-values increasing from 0.

LON, LAT = np.meshgrid(lon, lat)
LON = LON.flatten()
LAT = LAT.flatten()

LONhires, LAThires = np.meshgrid(lonhires, lathires)
LONhires = LONhires.flatten()
LAThires = LAThires.flatten()

ang_lcn_cubic = p21c.AngularLightconer.with_equal_cdist_slices(
    min_redshift=7.0,
    max_redshift=7.2,
    quantities=('brightness_temp', 'density','velocity_x', 'velocity_y', 'velocity_z'),
    resolution=user_params.cell_size,
    latitude=LAT,
    longitude=LON,
    origin=-origin,
    rotation=rot,
    get_los_velocity=True,
    interpolation_order=3,
)

ang_lcn_hires = p21c.AngularLightconer.with_equal_cdist_slices(
    min_redshift=7.0,
    max_redshift=7.2,
    quantities=('brightness_temp', 'density','velocity_x', 'velocity_y', 'velocity_z'),
    resolution=user_params.cell_size,
    latitude=LAThires,
    longitude=LONhires,
    origin=-origin,
    rotation=rot,
    get_los_velocity=True
)

ang_lightcone_hires = p21c.run_lightcone(
    lightconer=ang_lcn_hires,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=user_params,
    flag_options={"APPLY_RSDS": False}
)

ang_lightcone_cubic = p21c.run_lightcone(
    lightconer=ang_lcn_cubic,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=user_params,
    flag_options={"APPLY_RSDS": False}
)

fig, ax = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12,20), constrained_layout=True)

bt= ang_lightcone.brightness_temp[:, 0].reshape(user_params.HII_DIM, user_params.HII_DIM)
bthires= ang_lightcone_hires.brightness_temp[:, 0].reshape(user_params.HII_DIM*4, user_params.HII_DIM*4)

n = user_params.HII_DIM

ax[0].imshow(lightcone.brightness_temp[n//2:, :n//2, 0], origin='lower', extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[0].set_title("Rectilinear Lightcone Face")
ax[1].imshow(bt[n//2:,:n//2], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[1].set_title("Angular (Same Res)")

ax[2].imshow(bthires[n*2:,:n*2], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[2].set_title("Angular (Hi Res)")


btcubic= ang_lightcone_cubic.brightness_temp[:, 0].reshape(user_params.HII_DIM, user_params.HII_DIM)

ax[3].imshow(btcubic[n//2:,:n//2], origin='lower',extent=(0, user_params.BOX_LEN, 0, user_params.BOX_LEN), cmap='EoR', vmin=-150, vmax=30)
ax[3].set_title("Angular (cubic interp)")


ax[0].set_xlabel("Transverse (x) dimension [cMpc]")
ax[1].set_xlabel("Transverse (x) dimension [cMpc]")
ax[0].set_ylabel("Transverse (y) dimension [cMpc]")

Text(0, 0.5, 'Transverse (y) dimension [cMpc]')

Higher Resolution Run For Visualization

if os.environ.get("TEST_RUN", "False") == "True":
    print("Skipping the rest of the notebook because TEST_RUN is True")
    raise SystemExit

So far, we have constructed our angular lightcone to be similar in size and resolution to the rectilinear lightcone. But as we’ve said, there is no restriction on the angular coordinates that can be put in. You could input the coordinates from a full-sky healpix map. To illustrate using a wider set of coordinates, let’s create a “lightplane” taking a set of regular coordinates on the horizon, extending back into the line of sight. We do this for a wide angle (~22.5 degrees, or about 6 of our boxes stacked). In this case, we use no rotation or origin offset:

lath = np.zeros(3000)
lonh = np.linspace(0, np.pi/8, 3000)

For this last run, let’s setup a much more realistically-sized simulation, to get a nicer visualization:

hires_user_params = p21c.UserParams(
    BOX_LEN=900.0, HII_DIM=600, DIM=1200, KEEP_3D_VELOCITIES=True
)
hires_flags = p21c.FlagOptions(
    USE_TS_FLUCT=True,
    APPLY_RSDS=False
)

ang_lcn_horiz = p21c.AngularLightconer.with_equal_cdist_slices(
    min_redshift=6.0,
    max_redshift=30.0,
    quantities=('brightness_temp', 'density'),
    resolution=hires_user_params.cell_size,
    latitude=lath,
    longitude=lonh,
    get_los_velocity=True
)

lc_horiz = ang_lightcone = p21c.run_lightcone(
    lightconer=ang_lcn_horiz,
    global_quantities=("brightness_temp", 'density', 'xH_box'),
    direc='_cache',
    user_params=hires_user_params,
    flag_options=hires_flags,
    regenerate=False
)

/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/inputs.py:510: UserWarning: The USE_INTERPOLATION_TABLES setting has changed in v3.1.2 to be default True. You can likely ignore this warning, but if you relied onhaving USE_INTERPOLATION_TABLES=False by *default*, please set it explicitly. To silence this warning, set it explicitly to True. Thiswarning will be removed in v4.
  warnings.warn(
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/_utils.py:815: UserWarning: Trying to remove array that isn't yet created: hires_vx
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/_utils.py:815: UserWarning: Trying to remove array that isn't yet created: hires_vy
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/_utils.py:815: UserWarning: Trying to remove array that isn't yet created: hires_vz
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/_utils.py:815: UserWarning: Trying to remove array that isn't yet created: hires_vx_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/_utils.py:815: UserWarning: Trying to remove array that isn't yet created: hires_vy_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/py21cmfast/_utils.py:815: UserWarning: Trying to remove array that isn't yet created: hires_vz_2LPT
  warnings.warn(f"Trying to remove array that isn't yet created: {k}")
/home/smurray/miniconda3/envs/21cmfast/lib/python3.10/site-packages/cosmotile/__init__.py:391: UserWarning: Install numba for a speedup of cloud_in_cell
  fine_field = cloud_in_cell_los(fine_field, fine_rsd)

When plotting, we add some grid lines identifying roughly the positions of the stacked underlying boxes (the one closest to the actual stacked coeval position is the one closest to \(\theta=0\)).

from astropy.cosmology import z_at_value
from astropy import units

def get_eor_cmap(vmin=-150, vmax=30):
    name = f"EoR-{vmin}-{vmax}"
    negative_segments = 4
    positive_segments = 2
    neg_frac = abs(vmin) / (vmax - vmin)
    neg_seg_size = neg_frac / negative_segments
    pos_frac = abs(vmax) / (vmax - vmin)
    pos_seg_size = pos_frac / positive_segments
    eor_colour = colors.LinearSegmentedColormap.from_list(
        name,
        [
            (0, "white"),
            (neg_seg_size, "yellow"),
            (neg_seg_size*2, "orange"),
            (neg_seg_size*3, "maroon"),
            (neg_seg_size*4, "black"),
            (neg_seg_size*4 + pos_seg_size, "#1f77b4"),
            (1, "skyblue"),
        ],
    )
    try:
        plt.register_cmap(cmap=eor_colour)
    except ValueError:
        cm.unregister_cmap(name)
        plt.register_cmap(cmap=eor_colour)

    return name

cmap = get_eor_cmap(vmin, vmax)

plt.style.use('dark_background')

fig, ax = plt.subplots(1,1,figsize=[12, 10], constrained_layout=True,subplot_kw={'projection': 'polar'})

vmin, vmax = lc_horiz.brightness_temp.min(), 30.0
cmap = get_eor_cmap(vmin, vmax)

theta = hires_user_params.BOX_LEN / ang_lcn_horiz.cosmo.comoving_distance(ang_lcn_horiz.lc_redshifts.min()).value
nchunks = int(lonh.max() / theta) + 1
ax.set_thetamax(lonh.max() * 180/np.pi)
ax.set_ylim(ang_lcn_horiz.lc_distances.min().value, ang_lcn_horiz.lc_distances.max().value)
ax.set_rorigin(0)
H, D = np.meshgrid(lonh, ang_lcn_horiz.lc_distances.value)
img = ax.pcolormesh(H, D, lc_horiz.brightness_temp.T, edgecolors='face', cmap=cmap, vmin=vmin, vmax=vmax)
cax = plt.colorbar(img)
cax.set_label(r"$T_{\rm B}$ [mK]")
ax.set_thetagrids(np.arange(nchunks)*theta * 180/np.pi)
rgrid = np.arange(0, D.max(), hires_user_params.BOX_LEN)
rgrid = rgrid[rgrid > D.min()]
rgrid = np.append(rgrid, D.max())
ax.set_rgrids(rgrid, labels=[f"{z_at_value(ang_lcn_horiz.cosmo.comoving_distance, d*units.Mpc).value:.1f}" for d in rgrid])
ax.grid(True)
#ax.set_xlabel("Redshift")
ax.text(0.5, -0.06, "Redshift", fontsize=14, transform=ax.transAxes)

/tmp/ipykernel_56332/325665923.py:14: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  img = ax.pcolormesh(H, D, lc_horiz.brightness_temp.T, edgecolors='face', cmap=cmap, vmin=vmin, vmax=vmax)

Text(0.5, -0.06, 'Redshift')