diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/README b/examples/RadiativeTransferTests/CosmoAdvection_1D/README
new file mode 100644
index 0000000000000000000000000000000000000000..44b286830c0ca8bf3d1b599582dfd5835b898ba1
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/README
@@ -0,0 +1,29 @@
+1D advection test for cosmological radiative transfer.
+
+Test that your method is TVD and the propagation speed of the photons is
+correct. The ICs set up three photon groups: 
+- The first is a top hat function initial distribution where outside values
+  are zero.
+- The second is a top hat function initial distribution where outside values
+  are nonzero. This distinction is important to test because photon energies 
+  can't be negative, so these cases need to be tested individually.
+- the third is a smooth Gaussian. 
+
+This way, you can test multiple initial condition scenarios simultaneously. 
+There are no stars to act as sources. Also make sure that you choose your
+photon frequencies in a way that doesn't interact with gas!
+
+There are 3 parameter files available:
+- rt_advection1D_low_redshift.yml for a test at very low redshift (almost no cosmological expansion)
+- rt_advection1D_medium_redshift.yml for a test at moderate redshift (z=10)
+- rt_advection1D_high_redshift.yml for a test at high redshift (z=110)
+
+Select one and give it to the run.sh script via command line argument:
+    ./run.sh rt_advection1D_[redshift_of_your_choice].yml
+
+
+The ICs are created to be compatible with GEAR_RT and SPHM1RT. Recommended configuration:
+GEAR_RT:
+    --with-rt=GEAR_3 --with-rt-riemann-solver=GLF --with-hydro-dimension=1 --with-hydro=gizmo-mfv \
+     --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
+
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/makeIC.py b/examples/RadiativeTransferTests/CosmoAdvection_1D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..d4428ba9c7ba27b4dbfec93cebcea68355822287
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/makeIC.py
@@ -0,0 +1,213 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+#               2022 Tsang Keung Chan (chantsangkeung@gmail.com)
+#               2024 Stan Verhoeve (s06verhoeve@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+
+# -----------------------------------------------------------
+# Add initial conditions for photon energies and fluxes
+# for 1D advection of photons.
+# First photon group: Top hat function with zero as the
+#       baseline.
+# Second photon group: Top hat function with nonzero value
+#       as the baseline.
+# Third photon group: Gaussian.
+# -----------------------------------------------------------
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+# define unit system to use
+unitsystem = cosmo_units
+
+# Box is 260 Mpc
+boxsize = 260 * unyt.Mpc
+boxsize = boxsize.to(unitsystem["length"])
+
+reduced_speed_of_light_fraction = 1.0
+
+
+# number of photon groups
+nPhotonGroups = 3
+
+# number of particles in each dimension
+n_p = 1000
+
+# filename of ICs to be generated
+outputfilename = "advection_1D.hdf5"
+
+
+def initial_condition(x, unitsystem):
+    """
+    The initial conditions that will be advected
+
+    x: particle position. 3D unyt array
+    unitsystem: Currently used unitsystem.
+
+    returns:
+    E: photon energy density for each photon group. List of scalars with size of nPhotonGroups
+    F: photon flux for each photon group. List with size of nPhotonGroups of numpy arrays of shape (3,)
+    """
+
+    # you can make the photon quantities unitless, the units will
+    # already have been written down in the writer.
+    # However, that means that you need to convert them manually.
+
+    unit_energy = (
+        unitsystem["mass"] * unitsystem["length"] ** 2 / unitsystem["time"] ** 2
+    )
+    unit_velocity = unitsystem["length"] / unitsystem["time"]
+    unit_flux = unit_energy * unit_velocity
+
+    c_internal = (unyt.c * reduced_speed_of_light_fraction).to(unit_velocity)
+
+    E_list = []
+    F_list = []
+
+    # Group 1 Photons:
+    # -------------------
+
+    if x[0] < 0.33 * boxsize:
+        E = 0.0 * unit_energy
+    elif x[0] < 0.66 * boxsize:
+        E = 1.0 * unit_energy
+    else:
+        E = 0.0 * unit_energy
+
+    # Assuming all photons flow in only one direction
+    # (optically thin regime, "free streaming limit"),
+    #  we have that |F| = c * E
+    F = np.zeros(3, dtype=np.float64)
+    F[0] = (E * c_internal).to(unit_flux)
+
+    E1 = E
+    F1 = F[0]
+
+    E_list.append(E)
+    F_list.append(F)
+
+    # Group 2 Photons:
+    # -------------------
+
+    if x[0] < 0.33 * boxsize:
+        E = 1.0 * unit_energy
+    elif x[0] < 0.66 * boxsize:
+        E = 3.0 * unit_energy
+    else:
+        E = 1.0 * unit_energy
+
+    F = np.zeros(3, dtype=np.float64)
+
+    F[0] = (E * c_internal).to(unit_flux)
+
+    E_list.append(E)
+    F_list.append(F)
+
+    E2 = E
+    F2 = F[0]
+
+    # Group 3 Photons:
+    # -------------------
+    sigma = 0.1 * boxsize
+    mean = 0.5 * boxsize
+    amplitude = 2.0
+
+    E = amplitude * np.exp(-(x[0] - mean) ** 2 / (2 * sigma ** 2)) * unit_energy
+
+    F = np.zeros(3, dtype=np.float64)
+    F[0] = (E * c_internal).to(unit_flux)
+
+    E_list.append(E)
+    F_list.append(F)
+
+    E3 = E
+    F3 = F[0]
+
+    return E_list, F_list
+
+
+if __name__ == "__main__":
+
+    xp = unyt.unyt_array(np.zeros((n_p, 3), dtype=np.float64), boxsize.units)
+
+    dx = boxsize / n_p
+
+    for i in range(n_p):
+        xp[i, 0] = (i + 0.5) * dx
+
+    w = Writer(unitsystem, boxsize, dimension=1)
+
+    w.gas.coordinates = xp
+    w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
+    mpart = 1e20 * unyt.M_Sun
+    mpart = mpart.to(unitsystem["mass"])
+
+    w.gas.masses = np.ones(xp.shape[0], dtype=np.float64) * mpart
+    w.gas.internal_energy = (
+        np.ones(xp.shape[0], dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g
+    )
+
+    # Generate initial guess for smoothing lengths based on MIPS
+    w.gas.generate_smoothing_lengths(boxsize=boxsize, dimension=1)
+
+    # If IDs are not present, this automatically generates
+    w.write(outputfilename)
+
+    # Now open file back up again and add photon groups
+    # you can make them unitless, the units have already been
+    # written down in the writer. In this case, it's in cgs.
+
+    F = h5py.File(outputfilename, "r+")
+    header = F["Header"]
+    nparts = header.attrs["NumPart_ThisFile"][0]
+    parts = F["/PartType0"]
+
+    for grp in range(nPhotonGroups):
+        dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+        energydata = np.zeros((nparts), dtype=np.float32)
+        parts.create_dataset(dsetname, data=energydata)
+
+        dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+        fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+        parts.create_dataset(dsetname, data=fluxdata)
+
+    for p in range(nparts):
+        E, Flux = initial_condition(xp[p], unitsystem)
+        for g in range(nPhotonGroups):
+            Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1)
+            parts[Esetname][p] = E[g]
+            Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+            parts[Fsetname][p] = Flux[g]
+
+    #  from matplotlib import pyplot as plt
+    #  plt.figure()
+    #  for g in range(nPhotonGroups):
+    #      #  Esetname = "PhotonEnergiesGroup{0:d}".format(g+1)
+    #      #  plt.plot(xp[:,0], parts[Esetname], label="E "+str(g+1))
+    #      Fsetname = "PhotonFluxesGroup{0:d}".format(g+1)
+    #      plt.plot(xp[:,0], parts[Fsetname][:,0], label="F "+str(g+1))
+    #  plt.legend()
+    #  plt.show()
+
+    F.close()
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/plotEnergy.py b/examples/RadiativeTransferTests/CosmoAdvection_1D/plotEnergy.py
new file mode 100755
index 0000000000000000000000000000000000000000..7dc1bff0186f740a01a5a2d6ddfe547a757f9e79
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/plotEnergy.py
@@ -0,0 +1,211 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 Stan Verhoeve (s06verhoeve@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+import os
+import sys
+import argparse
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+import unyt
+from matplotlib import pyplot as plt
+
+# Parameters users should/may tweak
+snapshot_base = "output"  # Snapshot basename
+plot_physical_quantities = True  # Plot physical or comoving quantities
+
+time_first = 0  # Time of first snapshot
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-z", "--redshift", help="Redshift domain to plot advection for", default="high"
+    )
+
+    args = parser.parse_args()
+    return args
+
+
+def get_snapshot_list(snapshot_basename="output"):
+    """
+    Find the snapshot(s) that are to be plotted
+    and return their names as list
+    """
+
+    snaplist = []
+
+    dirlist = os.listdir()
+    for f in dirlist:
+        if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+            snaplist.append(f)
+
+    snaplist = sorted(snaplist)
+    if len(snaplist) == 0:
+        print(f"No snapshots with base {snapshot_basename} found!")
+        sys.exit(1)
+
+    return snaplist
+
+
+def plot_param_over_time(snapshot_list, param="energy density", redshift_domain="high"):
+    print(f"Now plotting {param} over time")
+    # Grab number of photon groups
+    data = swiftsimio.load(snapshot_list[0])
+    meta = data.metadata
+    ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"][0])
+
+    # Number of rows and columns
+    nrows = 1 + int(plot_physical_quantities)
+    ncols = ngroups
+    # Create figure and axes
+    fig, axs = plt.subplots(nrows, ncols, figsize=(5.04 * ncols, 5.4 * nrows), dpi=200)
+
+    # Iterate over all photon groups
+    for n in range(ngroups):
+        # Arrays to keep track of plot_param and scale factor
+        plot_param = [[], []]
+        scale_factor = []
+        analytic_exponent = [0, 0]
+
+        # Functions to convert between scale factor and redshift
+        a2z = lambda a: 1 / a - 1
+        z2a = lambda z: 1 / (z + 1)
+
+        for file in snapshot_list:
+            data = swiftsimio.load(file)
+            meta = data.metadata
+
+            # Read comoving variables
+            energy = getattr(data.gas.photon_energies, f"group{n+1}")
+            mass = data.gas.masses
+            rho = data.gas.densities
+            volume = mass / rho
+
+            energy_density = energy / volume
+
+            if plot_physical_quantities:
+                # The SWIFT cosmology module assumes 3-dimensional lengths and volumes,
+                # so multiply by a**2 to get the correct relations
+                physical_energy_density = (
+                    energy_density.to_physical() * meta.scale_factor ** 2
+                )
+                physical_energy = energy.to_physical()
+
+                if param == "energy density":
+                    plot_param[1].append(
+                        1
+                        * physical_energy_density.sum()
+                        / physical_energy_density.shape[0]
+                    )
+                    analytic_exponent[1] = -1.0
+                elif param == "total energy":
+                    plot_param[1].append(1 * physical_energy.sum())
+                    analytic_exponent[1] = 0.0
+
+            if param == "energy density":
+                plot_param[0].append(1 * energy_density.sum() / energy_density.shape[0])
+                analytic_exponent[0] = 0.0
+            elif param == "total energy":
+                plot_param[0].append(1 * energy.sum())
+                analytic_exponent[0] = 0.0
+
+            scale_factor.append(meta.scale_factor)
+
+        if param == "energy density":
+            titles = [
+                "Comoving energy density",
+                "Physical energy density $\\times a^2$",
+            ]
+            ylabel = "Average energy density"
+            figname = "output_energy_density_over_time"
+        elif param == "total energy":
+            titles = ["Comoving total energy", "Physical total energy"]
+            ylabel = "Total energy"
+            figname = "output_total_energy_over_time"
+
+        # Analytic scale factor
+        analytic_scale_factor = np.linspace(min(scale_factor), max(scale_factor), 1000)
+
+        for i in range(nrows):
+            ax = axs[i, n]
+            ax.scatter(scale_factor, plot_param[i], label="Simulation")
+
+            # Analytic scale factor relation
+            analytic = analytic_scale_factor ** analytic_exponent[i]
+
+            # Scale solution to correct offset
+            analytic = analytic / analytic[0] * plot_param[i][0]
+            ax.plot(
+                analytic_scale_factor,
+                analytic,
+                c="r",
+                label=f"Analytic solution $\propto a^{{{analytic_exponent[i]}}}$",
+            )
+
+            ax.legend()
+            ax.set_title(titles[i] + f" group {n+1}")
+
+            ax.set_xlabel("Scale factor")
+            secax = ax.secondary_xaxis("top", functions=(a2z, z2a))
+            secax.set_xlabel("Redshift")
+
+            ax.yaxis.get_offset_text().set_position((-0.05, 1))
+
+            if analytic_exponent[i] == 0.0:
+                ax.set_ylim(plot_param[i][0] * 0.95, plot_param[i][0] * 1.05)
+                ylabel_scale = ""
+            else:
+                ylabel_scale = "$\\times a^2$"
+            if n == 0:
+                units = plot_param[i][0].units.latex_representation()
+                ax.set_ylabel(f"{ylabel} [${units}$] {ylabel_scale}")
+    plt.tight_layout()
+    plt.savefig(f"{figname}-{redshift_domain}.png")
+    plt.close()
+
+
+if __name__ in ("__main__"):
+    # Get command line args
+    args = parse_args()
+    domain = args.redshift.lower()
+    if domain in ("low", "l", "low_redshift", "low redshift", "low-redshift"):
+        redshift_domain = "low_redshift"
+    elif domain in (
+        "medium",
+        "m",
+        "medium_redshift",
+        "medium redshift",
+        "medium-redshift",
+    ):
+        redshift_domain = "medium_redshift"
+    elif domain in ("high", "h", "high_redshift", "high redshift", "high-redshift"):
+        redshift_domain = "high_redshift"
+    else:
+        print("Redshift domain not recognised!")
+        sys.exit(1)
+
+    print(snapshot_base + f"_{redshift_domain}")
+    snaplist = get_snapshot_list(snapshot_base + f"_{redshift_domain}")
+
+    for param in ["energy density", "total energy"]:
+        plot_param_over_time(snaplist, param, redshift_domain)
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/plotSolution.py b/examples/RadiativeTransferTests/CosmoAdvection_1D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..cfd051d812245ec6318ae5dad324bd65f12a779c
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/plotSolution.py
@@ -0,0 +1,421 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+#               2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+
+# ----------------------------------------------
+# plot photon data assuming a 1D problem
+# give snapshot number as cmdline arg to plot
+# single snapshot, otherwise this script plots
+# all snapshots available in the workdir
+# ----------------------------------------------
+
+import os
+import sys
+import argparse
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+import unyt
+from matplotlib import pyplot as plt
+
+# Parameters users should/may tweak
+plot_all_data = True  # plot all groups and all photon quantities
+snapshot_base = "output"  # snapshot basename
+fancy = True  # fancy up the plots a bit
+plot_analytical_solutions = True  # overplot analytical solution
+plot_physical_quantities = True  # Plot physiical or comoving energy densities
+
+time_first = 0
+
+# properties for all scatterplots
+scatterplot_kwargs = {
+    "facecolor": "red",
+    "s": 4,
+    "alpha": 0.6,
+    "linewidth": 0.0,
+    "marker": ".",
+}
+
+# properties for all analytical solution plots
+analytical_solution_kwargs = {"linewidth": 1.0, "ls": "--", "c": "k", "alpha": 0.5}
+
+# -----------------------------------------------------------------------
+
+if plot_analytical_solutions:
+    from makeIC import initial_condition
+
+mpl.rcParams["text.usetex"] = True
+plot_all = False  # plot all snapshots
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-n", "--snapshot-number", help="Number of snaphot to plot", type=int
+    )
+    parser.add_argument(
+        "-z",
+        "--redshift",
+        help="Redshift domain to plot advection for",
+        default="high_redshift",
+    )
+
+    args = parser.parse_args()
+    return args
+
+
+def get_snapshot_list(snapshot_basename="output"):
+    """
+    Find the snapshot(s) that are to be plotted 
+    and return their names as list
+    """
+
+    snaplist = []
+
+    if plot_all:
+        dirlist = os.listdir()
+        for f in dirlist:
+            if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+                snaplist.append(f)
+
+        snaplist = sorted(snaplist)
+        if len(snaplist) == 0:
+            print(f"No snapshots with base {snapshot_basename} found!")
+            sys.exit(1)
+
+    else:
+        fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+        if not os.path.exists(fname):
+            print("Didn't find file", fname)
+            quit(1)
+        snaplist.append(fname)
+
+    return snaplist
+
+
+def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
+    """
+    Create the actual plot.
+
+    filename: file to work with
+    energy_boundaries:  list of [E_min, E_max] for each photon group. 
+                        If none, limits are set automatically.
+    flux_boundaries:    list of [F_min, F_max] for each photon group. 
+                        If none, limits are set automatically.
+    """
+    global time_first
+    print("working on", filename)
+
+    # Read in data firt
+    data = swiftsimio.load(filename)
+    meta = data.metadata
+    cosmology = meta.cosmology_raw
+
+    scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+
+    boxsize = meta.boxsize[0]
+
+    ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"][0])
+    # Currently, SPHM1RT only works for frequency group = 4 in the code
+    # However, we only plot 3 frequency groups here, so
+    # we set ngroups = 3:
+    if scheme.startswith("SPH M1closure"):
+        ngroups = 3
+
+    for g in range(ngroups):
+        # workaround to access named columns data with swiftsimio visualisaiton
+        new_attribute_str = "radiation_energy" + str(g + 1)
+        en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+        setattr(data.gas, new_attribute_str, en)
+
+        if plot_all_data:
+            # prepare also the fluxes
+            for direction in ["X"]:
+                new_attribute_str = "radiation_flux" + str(g + 1) + direction
+                f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+                setattr(data.gas, new_attribute_str, f)
+
+    part_positions = data.gas.coordinates[:, 0].copy()
+
+    # get analytical solutions
+    if plot_analytical_solutions:
+        # Grab unit system used in IC
+        IC_units = swiftsimio.units.cosmo_units
+
+        # Grab unit system used in snapshot
+        snapshot_units = meta.units
+
+        snapshot_unit_energy = (
+            snapshot_units.mass * snapshot_units.length ** 2 / snapshot_units.time ** 2
+        )
+
+        time = meta.time.copy()
+
+        # Take care of simulation time not starting at 0
+        if time_first == 0:
+            time_first = time
+            time = 0 * snapshot_units.time
+        else:
+            time -= time_first
+
+        speed = meta.reduced_lightspeed
+
+        advected_positions = data.gas.coordinates[:].copy()
+        advected_positions[:, 0] -= speed * time
+        nparts = advected_positions.shape[0]
+        # add periodicity corrections
+        negatives = advected_positions < 0.0
+        if negatives.any():
+            while advected_positions.min() < 0.0:
+                advected_positions[negatives] += boxsize
+        overshooters = advected_positions > boxsize
+        if overshooters.any():
+            while advected_positions.max() > boxsize:
+                advected_positions[overshooters] -= boxsize
+
+        analytical_solutions = np.zeros((nparts, ngroups), dtype=np.float64)
+        for p in range(part_positions.shape[0]):
+            E, F = initial_condition(advected_positions[p], IC_units)
+            for g in range(ngroups):
+                analytical_solutions[p, g] = E[g] / snapshot_unit_energy
+
+    # Plot plot plot!
+    if plot_all_data:
+
+        fig = plt.figure(figsize=(5.05 * ngroups, 5.4), dpi=200)
+        figname = filename[:-5] + f"-all-quantities.png"
+
+        for g in range(ngroups):
+
+            # plot energy
+            new_attribute_str = "radiation_energy" + str(g + 1)
+            photon_energy = getattr(data.gas, new_attribute_str)
+
+            # physical_photon_energy = photon_energy.to_physical()
+            ax = fig.add_subplot(2, ngroups, g + 1)
+            s = np.argsort(part_positions)
+            if plot_analytical_solutions:
+                ax.plot(
+                    part_positions[s],
+                    analytical_solutions[s, g],
+                    **analytical_solution_kwargs,
+                    label="analytical solution",
+                )
+            ax.scatter(
+                part_positions, photon_energy, **scatterplot_kwargs, label="simulation"
+            )
+            ax.legend()
+
+            ax.set_title("Group {0:2d}".format(g + 1))
+            if g == 0:
+                ax.set_ylabel(
+                    "Energies [$" + photon_energy.units.latex_representation() + "$]"
+                )
+            ax.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+            if energy_boundaries is not None:
+
+                if abs(energy_boundaries[g][1]) > abs(energy_boundaries[g][0]):
+                    fixed_min = energy_boundaries[g][0] - 0.1 * abs(
+                        energy_boundaries[g][1]
+                    )
+                    fixed_max = energy_boundaries[g][1] * 1.1
+                else:
+                    fixed_min = energy_boundaries[g][0] * 1.1
+                    fixed_max = energy_boundaries[g][1] + 0.1 * abs(
+                        energy_boundaries[g][0]
+                    )
+                ax.set_ylim(fixed_min, fixed_max)
+
+            # plot flux X
+            new_attribute_str = "radiation_flux" + str(g + 1) + "X"
+            photon_flux = getattr(data.gas, new_attribute_str)
+            physical_photon_flux = photon_flux.to_physical()
+
+            if scheme.startswith("GEAR M1closure"):
+                photon_flux = photon_flux.to("erg/cm**2/s")
+            elif scheme.startswith("SPH M1closure"):
+                photon_flux = photon_flux.to("erg*cm/s")
+            else:
+                print("Error: Unknown RT scheme " + scheme)
+                exit()
+            ax = fig.add_subplot(2, ngroups, g + 1 + ngroups)
+            ax.scatter(part_positions, photon_flux, **scatterplot_kwargs)
+
+            if g == 0:
+                ax.set_ylabel(
+                    "Flux X [$" + photon_flux.units.latex_representation() + "$]"
+                )
+            ax.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+            if flux_boundaries is not None:
+
+                if abs(flux_boundaries[g][1]) > abs(flux_boundaries[g][0]):
+                    fixed_min = flux_boundaries[g][0] - 0.1 * abs(flux_boundaries[g][1])
+                    fixed_max = flux_boundaries[g][1] * 1.1
+                else:
+                    fixed_min = flux_boundaries[g][0] * 1.1
+                    fixed_max = flux_boundaries[g][1] + 0.1 * abs(flux_boundaries[g][0])
+
+                ax.set_ylim(fixed_min, fixed_max)
+
+    else:  # plot just energies
+
+        fig = plt.figure(figsize=(5 * ngroups, 5), dpi=200)
+        figname = filename[:-5] + ".png"
+
+        for g in range(ngroups):
+
+            ax = fig.add_subplot(1, ngroups, g + 1)
+
+            new_attribute_str = "radiation_energy" + str(g + 1)
+            photon_energy = getattr(data.gas, new_attribute_str).to_physical()
+
+            s = np.argsort(part_positions)
+            if plot_analytical_solutions:
+                ax.plot(
+                    part_positions[s],
+                    analytical_solutions[s, g],
+                    **analytical_solution_kwargs,
+                    label="analytical solution",
+                )
+            ax.scatter(
+                part_positions, photon_energy, **scatterplot_kwargs, label="simulation"
+            )
+            ax.set_title("Group {0:2d}".format(g + 1))
+            if g == 0:
+                ax.set_ylabel(
+                    "Energies [$" + photon_energy.units.latex_representation() + "$]"
+                )
+            ax.set_xlabel("x [$" + part_positions.units.latex_representation() + "$]")
+
+            if energy_boundaries is not None:
+
+                if abs(energy_boundaries[g][1]) > abs(energy_boundaries[g][0]):
+                    fixed_min = energy_boundaries[g][0] - 0.1 * abs(
+                        energy_boundaries[g][1]
+                    )
+                    fixed_max = energy_boundaries[g][1] * 1.1
+                else:
+                    fixed_min = energy_boundaries[g][0] * 1.1
+                    fixed_max = energy_boundaries[g][1] + 0.1 * abs(
+                        energy_boundaries[g][0]
+                    )
+                ax.set_ylim(fixed_min, fixed_max)
+
+    # add title
+    title = filename.replace("_", r"\_")  # exception handle underscore for latex
+    if meta.cosmology is not None:
+        title += ", $z$ = {0:.2e}".format(meta.z)
+    title += ", $t$ = {0:.3e}".format(1 * meta.time)
+    fig.suptitle(title)
+
+    plt.tight_layout()
+    plt.savefig(figname)
+    plt.close()
+    return
+
+
+def get_minmax_vals(snaplist):
+    """
+    Find minimal and maximal values for energy and flux,
+    so you can fix axes limits over all snapshots
+
+    snaplist: list of snapshot filenames
+
+    returns:
+
+    energy_boundaries: list of [E_min, E_max] for each photon group
+    flux_boundaries: list of [Fx_min, Fy_max] for each photon group
+    """
+
+    emins = []
+    emaxs = []
+    fmins = []
+    fmaxs = []
+
+    for filename in snaplist:
+
+        data = swiftsimio.load(filename)
+        meta = data.metadata
+
+        ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"])
+        emin_group = []
+        emax_group = []
+        fluxmin_group = []
+        fluxmax_group = []
+
+        for g in range(ngroups):
+            en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+            emin_group.append((1 * en.min()).value)
+            emax_group.append((1 * en.max()).value)
+
+            for direction in ["X"]:
+                #  for direction in ["X", "Y", "Z"]:
+                f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+                fluxmin_group.append(
+                    (1 * f.min()).to(unyt.erg / unyt.cm ** 2 / unyt.s).value
+                )
+                fluxmax_group.append(
+                    (1 * f.max()).to(unyt.erg / unyt.cm ** 2 / unyt.s).value
+                )
+
+        emins.append(emin_group)
+        emaxs.append(emax_group)
+        fmins.append(fluxmin_group)
+        fmaxs.append(fluxmax_group)
+
+    energy_boundaries = []
+    flux_boundaries = []
+    for g in range(ngroups):
+        emin = min([emins[f][g] for f in range(len(snaplist))])
+        emax = max([emaxs[f][g] for f in range(len(snaplist))])
+        energy_boundaries.append([emin, emax])
+        fmin = min([fmins[f][g] for f in range(len(snaplist))])
+        fmax = max([fmaxs[f][g] for f in range(len(snaplist))])
+        flux_boundaries.append([fmin, fmax])
+
+    return energy_boundaries, flux_boundaries
+
+
+if __name__ == "__main__":
+    # Get command line arguments
+    args = parse_args()
+
+    if args.snapshot_number:
+        plot_all = False
+        snapnr = int(args.snapshot_number)
+    else:
+        plot_all = True
+
+    domain = args.redshift
+    snaplist = get_snapshot_list(snapshot_base + f"_{domain}")
+
+    if fancy:
+        energy_boundaries, flux_boundaries = get_minmax_vals(snaplist)
+    else:
+        energy_boundaries, flux_boundaries = (None, None)
+
+    for f in snaplist:
+        plot_photons(
+            f, energy_boundaries=energy_boundaries, flux_boundaries=flux_boundaries
+        )
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_high_redshift.yml b/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_high_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e99c1370c6f4d9988e8f08c4e4958f39b077983f
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_high_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: RT_advection-1D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # kpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_high_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_high_redshift
+  scale_factor_first:  0.00990099  # Time of the first output (in internal units)
+  delta_time:          1.06
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.00990099
+  delta_time:          1.06   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advection_1D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [0., 1., 2.]  # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [0., 0., 0.]  # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                  # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                       # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  skip_thermochemistry: 1                        # ignore thermochemistry.
+  set_initial_ionization_mass_fractions: 1       # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                        # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                        # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.00990099  # z~100
+  a_end:          0.01408451  # z~70
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_low_redshift.yml b/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_low_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..14da127f741675fbea8ad7825ec11e24b9df4885
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_low_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: RT_advection-1D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_low_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_low_redshift
+  scale_factor_first:  0.93  # Time of the first output (in internal units)
+  delta_time:          1.005
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.93
+  delta_time:          1.005   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advection_1D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [0., 1., 2.]  # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [0., 0., 0.]  # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                  # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                       # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  skip_thermochemistry: 1                        # ignore thermochemistry.
+  set_initial_ionization_mass_fractions: 1       # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                        # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                        # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.93
+  a_end:          1.
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_medium_redshift.yml b/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_medium_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9b5f742fd762cec2877f59f9c4b8acad6608b71d
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/rt_advection1D_medium_redshift.yml
@@ -0,0 +1,71 @@
+MetaData:
+  run_name: RT_advection-1D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2   # The maximal time-step size of the simulation (in internal units).
+  time_begin: 0
+  time_end:   1.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_medium_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_medium_redshift
+  scale_factor_first:  0.0909
+  delta_time:          1.03
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.0909
+  delta_time:          1.02   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advection_1D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [0., 1., 2.]  # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [0., 0., 0.]  # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                  # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                       # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17 # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  skip_thermochemistry: 1                        # ignore thermochemistry.
+  set_initial_ionization_mass_fractions: 1       # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                        # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                        # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.0909    # z=10
+  a_end:          0.1428571 # z=6
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/run.sh b/examples/RadiativeTransferTests/CosmoAdvection_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..964d3e80e2c665b45228f98fe3552bbd4a380e04
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/run.sh
@@ -0,0 +1,52 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -f advection_1D.hdf5 ]; then
+    echo "Generating ICs"
+    python3 makeIC.py
+fi
+
+# Default run
+ymlfile=rt_advection1D_high_redshift.yml
+zdomain="h"
+
+# Do we have a cmdline argument provided?
+if [ $# -gt 0 ]; then
+    case "$1" in 
+    -l | -low | --l | --low | l | ./rt_advection1D_low_redshift.yml | rt_advection1D_low_redshift | rt_advection1D_low_redshift.yml )
+        ymlfile=rt_advection1D_low_redshift.yml
+	zdomain="l"
+        ;;
+    -m | -mid | --m | --mid | m | ./rt_advection1D_medium_redshift.yml | rt_advection1D_medium_redshift | rt_advection1D_medium_redshift.yml )
+        ymlfile=rt_advection1D_medium_redshift.yml
+	zdomain="m"
+        ;;
+    -h | -high | -hi | --h | --hi | --high | h | ./rt_advection1D_high_redshift.yml | rt_advection1D_high_redshift | rt_advection1D_high_redshift.yml )
+        ymlfile=rt_advection1D_high_redshift.yml
+	zdomain="h"
+        ;;
+    *)
+        echo unknown cmdline param, running default $ymlfile
+        ;;
+    esac
+fi
+
+
+
+# Run SWIFT with RT
+../../../swift \
+    --hydro \
+    --cosmology \
+    --threads=4 \
+    --verbose=0  \
+    --radiation \
+    --stars \
+    --feedback \
+    --external-gravity \
+    $ymlfile 2>&1 | tee output.log
+
+python3 ./plotSolution.py -z $zdomain
+python3 ./plotEnergy.py -z $zdomain
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_high_redshift b/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_high_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e504014177837a7b2c1d134d1c4a0d5a785d0044
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_high_redshift
@@ -0,0 +1,12 @@
+# Redshift
+99.98460572382712
+96.47031888378449
+93.23367859034735
+90.25299848922855
+87.49735694548497
+84.9407463291136
+82.56112036003019
+80.33965676782498
+78.26018022924342
+76.30870609169932
+74.47307621234665
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_low_redshift b/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_low_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e58febac09ba12e5677c6fc8abafe94cc31824ad
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_low_redshift
@@ -0,0 +1,12 @@
+# Redshift
+0.0751932785953775
+0.06869649828543634
+0.06225945012746337
+0.055881161353223074
+0.04956068264358926
+0.04329708741463545
+0.03708947113706218
+0.030936950674366415
+0.02483866364616727
+0.018793767817153695
+0.01280144050017884
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_medium_redshift b/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_medium_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..958165894af540df676443cfdf3fbf20ce46dd1c
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_1D/snapshot_times_medium_redshift
@@ -0,0 +1,12 @@
+# Redshift
+9.983090485639307
+9.499283920730043
+8.92297287774699
+8.419654113369706
+7.975570397565402
+7.580294600085116
+7.225767198868688
+6.905651913777566
+6.614892571630509
+6.349401127100822
+6.1058334302500645
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/README b/examples/RadiativeTransferTests/CosmoAdvection_2D/README
new file mode 100644
index 0000000000000000000000000000000000000000..d08a30340bc6122bc675fba4e927365ae5c5378d
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/README
@@ -0,0 +1,26 @@
+2D advection test for radiative transfer.
+
+Test that your method is TVD and the propagation speed of the photons is
+correct. The ICs set up four photon groups: 
+- The first is a top hat function initial distribution where outside values
+  are zero, advecting along the x direction
+- The second is a top hat function initial distribution where outside values
+  are nonzero. This distinction is important to test because photon energies 
+  can't be negative, so these cases need to be tested individually. This
+  group advects along the y direction
+- the third is a smooth Gaussian advecting diagonally.
+- the fourth is a circle in the center advecting radially.
+
+This way, you can test multiple initial condition scenarios simultaneously. 
+There are no stars to act as sources. Also make sure that you choose your
+photon frequencies in a way that doesn't interact with gas!
+
+The ICs are created to be compatible with GEAR_RT. Recommended configuration:
+    --with-rt=GEAR_4 --with-rt-riemann-solver=GLF --with-hydro-dimension=2 --with-hydro=gizmo-mfv \
+     --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
+
+    
+
+Note that if you want to use a reduced speed of light for this test, you also 
+need to adapt the fluxes in the initial conditions! They are generated assuming
+that the speed of light is not reduced.
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/getGlass.sh b/examples/RadiativeTransferTests/CosmoAdvection_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ae3c977064f5e7a408aa249c5fd9089b3c52ecb1
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/makeIC.py b/examples/RadiativeTransferTests/CosmoAdvection_2D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..d3769ab5498cc445e7e5822836a2cc6df8de3c7f
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/makeIC.py
@@ -0,0 +1,220 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+#               2022 Tsang Keung Chan (chantsangkeung@gmail.com)
+#               2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+
+# -------------------------------------------------------------
+# Add initial conditions for photon energies and fluxes
+# for 2D advection of photons.
+# First photon group: Top hat function with zero as the
+#       baseline, advects in x direction
+# Second photon group: Top hat function with nonzero value
+#       as the baseline, advcts in y direction.
+# Third photon group: Gaussian advecting diagonally
+# Fourth photon group: Circle moving radially from the center
+# -------------------------------------------------------------
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+# define unit system to use
+unitsystem = cosmo_units
+
+# define box size
+boxsize = 260 * unyt.Mpc
+boxsize = boxsize.to(unitsystem["length"])
+
+reduced_speed_of_light_fraction = 1.0
+
+# number of photon groups
+nPhotonGroups = 4
+
+# filename of ICs to be generated
+outputfilename = "advection_2D.hdf5"
+
+
+def initial_condition(x, unitsystem):
+    """
+    The initial conditions that will be advected
+
+    x: particle position. 3D unyt array
+
+    returns: 
+    E: photon energy for each photon group. List of scalars with size of nPhotonGroups
+    F: photon flux for each photon group. List with size of nPhotonGroups of numpy arrays of shape (3,)
+    """
+
+    # you can make the photon quantities unitless, the units will
+    # already have been written down in the writer.
+    # However, that means that you need to convert them manually.
+
+    unit_energy = (
+        unitsystem["mass"] * unitsystem["length"] ** 2 / unitsystem["time"] ** 2
+    )
+    unit_velocity = unitsystem["length"] / unitsystem["time"]
+    unit_flux = unit_energy * unit_velocity
+
+    c_internal = (unyt.c * reduced_speed_of_light_fraction).to(unit_velocity)
+
+    E_list = []
+    F_list = []
+
+    # Group 1 Photons:
+    # -------------------
+
+    in_x = 0.33 * boxsize < x[0] < 0.66 * boxsize
+    in_y = 0.33 * boxsize < x[1] < 0.66 * boxsize
+    if in_x and in_y:
+        E = 1.0 * unit_energy
+    else:
+        E = 0.0 * unit_energy
+
+    # Assuming all photons flow in only one direction
+    # (optically thin regime, "free streaming limit"),
+    #  we have that |F| = c * E
+    F = np.zeros(3, dtype=np.float64)
+    F[0] = (c_internal * E).to(unit_flux)
+
+    E_list.append(E)
+    F_list.append(F)
+
+    # Group 2 Photons:
+    # -------------------
+
+    in_x = 0.33 * boxsize < x[0] < 0.66 * boxsize
+    in_y = 0.33 * boxsize < x[1] < 0.66 * boxsize
+    if in_x and in_y:
+        E = 2.0 * unit_energy
+    else:
+        E = 1.0 * unit_energy
+
+    F = np.zeros(3, dtype=np.float64)
+    F[1] = (c_internal * E).to(unit_flux)
+
+    E_list.append(E)
+    F_list.append(F)
+
+    # Group 3 Photons:
+    # -------------------
+    sigma = 0.1 * boxsize
+    mean = 0.5 * boxsize
+    amplitude = 2.0
+    baseline = 1.0
+
+    E = (
+        amplitude
+        * np.exp(-((x[0] - mean) ** 2 + (x[1] - mean) ** 2) / (2 * sigma ** 2))
+        + baseline
+    ) * unit_energy
+
+    F = np.zeros(3, dtype=np.float64)
+    F[0] = (c_internal * E / 1.414213562).to(unit_flux)  # sqrt(2)
+    F[1] = (c_internal * E / 1.414213562).to(unit_flux)  # sqrt(2)
+
+    E_list.append(E)
+    F_list.append(F)
+
+    # Group 4 Photons:
+    # -------------------
+
+    circle_radius = 0.15 * boxsize
+    center = 0.5 * boxsize
+    dx = x[0] - center
+    dy = x[1] - center
+    r = np.sqrt(dx ** 2 + dy ** 2)
+    if r <= circle_radius:
+        unit_vector = (dx / r, dy / r)
+
+        E = 1.0 * unit_energy
+        F = np.zeros(3, dtype=np.float64)
+        F[0] = (unit_vector[0] * c_internal * E).to(unit_flux)
+        F[1] = (unit_vector[1] * c_internal * E).to(unit_flux)
+
+    else:
+        E = 0.0 * unit_energy
+        F = np.zeros(3, dtype=np.float64)
+
+    E_list.append(E)
+    F_list.append(F)
+
+    return E_list, F_list
+
+
+if __name__ == "__main__":
+    glass = h5py.File("glassPlane_128.hdf5", "r")
+
+    # Read particle positions and h from the glass
+    pos = glass["/PartType0/Coordinates"][:, :]
+    h = glass["/PartType0/SmoothingLength"][:]
+    glass.close()
+
+    pos *= boxsize
+    h *= boxsize
+    numPart = np.size(h)
+
+    w = Writer(unitsystem, boxsize, dimension=2)
+
+    w.gas.coordinates = pos
+    w.gas.velocities = np.zeros((numPart, 3)) * (unyt.cm / unyt.s)
+    mpart = 1e20 * unyt.M_Sun
+    mpart = mpart.to(unitsystem["mass"])
+    w.gas.masses = np.ones(numPart, dtype=np.float64) * mpart
+    w.gas.internal_energy = (
+        np.ones(numPart, dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g
+    )
+
+    # Generate initial guess for smoothing lengths based on MIPS
+    w.gas.smoothing_length = h
+
+    # If IDs are not present, this automatically generates
+    w.write(outputfilename)
+
+    # Now open file back up again and add photon groups
+    # you can make them unitless, the units have already been
+    # written down in the writer. In this case, it's in cgs.
+
+    F = h5py.File(outputfilename, "r+")
+    header = F["Header"]
+    nparts = header.attrs["NumPart_ThisFile"][0]
+    parts = F["/PartType0"]
+
+    for grp in range(nPhotonGroups):
+        dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1)
+        energydata = np.zeros(nparts, dtype=np.float32)
+        parts.create_dataset(dsetname, data=energydata)
+
+        dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1)
+        #  if dsetname not in parts.keys():
+        fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+        parts.create_dataset(dsetname, data=fluxdata)
+
+    for p in range(nparts):
+        E, Flux = initial_condition(pos[p], unitsystem)
+        for g in range(nPhotonGroups):
+            Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1)
+            parts[Esetname][p] = E[g]
+            Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1)
+            parts[Fsetname][p] = Flux[g]
+
+    F.close()
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/plotEnergy.py b/examples/RadiativeTransferTests/CosmoAdvection_2D/plotEnergy.py
new file mode 100755
index 0000000000000000000000000000000000000000..c32ad26515e698d401003621cc47c4020c95b3d7
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/plotEnergy.py
@@ -0,0 +1,205 @@
+#!/usr/bin/env python3
+
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+import os
+import sys
+import argparse
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+import unyt
+from matplotlib import pyplot as plt
+
+# Parameters users should/may tweak
+snapshot_base = "output"  # Snapshot basename
+plot_physical_quantities = True  # Plot physical or comoving quantities
+
+time_first = 0  # Time of first snapshot
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-z", "--redshift", help="Redshift domain to plot advection for", default="high"
+    )
+
+    args = parser.parse_args()
+    return args
+
+
+def get_snapshot_list(snapshot_basename="output"):
+    """
+    Find the snapshot(s) that are to be plotted
+    and return their names as list
+    """
+
+    snaplist = []
+
+    dirlist = os.listdir()
+    for f in dirlist:
+        if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+            snaplist.append(f)
+
+    snaplist = sorted(snaplist)
+    if len(snaplist) == 0:
+        print(f"No snapshots with base {snapshot_basename} found!")
+        sys.exit(1)
+    return snaplist
+
+
+def plot_param_over_time(snapshot_list, param="energy density", redshift_domain="high"):
+    print(f"Now plotting {param} over time")
+    # Grab number of photon groups
+    data = swiftsimio.load(snapshot_list[0])
+    meta = data.metadata
+    ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"][0])
+
+    # Number of rows and columns
+    nrows = 1 + int(plot_physical_quantities)
+    ncols = ngroups
+    # Create figure and axes
+    fig, axs = plt.subplots(nrows, ncols, figsize=(5.04 * ncols, 5.4 * nrows), dpi=200)
+
+    # Iterate over all photon groups
+    for n in range(ngroups):
+        # Arrays to keep track of plot_param and scale factor
+        plot_param = [[], []]
+        scale_factor = []
+        analytic_exponent = [0, 0]
+
+        # Functions to convert between scale factor and redshift
+        a2z = lambda a: 1 / a - 1
+        z2a = lambda z: 1 / (z + 1)
+
+        for file in snapshot_list:
+            data = swiftsimio.load(file)
+            meta = data.metadata
+
+            # Read comoving variables
+            energy = getattr(data.gas.photon_energies, f"group{n+1}")
+            mass = data.gas.masses
+            rho = data.gas.densities
+            volume = mass / rho
+
+            energy_density = energy / volume
+
+            if plot_physical_quantities:
+                # The SWIFT cosmology module assumes 3-dimensional lengths and volumes,
+                # so multiply by a to get the correct relations
+                physical_energy_density = (
+                    energy_density.to_physical() * meta.scale_factor
+                )
+                physical_energy = energy.to_physical()
+
+                if param == "energy density":
+                    plot_param[1].append(
+                        1
+                        * physical_energy_density.sum()
+                        / physical_energy_density.shape[0]
+                    )
+                    analytic_exponent[1] = -2.0
+                elif param == "total energy":
+                    plot_param[1].append(1 * physical_energy.sum())
+                    analytic_exponent[1] = 0.0
+
+            if param == "energy density":
+                plot_param[0].append(1 * energy_density.sum() / energy_density.shape[0])
+                analytic_exponent[0] = 0.0
+            elif param == "total energy":
+                plot_param[0].append(1 * energy.sum())
+                analytic_exponent[0] = 0.0
+
+            scale_factor.append(meta.scale_factor)
+
+        if param == "energy density":
+            titles = ["Comoving energy density", "Physical energy density $\\times a$"]
+            ylabel = "Average energy density"
+            figname = "output_energy_density_over_time"
+        elif param == "total energy":
+            titles = ["Comoving total energy", "Physical total energy"]
+            ylabel = "Total energy"
+            figname = "output_total_energy_over_time"
+
+        # Analytic scale factor
+        analytic_scale_factor = np.linspace(min(scale_factor), max(scale_factor), 1000)
+
+        for i in range(nrows):
+            ax = axs[i, n]
+            ax.scatter(scale_factor, plot_param[i], label="Simulation")
+
+            # Analytic scale factor relation
+            analytic = analytic_scale_factor ** analytic_exponent[i]
+
+            # Scale solution to correct offset
+            analytic = analytic / analytic[0] * plot_param[i][0]
+            ax.plot(
+                analytic_scale_factor,
+                analytic,
+                c="r",
+                label=f"Analytic solution $\propto a^{{{analytic_exponent[i]}}}$",
+            )
+
+            ax.legend()
+            ax.set_title(titles[i] + f" group {n+1}")
+
+            ax.set_xlabel("Scale factor")
+            secax = ax.secondary_xaxis("top", functions=(a2z, z2a))
+            secax.set_xlabel("Redshift")
+
+            ax.yaxis.get_offset_text().set_position((-0.05, 1))
+
+            if analytic_exponent[i] == 0.0:
+                ax.set_ylim(plot_param[i][0] * 0.95, plot_param[i][0] * 1.05)
+                ylabel_scale = ""
+            else:
+                ylabel_scale = "$\\times a$"
+            if n == 0:
+                units = plot_param[i][0].units.latex_representation()
+                ax.set_ylabel(f"{ylabel} [${units}$] {ylabel_scale}")
+    plt.tight_layout()
+    plt.savefig(f"{figname}-{redshift_domain}.png")
+    plt.close()
+
+
+if __name__ in ("__main__"):
+    # Get command line args
+    args = parse_args()
+    domain = args.redshift.lower()
+    if domain in ("low", "l", "low_redshift", "low redshift", "low-redshift"):
+        redshift_domain = "low_redshift"
+    elif domain in (
+        "medium",
+        "m",
+        "medium_redshift",
+        "medium redshift",
+        "medium-redshift",
+    ):
+        redshift_domain = "medium_redshift"
+    elif domain in ("high", "h", "high_redshift", "high redshift", "high-redshift"):
+        redshift_domain = "high_redshift"
+    else:
+        print("Redshift domain not recognised!")
+        sys.exit(1)
+
+    snaplist = get_snapshot_list(snapshot_base + f"_{redshift_domain}")
+
+    for param in ["energy density", "total energy"]:
+        plot_param_over_time(snaplist, param, redshift_domain)
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/plotSolution.py b/examples/RadiativeTransferTests/CosmoAdvection_2D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..6b7eb6a3918937cdc29badf72973f7fa51b2e18e
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/plotSolution.py
@@ -0,0 +1,354 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+#               2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+
+# ----------------------------------------------------
+# plot photon data for 2D problems
+# give snapshot number as cmdline arg to plot
+# single snapshot, otherwise this script plots
+# all snapshots available in the workdir
+# ----------------------------------------------------
+
+import gc
+import os
+import sys
+import argparse
+
+import unyt
+import numpy as np
+import matplotlib as mpl
+import swiftsimio
+from matplotlib import pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+
+# Parameters users should/may tweak
+plot_all_data = True  # plot all groups and all photon quantities
+snapshot_base = "output"  # snapshot basename
+fancy = True  # fancy up the plots a bit?
+plot_physical_quantities = True
+
+# parameters for imshow plots
+imshow_kwargs = {"origin": "lower", "cmap": "viridis"}
+
+
+projection_kwargs = {"resolution": 1024, "parallel": True}
+# -----------------------------------------------------------------------
+
+
+plot_all = False
+mpl.rcParams["text.usetex"] = True
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-n", "--snapshot-number", help="Number of snapshot to plot", type=int
+    )
+    parser.add_argument(
+        "-z",
+        "--redshift",
+        help="Redshift domain to plot advection for",
+        default="high_redshift",
+    )
+
+    args = parser.parse_args()
+    return args
+
+
+def get_snapshot_list(snapshot_basename="output"):
+    """
+    Find the snapshot(s) that are to be plotted 
+    and return their names as list
+    """
+
+    snaplist = []
+
+    if plot_all:
+        dirlist = os.listdir()
+        for f in dirlist:
+            if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+                snaplist.append(f)
+
+        snaplist = sorted(snaplist)
+        if len(snaplist) == 0:
+            print(f"No snapshots with base {snapshot_basename} found!")
+            sys.exit(1)
+
+    else:
+        fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+        if not os.path.exists(fname):
+            print("Didn't find file", fname)
+            quit(1)
+        snaplist.append(fname)
+
+    return snaplist
+
+
+def set_colorbar(ax, im):
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size="5%", pad=0.05)
+    plt.colorbar(im, cax=cax)
+    return
+
+
+def plot_photons(filename, energy_boundaries=None, flux_boundaries=None):
+    """
+    Create the actual plot.
+
+    filename: file to work with
+    energy_boundaries:  list of [E_min, E_max] for each photon group. 
+                        If none, limits are set automatically.
+    flux_boundaries:    list of [F_min, F_max] for each photon group. 
+                        If none, limits are set automatically.
+    """
+
+    print("working on", filename)
+
+    # Read in data first
+    data = swiftsimio.load(filename)
+    meta = data.metadata
+
+    ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"][0])
+    xlabel_units_str = meta.boxsize.units.latex_representation()
+
+    global imshow_kwargs
+    imshow_kwargs["extent"] = [0, meta.boxsize[0].v, 0, meta.boxsize[1].v]
+
+    for g in range(ngroups):
+        # workaround to access named columns data with swiftsimio visualisaiton
+        # add mass weights to remove surface density dependence in images
+        new_attribute_str = "mass_weighted_radiation_energy" + str(g + 1)
+        en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+        en = en * data.gas.masses
+        setattr(data.gas, new_attribute_str, en)
+
+        if plot_all_data:
+            # prepare also the fluxes
+            #  for direction in ["X", "Y", "Z"]:
+            for direction in ["X", "Y"]:
+                new_attribute_str = (
+                    "mass_weighted_radiation_flux" + str(g + 1) + direction
+                )
+                f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+                # f *= data.gas.masses
+                f = f * data.gas.masses
+                setattr(data.gas, new_attribute_str, f)
+
+    # get mass surface density projection that we'll use to remove density dependence in image
+    mass_map = swiftsimio.visualisation.projection.project_gas(
+        data, project="masses", **projection_kwargs
+    )
+
+    if plot_all_data:
+        fig = plt.figure(figsize=(5 * 3, 5.05 * ngroups), dpi=200)
+        figname = filename[:-5] + f"-all-quantities.png"
+
+        for g in range(ngroups):
+
+            # get energy projection
+            new_attribute_str = "mass_weighted_radiation_energy" + str(g + 1)
+            photon_map = swiftsimio.visualisation.projection.project_gas(
+                data, project=new_attribute_str, **projection_kwargs
+            )
+            photon_map = photon_map / mass_map
+
+            ax = fig.add_subplot(ngroups, 3, g * 3 + 1)
+            if energy_boundaries is not None:
+                imshow_kwargs["vmin"] = energy_boundaries[g][0]
+                imshow_kwargs["vmax"] = energy_boundaries[g][1]
+            im = ax.imshow(photon_map.T, **imshow_kwargs)
+            set_colorbar(ax, im)
+            ax.set_ylabel("Group {0:2d}".format(g + 1))
+            ax.set_xlabel("x [$" + xlabel_units_str + "$]")
+            if g == 0:
+                ax.set_title("Energies")
+
+            # get flux X projection
+            new_attribute_str = "mass_weighted_radiation_flux" + str(g + 1) + "X"
+            photon_map = swiftsimio.visualisation.projection.project_gas(
+                data, project=new_attribute_str, **projection_kwargs
+            )
+            photon_map = photon_map / mass_map
+
+            ax = fig.add_subplot(ngroups, 3, g * 3 + 2)
+            if flux_boundaries is not None:
+                imshow_kwargs["vmin"] = flux_boundaries[g][0]
+                imshow_kwargs["vmax"] = flux_boundaries[g][1]
+            im = ax.imshow(photon_map.T, **imshow_kwargs)
+            set_colorbar(ax, im)
+            ax.set_xlabel("x [$" + xlabel_units_str + "$]")
+            ax.set_ylabel("y [$" + xlabel_units_str + "$]")
+            if g == 0:
+                ax.set_title("Flux X")
+
+            # get flux Y projection
+            new_attribute_str = "mass_weighted_radiation_flux" + str(g + 1) + "Y"
+            photon_map = swiftsimio.visualisation.projection.project_gas(
+                data, project=new_attribute_str, **projection_kwargs
+            )
+            photon_map = photon_map / mass_map
+
+            ax = fig.add_subplot(ngroups, 3, g * 3 + 3)
+            im = ax.imshow(photon_map.T, **imshow_kwargs)
+            set_colorbar(ax, im)
+            ax.set_xlabel("x [$" + xlabel_units_str + "$]")
+            ax.set_ylabel("y [$" + xlabel_units_str + "$]")
+            if g == 0:
+                ax.set_title("Flux Y")
+
+    else:  # plot just energies
+
+        fig = plt.figure(figsize=(5 * ngroups, 5), dpi=200)
+        figname = filename[:-5] + ".png"
+
+        for g in range(ngroups):
+
+            # get projection
+            new_attribute_str = "mass_weighted_radiation_energy" + str(g + 1)
+            photon_map = swiftsimio.visualisation.projection.project_gas(
+                data, project=new_attribute_str, **projection_kwargs
+            )
+            photon_map = photon_map / mass_map
+
+            ax = fig.add_subplot(1, ngroups, g + 1)
+            if energy_boundaries is not None:
+                imshow_kwargs["vmin"] = energy_boundaries[g][0]
+                imshow_kwargs["vmax"] = energy_boundaries[g][1]
+            im = ax.imshow(photon_map.T, **imshow_kwargs)
+            set_colorbar(ax, im)
+            ax.set_title("Group {0:2d}".format(g + 1))
+            if g == 0:
+                ax.set_ylabel("Energies")
+
+    # Add title
+    title = filename.replace("_", r"\_")  # exception handle underscore for latex
+    if meta.cosmology is not None:
+        title += ", $z$ = {0:.2e}".format(meta.z)
+    title += ", $t$ = {0:.2e}".format(1 * meta.time)
+    fig.suptitle(title)
+
+    plt.tight_layout()
+    plt.savefig(figname)
+    plt.close()
+    gc.collect()
+
+    return
+
+
+def get_minmax_vals(snaplist):
+    """
+    Find minimal and maximal values for energy and flux,
+    so you can fix axes limits over all snapshots
+
+    snaplist: list of snapshot filenames
+
+    returns:
+
+    energy_boundaries: list of [E_min, E_max] for each photon group
+    flux_boundaries: list of [Fx_min, Fy_max] for each photon group
+    """
+
+    emins = []
+    emaxs = []
+    fmins = []
+    fmaxs = []
+
+    for filename in snaplist:
+
+        data = swiftsimio.load(filename)
+        meta = data.metadata
+
+        ngroups = int(meta.subgrid_scheme["PhotonGroupNumber"][0])
+        emin_group = []
+        emax_group = []
+        fluxmin_group = []
+        fluxmax_group = []
+
+        for g in range(ngroups):
+            en = getattr(data.gas.photon_energies, "group" + str(g + 1))
+            emin_group.append((1 * en.min()).value)
+            emax_group.append((1 * en.max()).value)
+
+            dirmin = []
+            dirmax = []
+            for direction in ["X", "Y"]:
+                f = getattr(data.gas.photon_fluxes, "Group" + str(g + 1) + direction)
+                dirmin.append((1 * f.min()).value)
+                dirmax.append((1 * f.max()).value)
+            fluxmin_group.append(min(dirmin))
+            fluxmax_group.append(max(dirmax))
+
+        emins.append(emin_group)
+        emaxs.append(emax_group)
+        fmins.append(fluxmin_group)
+        fmaxs.append(fluxmax_group)
+
+    energy_boundaries = []
+    flux_boundaries = []
+    for g in range(ngroups):
+        emin = min([emins[f][g] for f in range(len(snaplist))])
+        emax = max([emaxs[f][g] for f in range(len(snaplist))])
+        energy_boundaries.append([emin, emax])
+        fmin = min([fmins[f][g] for f in range(len(snaplist))])
+        fmax = max([fmaxs[f][g] for f in range(len(snaplist))])
+        flux_boundaries.append([fmin, fmax])
+
+    return energy_boundaries, flux_boundaries
+
+
+if __name__ == "__main__":
+    # Get command line arguments
+    args = parse_args()
+
+    if args.snapshot_number:
+        plot_all = False
+        snapnr = int(args.snapshot_number)
+    else:
+        plot_all = True
+
+    domain = args.redshift
+    if domain in ("low", "l", "low_redshift", "low redshift", "low-redshift"):
+        redshift_domain = "low_redshift"
+    elif domain in (
+        "medium",
+        "m",
+        "medium_redshift",
+        "medium redshift",
+        "medium-redshift",
+    ):
+        redshift_domain = "medium_redshift"
+    elif domain in ("high", "h", "high_redshift", "high redshift", "high-redshift"):
+        redshift_domain = "high_redshift"
+    else:
+        print("Redshift domain not recognised!")
+        sys.exit(1)
+
+    snaplist = get_snapshot_list(snapshot_base + f"_{domain}")
+    if fancy:
+        energy_boundaries, flux_boundaries = get_minmax_vals(snaplist)
+    else:
+        energy_boundaries = None
+        flux_boundaries = None
+
+    for f in snaplist:
+        plot_photons(
+            f, energy_boundaries=energy_boundaries, flux_boundaries=flux_boundaries
+        )
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_high_redshift.yml b/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_high_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..1e257683a84c1379eb623300e06c592360701b39
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_high_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: cosmo_RT_advection-2D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # kpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_high_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_high_redshift
+  scale_factor_first:  0.00990099  # Time of the first output (in internal units)
+  delta_time:          1.06
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.00990099
+  delta_time:          1.06   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advection_2D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0                       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99                             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [1., 2., 3., 4.]              # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                 # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [1e-32, 1e-32, 1e-32, 1e-32]   # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                   # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                        # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17  # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1        # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                         # Skip thermochemsitry.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.00990099  # z~100
+  a_end:          0.01408451  # z~70
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_low_redshift.yml b/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_low_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..64db5443d404607066faea7513f188dd7cb6fa03
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_low_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: cosmo_RT_advection-2D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_low_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_low_redshift
+  scale_factor_first:  0.93  # Time of the first output (in internal units)
+  delta_time:          1.005
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.93
+  delta_time:          1.005   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advection_2D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0                       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99                             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [1., 2., 3., 4.]              # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                 # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [1e-32, 1e-32, 1e-32, 1e-32]   # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                   # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                        # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17  # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1        # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                         # Skip thermochemsitry.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.93
+  a_end:          1.
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_medium_redshift.yml b/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_medium_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e4cc8b2d0f227cad0e8211ab5833a129e28962fe
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/rt_advection2D_medium_redshift.yml
@@ -0,0 +1,71 @@
+MetaData:
+  run_name: cosmo_RT_advection-2D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2   # The maximal time-step size of the simulation (in internal units).
+  time_begin: 0
+  time_end:   1.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_medium_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_medium_redshift
+  scale_factor_first:  0.0909
+  delta_time:          1.01
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.0909
+  delta_time:          1.02   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./advection_2D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0                       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99                             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [1., 2., 3., 4.]              # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                 # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [1e-32, 1e-32, 1e-32, 1e-32]   # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                   # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                        # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17  # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1        # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                         # Skip thermochemsitry.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.0909    # z=10
+  a_end:          0.1428571 # z=6
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/run.sh b/examples/RadiativeTransferTests/CosmoAdvection_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c0c7bacab484cabd06fa0c2eefbce8357c4e74ad
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/run.sh
@@ -0,0 +1,57 @@
+#!/bin/bash
+
+# exit if anything fails
+set -e
+set -o pipefail
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+    echo "Fetching initial glass file for the 2D RT advection example..."
+    ./getGlass.sh
+fi
+if [ ! -e advection_2D.hdf5 ]
+then
+    echo "Generating initial conditions for the 2D RT advection example..."
+    python3 makeIC.py
+fi
+
+# Default run
+ymlfile=rt_advection2D_high_redshift.yml
+zdomain="h"
+
+# Do we have a cmdline argument provided?
+if [ $# -gt 0 ]; then
+    case "$1" in
+    -l | -low | --l | --low | l | ./rt_advection2D_low_redshift.yml | rt_advection2D_low_redshift | rt_advection2D_low_redshift.yml )
+        ymlfile=rt_advection2D_low_redshift.yml
+	zdomain="l"
+        ;;
+    -m | -mid | --m | --mid | m | ./rt_advection2D_medium_redshift.yml | rt_advection2D_medium_redshift | rt_advection2D_medium_redshift.yml )
+        ymlfile=rt_advection2D_medium_redshift.yml
+	zdomain="m"
+        ;;
+    -h | -high | -hi | --h | --hi | --high | h | ./rt_advection2D_high_redshift.yml | rt_advection2D_high_redshift | rt_advection2D_high_redshift.yml )
+        ymlfile=rt_advection2D_high_redshift.yml
+	zdomain="h"
+        ;;
+    *)
+        echo unknown cmdline param, running default $ymlfile
+        ;;
+    esac
+fi
+
+# Run SWIFT with RT
+../../../swift \
+    --hydro \
+    --cosmology \
+    --threads=4 \
+    --verbose=0  \
+    --radiation \
+    --stars \
+    --feedback \
+    --external-gravity \
+    $ymlfile 2>&1 | tee output.log
+
+python3 ./plotSolution.py -z $zdomain
+python3 ./plotEnergy.py -z $zdomain
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_high_redshift b/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_high_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e504014177837a7b2c1d134d1c4a0d5a785d0044
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_high_redshift
@@ -0,0 +1,12 @@
+# Redshift
+99.98460572382712
+96.47031888378449
+93.23367859034735
+90.25299848922855
+87.49735694548497
+84.9407463291136
+82.56112036003019
+80.33965676782498
+78.26018022924342
+76.30870609169932
+74.47307621234665
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_low_redshift b/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_low_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e58febac09ba12e5677c6fc8abafe94cc31824ad
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_low_redshift
@@ -0,0 +1,12 @@
+# Redshift
+0.0751932785953775
+0.06869649828543634
+0.06225945012746337
+0.055881161353223074
+0.04956068264358926
+0.04329708741463545
+0.03708947113706218
+0.030936950674366415
+0.02483866364616727
+0.018793767817153695
+0.01280144050017884
diff --git a/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_medium_redshift b/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_medium_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..958165894af540df676443cfdf3fbf20ce46dd1c
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoAdvection_2D/snapshot_times_medium_redshift
@@ -0,0 +1,12 @@
+# Redshift
+9.983090485639307
+9.499283920730043
+8.92297287774699
+8.419654113369706
+7.975570397565402
+7.580294600085116
+7.225767198868688
+6.905651913777566
+6.614892571630509
+6.349401127100822
+6.1058334302500645
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/README b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/README
new file mode 100644
index 0000000000000000000000000000000000000000..5ce885e6fe84b375e91521fa6899a91881e52bc7
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/README
@@ -0,0 +1,14 @@
+Strömgen Sphere example in 3D
+-----------------------------
+
+This directory contains the following example:
+    -   run a propagation test of photons emitted from a single 
+        central source in an otherwise uniform box.
+        To run this example, use the provided `runPropagationTest.sh` script.
+        This script will then make use of the `makePropagationTestIC.py` and
+        `plotPhotonPropagationCheck.py` script.
+             
+
+To use the GEAR RT model, compile with :
+    --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro=gizmo-mfv --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
+
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/getGlass.sh b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..50f27f4e9f9981da9449608bd70669c56d9d3985
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/getGlass.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
+#wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_128.hdf5
\ No newline at end of file
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/makePropagationTestIC.py b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/makePropagationTestIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..49e24b65e7d74649fc151dab17169b030b629ad2
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/makePropagationTestIC.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+#               2022 Tsang Keung Chan (chantsangkeung@gmail.com)
+#               2024 Stan Verhoeve (s06verhoeve@gmail.com)
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+
+# ---------------------------------------------------------------------
+# Add a single star in the center of a glass distribution
+# Intended for the photon propagation test.
+# ---------------------------------------------------------------------
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+glass = h5py.File("glassCube_64.hdf5", "r")
+parts = glass["PartType0"]
+xp = parts["Coordinates"][:]
+h = parts["SmoothingLength"][:]
+glass.close()
+
+# replace the particle closest to the center
+# by the star
+r = np.sqrt(np.sum((0.5 - xp) ** 2, axis=1))
+rmin = np.argmin(r)
+mininds = np.argsort(r)
+center_parts = xp[mininds[:4]]
+xs = center_parts.sum(axis=0) / center_parts.shape[0]
+
+# Double-check all particles for boundaries
+for i in range(3):
+    mask = xp[:, i] < 0.0
+    xp[mask, i] += 1.0
+    mask = xp[:, i] > 1.0
+    xp[mask, i] -= 1.0
+
+unitL = cosmo_units["length"]
+edgelen = (2 * 260 * unyt.Mpc).to(unitL)
+boxsize = unyt.unyt_array([edgelen.v, edgelen.v, edgelen.v], unitL)
+
+xs = unyt.unyt_array(
+    [np.array([xs[0] * edgelen, xs[1] * edgelen, xs[2] * edgelen])], unitL
+)
+xp *= edgelen
+h *= edgelen
+
+w = Writer(unit_system=cosmo_units, box_size=boxsize, dimension=3)
+
+w.gas.coordinates = xp
+w.stars.coordinates = xs
+w.gas.velocities = np.zeros(xp.shape) * (unyt.cm / unyt.s)
+w.stars.velocities = np.zeros(xs.shape) * (unyt.cm / unyt.s)
+w.gas.masses = np.ones(xp.shape[0], dtype=float) * 1e1 * unyt.Msun
+w.stars.masses = np.ones(xs.shape[0], dtype=float) * 100.0 * unyt.Msun
+w.gas.internal_energy = (
+    np.ones(xp.shape[0], dtype=float) * (300.0 * unyt.kb * unyt.K) / unyt.g
+)
+
+w.gas.smoothing_length = h
+w.stars.smoothing_length = w.gas.smoothing_length[:1]
+
+w.write("propagationTest-3D.hdf5")
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/plotPhotonPropagationCheck.py b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/plotPhotonPropagationCheck.py
new file mode 100755
index 0000000000000000000000000000000000000000..04c154314899b24a36f451775f413caec70c72e9
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/plotPhotonPropagationCheck.py
@@ -0,0 +1,598 @@
+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Mladen Ivkovic (mladen.ivkovic@hotmail.com)
+#               2022 Tsang Keung Chan (chantsangkeung@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+# ----------------------------------------------------------------------
+# plots
+#   - radiation energies of particles as function of radius
+#   - magnitude of radiation fluxes of particles as function of radius
+#   - total energy in radial bins
+#   - total vectorial sum of fluxes in radial bins
+# and compare with expected propagation speed solution.
+# Usage:
+#   give snapshot number as cmdline arg to plot
+#   single snapshot, otherwise this script plots
+#   all snapshots available in the workdir.
+#   Make sure to select the photon group to plot that
+#   doesn't interact with gas to check the *propagation*
+#   correctly.
+# ----------------------------------------------------------------------
+
+import gc
+import os
+import sys
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+import unyt
+from matplotlib import pyplot as plt
+from scipy import stats
+from scipy.optimize import curve_fit
+
+import stromgren_plotting_tools as spt
+
+# Parameters users should/may tweak
+
+# snapshot basename
+snapshot_base = "propagation_test"
+
+time_first = 0
+
+# additional anisotropy estimate plot?
+plot_anisotropy_estimate = False
+
+# which photon group to use.
+# NOTE: array index, not group number (which starts at 1 for GEAR)
+group_index = 0
+
+scatterplot_kwargs = {
+    "alpha": 0.1,
+    "s": 1,
+    "marker": ".",
+    "linewidth": 0.0,
+    "facecolor": "blue",
+}
+
+lineplot_kwargs = {"linewidth": 2}
+
+# -----------------------------------------------------------------------
+
+# Read in cmdline arg: Are we plotting only one snapshot, or all?
+plot_all = False
+try:
+    snapnr = int(sys.argv[1])
+except IndexError:
+    plot_all = True
+
+mpl.rcParams["text.usetex"] = True
+
+
+def analytical_integrated_energy_solution(L, time, r, rmax):
+    """
+    Compute analytical solution for the sum of the energy
+    in bins for given injection rate <L> at time <time> 
+    at bin edges <r> and maximal radius <rmax>
+    """
+
+    r_center = 0.5 * (r[:-1] + r[1:])
+    r0 = r[0]
+    Etot = L * time
+
+    if rmax == 0:
+        return r_center, np.zeros(r.shape[0] - 1) * Etot.units
+
+    E = np.zeros(r.shape[0] - 1) * Etot.units
+    mask = r_center <= rmax
+    E[mask] = Etot / (rmax - r0) * (r[1:] - r[:-1])[mask]
+
+    return r_center, E
+
+
+def analytical_energy_solution(L, time, r, rmax):
+    """
+    Compute analytical solution for the energy distribution
+    for given injection rate <L> at time <time> at radii <r>
+    """
+
+    r_center = 0.5 * (r[:-1] + r[1:])
+    r0 = r[0]
+    Etot = L * time
+
+    if rmax == 0:
+        return r_center, np.zeros(r.shape[0] - 1) * Etot.units
+
+    E_fraction_bin = np.zeros(r.shape[0] - 1) * Etot.units
+    mask = r_center <= rmax
+    dr = r[1:] ** 2 - r[:-1] ** 2
+    E_fraction_bin[mask] = 1.0 / (rmax ** 2 - r0 ** 2) * dr[mask]
+    bin_surface = dr
+    total_weight = Etot / np.sum(E_fraction_bin / bin_surface)
+    E = E_fraction_bin / bin_surface * total_weight
+
+    return r_center, E
+
+
+def analytical_flux_magnitude_solution(L, time, r, rmax, scheme):
+    """
+    For radiation that doesn't interact with the gas, the
+    flux should correspond to the free streaming (optically
+    thin) limit. So compute and return that.
+    """
+    r, E = analytical_energy_solution(L, time, r, rmax)
+    if scheme.startswith("GEAR M1closure"):
+        F = unyt.c.to(r.units / time.units) * E / r.units ** 3
+    elif scheme.startswith("SPH M1closure"):
+        F = unyt.c.to(r.units / time.units) * E
+    else:
+        raise ValueError("Unknown scheme", scheme)
+
+    return r, F
+
+
+def x2(x, a, b):
+    return a * x ** 2 + b
+
+
+def get_snapshot_list(snapshot_basename="output"):
+    """
+    Find the snapshot(s) that are to be plotted 
+    and return their names as list
+    """
+
+    snaplist = []
+
+    if plot_all:
+        dirlist = os.listdir()
+        for f in dirlist:
+            if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+                snaplist.append(f)
+
+        snaplist = sorted(snaplist)
+
+    else:
+        fname = snapshot_basename + "_" + str(snapnr).zfill(4) + ".hdf5"
+        if not os.path.exists(fname):
+            print("Didn't find file", fname)
+            quit(1)
+        snaplist.append(fname)
+
+    return snaplist
+
+
+def plot_photons(filename, emin, emax, fmin, fmax):
+    """
+    Create the actual plot.
+
+    filename: file to work with
+    emin: list of minimal nonzero energy of all snapshots
+    emax: list of maximal energy of all snapshots
+    fmin: list of minimal flux magnitude of all snapshots
+    fmax: list of maximal flux magnitude of all snapshots
+    """
+    global time_first
+    print("working on", filename)
+
+    # Read in data first
+    data = swiftsimio.load(filename)
+    meta = data.metadata
+    scheme = str(meta.subgrid_scheme["RT Scheme"].decode("utf-8"))
+    boxsize = meta.boxsize
+    edgelen = min(boxsize[0], boxsize[1])
+
+    xstar = data.stars.coordinates
+    xpart = data.gas.coordinates
+    dxp = xpart - xstar
+    r = np.sqrt(np.sum(dxp ** 2, axis=1))
+
+    time = meta.time.copy()
+    # Take care of simulation time not starting at 0
+    if time_first == 0:
+        time_first = time
+        time = 0 * unyt.Gyr
+    else:
+        time -= time_first
+
+    r_expect = time * meta.reduced_lightspeed
+
+    L = None
+
+    use_const_emission_rates = False
+    if scheme.startswith("GEAR M1closure"):
+        luminosity_model = meta.parameters["GEARRT:stellar_luminosity_model"]
+        use_const_emission_rates = luminosity_model.decode("utf-8") == "const"
+    elif scheme.startswith("SPH M1closure"):
+        use_const_emission_rates = bool(
+            meta.parameters["SPHM1RT:use_const_emission_rates"]
+        )
+    else:
+        print("Error: Unknown RT scheme " + scheme)
+        exit()
+
+    if use_const_emission_rates:
+        # read emission rate parameter as string
+        if scheme.startswith("GEAR M1closure"):
+            const_emission_rates = (
+                spt.trim_paramstr(
+                    meta.parameters["GEARRT:const_stellar_luminosities_LSol"].decode(
+                        "utf-8"
+                    )
+                )
+                * unyt.L_Sun
+            )
+            L = const_emission_rates[group_index]
+        elif scheme.startswith("SPH M1closure"):
+            units = data.units
+            unit_l_in_cgs = units.length.in_cgs()
+            unit_v_in_cgs = (units.length / units.time).in_cgs()
+            unit_m_in_cgs = units.mass.in_cgs()
+            const_emission_rates = (
+                spt.trim_paramstr(
+                    meta.parameters["SPHM1RT:star_emission_rates"].decode("utf-8")
+                )
+                * unit_m_in_cgs
+                * unit_v_in_cgs ** 3
+                / unit_l_in_cgs
+            )
+            L = const_emission_rates[group_index]
+        else:
+            print("Error: Unknown RT scheme " + scheme)
+            exit()
+
+    if plot_anisotropy_estimate:
+        ncols = 4
+    else:
+        ncols = 3
+    fig = plt.figure(figsize=(5 * ncols, 5.5), dpi=200)
+
+    nbins = 100
+    r_bin_edges = np.linspace(0.5 * edgelen * 1e-3, 0.507 * edgelen, nbins + 1)
+    r_bin_centres = 0.5 * (r_bin_edges[1:] + r_bin_edges[:-1])
+    r_analytical_bin_edges = np.linspace(
+        0.5 * edgelen * 1e-6, 0.507 * edgelen, nbins + 1
+    )
+
+    # --------------------------
+    # Read in and process data
+    # --------------------------
+
+    energies = getattr(data.gas.photon_energies, "group" + str(group_index + 1))
+    Fx = getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "X")
+    Fy = getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "Y")
+    Fz = getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "Z")
+
+    fmag = np.sqrt(Fx ** 2 + Fy ** 2 + Fz ** 2)
+    particle_count, _ = np.histogram(
+        r,
+        bins=r_analytical_bin_edges,
+        range=(r_analytical_bin_edges[0], r_analytical_bin_edges[-1]),
+    )
+    L = L.to(energies.units / time.units)
+
+    xlabel_units_str = boxsize.units.latex_representation()
+    energy_units_str = energies.units.latex_representation()
+    flux_units_str = Fx.units.latex_representation()
+
+    # ------------------------
+    # Plot photon energies
+    # ------------------------
+    ax1 = fig.add_subplot(1, ncols, 1)
+    ax1.set_title("Particle Radiation Energies")
+    ax1.set_ylabel("Photon Energy [$" + energy_units_str + "$]")
+
+    # don't expect more than float precision
+    emin_to_use = max(emin, 1e-5 * emax)
+
+    if use_const_emission_rates:
+        # plot entire expected solution
+        rA, EA = analytical_energy_solution(L, time, r_analytical_bin_edges, r_expect)
+
+        mask = particle_count > 0
+        if mask.any():
+            EA = EA[mask].to(energies.units)
+            rA = rA[mask]
+            pcount = particle_count[mask]
+
+            # the particle bin counts will introduce noise.
+            # So use a linear fit for the plot. I assume here
+            # that the particle number per bin increases
+            # proprtional to r^2, which should roughly be the
+            # case for the underlying glass particle distribution.
+            popt, _ = curve_fit(x2, rA, pcount)
+            ax1.plot(
+                rA,
+                EA / x2(rA.v, popt[0], popt[1]),
+                **lineplot_kwargs,
+                linestyle="--",
+                c="red",
+                label="Analytical Solution",
+            )
+
+    else:
+        # just plot where photon front should be
+        ax1.plot(
+            [r_expect, r_expect],
+            [emin_to_use, emax * 1.1],
+            label="expected photon front",
+            color="red",
+        )
+
+    ax1.scatter(r, energies, **scatterplot_kwargs)
+    energies_binned, _, _ = stats.binned_statistic(
+        r,
+        energies,
+        statistic="mean",
+        bins=r_bin_edges,
+        range=(r_bin_edges[0], r_bin_edges[-1]),
+    )
+    ax1.plot(
+        r_bin_centres, energies_binned, **lineplot_kwargs, label="Mean Radiation Energy"
+    )
+    ax1.set_ylim(emin_to_use, emax * 1.1)
+
+    # ------------------------------
+    # Plot binned photon energies
+    # ------------------------------
+    ax2 = fig.add_subplot(1, ncols, 2)
+    ax2.set_title("Total Radiation Energy in radial bins")
+    ax2.set_ylabel("Total Photon Energy [$" + energy_units_str + "$]")
+
+    energies_summed_bin, _, _ = stats.binned_statistic(
+        r,
+        energies,
+        statistic="sum",
+        bins=r_bin_edges,
+        range=(r_bin_edges[0], r_bin_edges[-1]),
+    )
+    ax2.plot(
+        r_bin_centres,
+        energies_summed_bin,
+        **lineplot_kwargs,
+        label="Total Energy in Bin",
+    )
+    current_ylims = ax2.get_ylim()
+    ax2.set_ylim(emin_to_use, current_ylims[1])
+
+    if use_const_emission_rates:
+        # plot entire expected solution
+        # Note: you need to use the same bins as for the actual results
+        rA, EA = analytical_integrated_energy_solution(L, time, r_bin_edges, r_expect)
+
+        ax2.plot(
+            rA,
+            EA.to(energies.units),
+            **lineplot_kwargs,
+            linestyle="--",
+            c="red",
+            label="Analytical Solution",
+        )
+    else:
+        # just plot where photon front should be
+        ax2.plot(
+            [r_expect, r_expect],
+            ax2.get_ylim(r),
+            label="Expected Photon Front",
+            color="red",
+        )
+
+    # ------------------------------
+    # Plot photon fluxes
+    # ------------------------------
+    ax3 = fig.add_subplot(1, ncols, 3)
+    ax3.set_title("Particle Radiation Flux Magnitudes")
+    ax3.set_ylabel("Photon Flux Magnitude [$" + flux_units_str + "$]")
+
+    fmin_to_use = max(fmin, 1e-5 * fmax)
+    ax3.set_ylim(fmin_to_use, fmax * 1.1)
+
+    ax3.scatter(r, fmag, **scatterplot_kwargs)
+
+    fmag_mean_bin, _, _ = stats.binned_statistic(
+        r,
+        fmag,
+        statistic="mean",
+        bins=r_bin_edges,
+        range=(r_bin_edges[0], r_bin_edges[-1]),
+    )
+    ax3.plot(
+        r_bin_centres,
+        fmag_mean_bin,
+        **lineplot_kwargs,
+        label="Mean Radiation Flux of particles",
+    )
+
+    if use_const_emission_rates:
+        # plot entire expected solution
+        rA, FA = analytical_flux_magnitude_solution(
+            L, time, r_analytical_bin_edges, r_expect, scheme
+        )
+
+        mask = particle_count > 0
+        if mask.any():
+            FA = FA[mask].to(Fx.units)
+            rA = rA[mask]
+            pcount = particle_count[mask]
+
+            # the particle bin counts will introduce noise.
+            # So use a linear fit for the plot. I assume here
+            # that the particle number per bin increases
+            # proprtional to r^2, which should roughly be the
+            # case for the underlying glass particle distribution.
+            popt, _ = curve_fit(x2, rA, pcount)
+
+            ax3.plot(
+                rA,
+                FA / x2(rA.v, popt[0], popt[1]),
+                **lineplot_kwargs,
+                linestyle="--",
+                c="red",
+                label="analytical solution",
+            )
+
+    else:
+        # just plot where photon front should be
+        ax1.plot(
+            [r_expect, r_expect],
+            [emin_to_use, emax * 1.1],
+            label="expected photon front",
+            color="red",
+        )
+
+    # ------------------------------
+    # Plot photon flux sum
+    # ------------------------------
+
+    if plot_anisotropy_estimate:
+        ax4 = fig.add_subplot(1, ncols, 4)
+        ax4.set_title("Vectorial Sum of Radiation Flux in radial bins")
+        ax4.set_ylabel("[1]")
+
+        fmag_sum_bin, _, _ = stats.binned_statistic(
+            r,
+            fmag,
+            statistic="sum",
+            bins=r_bin_edges,
+            range=(r_bin_edges[0], r_bin_edges[-1]),
+        )
+        mask_sum = fmag_sum_bin > 0
+        fmag_max_bin, _, _ = stats.binned_statistic(
+            r,
+            fmag,
+            statistic="max",
+            bins=r_bin_edges,
+            range=(r_bin_edges[0], r_bin_edges[-1]),
+        )
+        mask_max = fmag_max_bin > 0
+        Fx_sum_bin, _, _ = stats.binned_statistic(
+            r,
+            Fx,
+            statistic="sum",
+            bins=r_bin_edges,
+            range=(r_bin_edges[0], r_bin_edges[-1]),
+        )
+        Fy_sum_bin, _, _ = stats.binned_statistic(
+            r,
+            Fy,
+            statistic="sum",
+            bins=r_bin_edges,
+            range=(r_bin_edges[0], r_bin_edges[-1]),
+        )
+        F_sum_bin = np.sqrt(Fx_sum_bin ** 2 + Fy_sum_bin ** 2)
+
+        ax4.plot(
+            r_bin_centres[mask_sum],
+            F_sum_bin[mask_sum] / fmag_sum_bin[mask_sum],
+            **lineplot_kwargs,
+            label=r"$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ "
+            + r"/ $\sum_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
+        )
+        ax4.plot(
+            r_bin_centres[mask_max],
+            F_sum_bin[mask_max] / fmag_max_bin[mask_max],
+            **lineplot_kwargs,
+            linestyle="--",
+            label=r"$\left| \sum_{i \in \mathrm{particles \ in \ bin}} \mathbf{F}_i \\right| $ "
+            + r" / $\max_{i \in \mathrm{particles \ in \ bin}} \left| \mathbf{F}_{i} \\right| $",
+        )
+
+    # -------------------------------------------
+    # Cosmetics that all axes have in common
+    # -------------------------------------------
+    for ax in fig.axes:
+        ax.set_xlabel("r [$" + xlabel_units_str + "$]")
+        ax.set_yscale("log")
+        ax.set_xlim(0.0, 0.501 * edgelen)
+        ax.legend(fontsize="x-small")
+
+    # Add title
+    title = filename.replace("_", r"\_")  # exception handle underscore for latex
+    if meta.cosmology is not None:
+        title += ", $z$ = {0:.2e}".format(meta.z)
+    title += ", $t$ = {0:.2e}".format(1 * meta.time)
+    fig.suptitle(title)
+
+    plt.tight_layout()
+    figname = filename[:-5]
+    figname += "-PhotonPropagation.png"
+    plt.savefig(figname)
+    plt.close()
+    gc.collect()
+
+    return
+
+
+def get_plot_boundaries(filenames):
+    """
+    Get minimal and maximal nonzero photon energy values
+    """
+
+    data = swiftsimio.load(filenames[0])
+    energies = (
+        1 * getattr(data.gas.photon_energies, "group" + str(group_index + 1))
+    ).value
+    emaxguess = energies.max()
+
+    emin = emaxguess
+    emax = 0.0
+    fmagmin = 1e30
+    fmagmax = -10.0
+
+    for f in filenames:
+        data = swiftsimio.load(f)
+
+        energies = (
+            1 * getattr(data.gas.photon_energies, "group" + str(group_index + 1))
+        ).value
+        mask = energies > 0.0
+
+        if mask.any():
+
+            nonzero_energies = energies[mask]
+            this_emin = nonzero_energies.min()
+            emin = min(this_emin, emin)
+
+            this_emax = energies.max()
+            emax = max(emax, this_emax)
+
+        fx = (
+            1 * getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "X")
+        ).value
+        fy = (
+            1 * getattr(data.gas.photon_fluxes, "Group" + str(group_index + 1) + "Y")
+        ).value
+        fmag = np.sqrt(fx ** 2 + fy ** 2)
+
+        fmagmin = min(fmagmin, fmag.min())
+        fmagmax = max(fmagmax, fmag.max())
+
+    return emin, emax, fmagmin, fmagmax
+
+
+if __name__ == "__main__":
+
+    print(
+        "REMINDER: Make sure you selected the correct photon group",
+        "to plot, which is hardcoded in this script.",
+    )
+    snaplist = get_snapshot_list(snapshot_base)
+    emin, emax, fmagmin, fmagmax = get_plot_boundaries(snaplist)
+
+    for f in snaplist:
+        plot_photons(f, emin, emax, fmagmin, fmagmax)
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/propagationTest-3D.yml b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/propagationTest-3D.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d64ea4c0a6f500084f5814281e8f2f7770505a1e
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/propagationTest-3D.yml
@@ -0,0 +1,67 @@
+MetaData:
+  run_name: StromgrenSpherePropagationTest3D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-12 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     2.e-02  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            propagation_test_low_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_low_redshift
+  scale_factor_first:  0.93
+  delta_time:          1.02
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.93
+  delta_time:          1.02 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./propagationTest-3D.hdf5     # The file to read
+  periodic:   1                             # peridioc ICs. Keep them periodic so we don't loose photon energy.
+
+GEARRT:
+  f_reduce_c: 1.                                    # reduce the speed of light for the RT solver by multiplying c with this factor
+  CFL_condition: 0.99
+  photon_groups_Hz: [0.]                            # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                   # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [1.]             # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                     # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                          # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17    # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1          # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                            # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                             # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                            # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                           # skip thermochemistry.
+  stars_max_timestep: 1.562500e-05                  # (Optional) restrict the maximal timestep of stars to this value (in internal units). Set to negative to turn off.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.93
+  a_end:          1.
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
+
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/runPropagationTest.sh b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/runPropagationTest.sh
new file mode 100755
index 0000000000000000000000000000000000000000..02756b7971e32084304af7359d1b22f34085bc66
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/runPropagationTest.sh
@@ -0,0 +1,32 @@
+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassCube_64.hdf5 ]
+then
+    echo "Fetching initial glass file for Strömgen Sphere 3D example ..."
+    ./getGlass.sh
+fi
+
+if [ ! -f 'propagationTest-3D.hdf5' ]; then
+    echo "Generating ICs"
+    python3 makePropagationTestIC.py
+fi
+
+# Run SWIFT with RT
+../../../swift \
+    --hydro  \
+    --cosmology \
+    --threads=4 \
+    --stars \
+    --external-gravity \
+    --feedback \
+    --radiation \
+    ./propagationTest-3D.yml 2>&1 | tee output.log
+
+# Plot the photon propagation checks.
+# Make sure you set the correct photon group to plot
+# inside the script
+python3 ./plotPhotonPropagationCheck.py
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/snapshot_times_low_redshift b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/snapshot_times_low_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e58febac09ba12e5677c6fc8abafe94cc31824ad
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/snapshot_times_low_redshift
@@ -0,0 +1,12 @@
+# Redshift
+0.0751932785953775
+0.06869649828543634
+0.06225945012746337
+0.055881161353223074
+0.04956068264358926
+0.04329708741463545
+0.03708947113706218
+0.030936950674366415
+0.02483866364616727
+0.018793767817153695
+0.01280144050017884
diff --git a/examples/RadiativeTransferTests/CosmoPropagationTest_3D/stromgren_plotting_tools.py b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/stromgren_plotting_tools.py
new file mode 120000
index 0000000000000000000000000000000000000000..2121a14d08c7cfcd516d98ec9338f0f46ee7358f
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoPropagationTest_3D/stromgren_plotting_tools.py
@@ -0,0 +1 @@
+../StromgrenSphere_3D/stromgren_plotting_tools.py
\ No newline at end of file
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/README b/examples/RadiativeTransferTests/CosmoUniformBox_3D/README
new file mode 100644
index 0000000000000000000000000000000000000000..a5dbcb88d1942dc08660d255d47e8c4813a1292e
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/README
@@ -0,0 +1,6 @@
+A simple box with static particles on a grid, and uniform radiation energy density everywhere. 
+Test that total radiation energy and radiation energy density dilute with the correct scale-factor relation.
+
+The ICs are created to be compatible with GEAR_RT. Recommended configuration:
+    --with-rt=GEAR_1 --with-rt-riemann-solver=GLF --with-hydro-dimension=3 --with-hydro=gizmo-mfv \
+     --with-riemann-solver=hllc --with-stars=GEAR --with-feedback=none --with-grackle=$GRACKLE_ROOT
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py b/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py
new file mode 100755
index 0000000000000000000000000000000000000000..647080dce6415abcfa37980eb434c296d241c37f
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/makeIC.py
@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+
+import h5py
+import numpy as np
+import unyt
+from swiftsimio import Writer
+from swiftsimio.units import cosmo_units
+
+# Unit system we're working with
+unitsystem = cosmo_units
+
+# Box is 100 Mpc in each direction
+boxsize = 260 * unyt.Mpc
+boxsize = boxsize.to(unitsystem["length"])
+
+reduced_speed_of_light_fraction = 1.0
+
+# Number of photon groups
+nPhotonGroups = 1
+
+# Number of particles in each dimension
+# Total number of particles is thus n_p^3
+n_p = 10
+
+# Filename of ICs to be generated
+outputfilename = "uniform_3D.hdf5"
+
+
+def initial_condition(unitsystem):
+    """
+    The initial conditions of the uniform box
+    
+    unitsystem: The unit system to use for IC
+
+    returns:
+    E: Photon energy density
+    F: Photon flux
+    """
+    # you can make the photon quantities unitless, the units will
+    # already have been written down in the writer.
+    # However, that means that you need to convert them manually.
+
+    unit_energy = (
+        unitsystem["mass"] * unitsystem["length"] ** 2 / unitsystem["time"] ** 2
+    )
+    unit_velocity = unitsystem["length"] / unitsystem["time"]
+    unit_flux = unit_energy * unit_velocity
+
+    c_internal = (unyt.c * reduced_speed_of_light_fraction).to(unit_velocity)
+
+    # Uniform energy
+    E = np.ones((n_p ** 3), dtype=np.float64) * unit_energy
+
+    # Assuming all photons flow in only one direction
+    # (optically thin regime, "free streaming limit"),
+    # we have that |F| = c * E
+    fluxes = np.zeros((3, n_p ** 3), dtype=np.float64)
+    fluxes[0] *= (E * c_internal / 1.73205).to(unit_flux)  # sqrt(3)
+    fluxes[1] *= (E * c_internal / 1.73205).to(unit_flux)  # sqrt(3)
+    fluxes[2] *= (E * c_internal / 1.73205).to(unit_flux)  # sqrt(3)
+
+    return E, fluxes.T
+
+
+if __name__ in ("__main__"):
+    # Coordinate array
+    coords = np.zeros((n_p ** 3, 3), dtype=np.float64)
+
+    # Calculate grid of evenly spaced coordinates
+    coords_per_dim = np.linspace(0.5, n_p - 0.5, n_p)
+    grid = np.meshgrid(coords_per_dim, coords_per_dim, coords_per_dim)
+
+    for i in range(3):
+        coords[:, i] = grid[i].flatten()
+
+    # Calculate and apply grid spacing
+    dx = boxsize / n_p
+    coords *= dx
+
+    w = Writer(unitsystem, boxsize, dimension=3)
+
+    w.gas.coordinates = coords
+    w.gas.velocities = np.zeros((n_p ** 3, 3)) * (unyt.cm / unyt.s)
+
+    mpart = 1e20 * unyt.M_sun
+    mpart = mpart.to(unitsystem["mass"])
+    w.gas.masses = np.ones(n_p ** 3, dtype=np.float64) * mpart
+    w.gas.internal_energy = (
+        np.ones(n_p ** 3, dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g
+    )
+
+    # Generate initial guess for smoothing lengths based on MIPS
+    w.gas.generate_smoothing_lengths(boxsize=boxsize, dimension=3)
+
+    # If IDs are not present, this automatically generates
+    w.write(outputfilename)
+
+    # Now open file back up again and add photon groups
+    # you can make them unitless, the units have already been
+    # written down in the writer, In this case, it's cosmo_units
+
+    F = h5py.File(outputfilename, "r+")
+    header = F["Header"]
+    nparts = header.attrs["NumPart_ThisFile"][0]
+    parts = F["/PartType0"]
+
+    # Generate initial conditions
+    E, fluxes = initial_condition(unitsystem)
+
+    # Create photon energy data entry
+    dsetname = "PhotonEnergiesGroup1"
+    energydata = np.zeros((nparts), dtype=np.float32)
+    parts.create_dataset(dsetname, data=E)
+
+    # Create photon fluxes data entry
+    dsetname = "PhotonFluxesGroup1"
+    fluxdata = np.zeros((nparts, 3), dtype=np.float32)
+    parts.create_dataset(dsetname, data=fluxes)
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py b/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py
new file mode 100755
index 0000000000000000000000000000000000000000..d55f837662bf5c5767f5ad1b1e2cbc8298fb61a7
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/plotSolution.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+import argparse
+
+import matplotlib as mpl
+import numpy as np
+import swiftsimio
+from matplotlib import pyplot as plt
+from scipy.optimize import curve_fit as cf
+
+# Parameters users should/may tweak
+snapshot_base = "output"  # snapshot basename
+plot_physical_quantities = True
+
+mpl.rcParams["text.usetex"] = True
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-z", "--redshift", help="Redshift domain to plot advection for", default="high"
+    )
+
+    args = parser.parse_args()
+    return args
+
+
+def get_snapshot_list(snapshot_basename="output"):
+    """
+    Find the snapshots that are to be plotted 
+    and return their names as list
+    """
+    snaplist = []
+
+    dirlist = os.listdir()
+    for f in dirlist:
+        if f.startswith(snapshot_basename) and f.endswith("hdf5"):
+            snaplist.append(f)
+
+    snaplist = sorted(snaplist)
+    if len(snaplist) == 0:
+        print(f"No snapshots with base {snapshot_basename} found!")
+        sys.exit(1)
+    return snaplist
+
+
+def plot_param_over_time(
+    snapshot_list, param="energy density", redshift_domain="high_redshift"
+):
+    print(f"Now plotting {param} over time")
+
+    # Arrays to keep track of plot_param and scale factor
+    plot_param = [[], []]
+    scale_factor = []
+    analytic_exponent = [0, 0]
+
+    # Functions to convert between scale factor and redshift
+    a2z = lambda a: 1 / a - 1
+    z2a = lambda z: 1 / (z + 1)
+
+    for file in snapshot_list:
+        data = swiftsimio.load(file)
+        meta = data.metadata
+
+        # Read comoving quantities
+        energy = getattr(data.gas.photon_energies, "group1")
+        mass = data.gas.masses
+        rho = data.gas.densities
+        vol = mass / rho
+
+        energy_density = energy / vol
+
+        if plot_physical_quantities:
+            physical_energy_density = energy_density.to_physical()
+            physical_mass = mass.to_physical()
+            physical_vol = vol.to_physical()
+            physical_energy = physical_energy_density * physical_vol
+            if param == "energy density":
+                plot_param[1].append(
+                    1
+                    * np.sum(physical_energy_density)
+                    / physical_energy_density.shape[0]
+                )
+                analytic_exponent[1] = -3.0
+            elif param == "volume":
+                plot_param[1].append(1 * np.sum(physical_vol) / physical_vol.shape[0])
+                analytic_exponent[1] = 3.0
+            elif param == "total energy":
+                plot_param[1].append(1 * np.sum(physical_energy))
+                analytic_exponent[1] = 0.0
+            elif param == "mass":
+                plot_param[1].append(1 * np.sum(physical_mass))
+                analytic_exponent[1] = 0.0
+
+        if param == "energy density":
+            plot_param[0].append(1 * np.sum(energy_density) / energy_density.shape[0])
+            analytic_exponent[0] = 0.0
+        elif param == "volume":
+            plot_param[0].append(1 * np.sum(vol) / vol.shape[0])
+            analytic_exponent[0] = 0.0
+        elif param == "total energy":
+            plot_param[0].append(1 * np.sum(energy))
+            analytic_exponent[0] = 0.0
+        elif param == "mass":
+            plot_param[0].append(1 * np.sum(mass))
+            analytic_exponent[0] = 0.0
+
+        scale_factor.append(meta.scale_factor)
+
+    fig = plt.figure(figsize=(5.05 * (1 + plot_physical_quantities), 5.4), dpi=200)
+
+    x = np.linspace(min(scale_factor), max(scale_factor), 1000)
+
+    if param == "energy density":
+        titles = ["Comoving energy density", "Physical energy density"]
+        ylabel = "Average energy density"
+        figname = f"output_energy_density_over_time-{redshift_domain}.png"
+    elif param == "volume":
+        titles = ["Comoving particle volume", "Physical particle volume"]
+        ylabel = "Average particle volume"
+        figname = f"output_volume_over_time-{redshift_domain}.png"
+    elif param == "total energy":
+        titles = ["Comoving total energy", "Physical total energy"]
+        ylabel = "Total energy"
+        figname = f"output_total_energy_over_time-{redshift_domain}.png"
+    elif param == "mass":
+        titles = ["Comoving total mass", "Physical total mass"]
+        ylabel = "Total mass"
+        figname = f"output_total_mass_over_time-{redshift_domain}.png"
+
+    for i in range(1 + plot_physical_quantities):
+        ax = fig.add_subplot(1, (1 + plot_physical_quantities), (1 + i))
+        ax.scatter(scale_factor, plot_param[i], label="Simulation")
+
+        # Analytic scale-factor relation
+        analytic = x ** analytic_exponent[i]
+
+        # Scale solution to correct offset
+        analytic = analytic / analytic[0] * plot_param[i][0]
+        ax.plot(
+            x,
+            analytic,
+            c="r",
+            label=f"Analytic solution $\propto a^{{{analytic_exponent[i]}}}$",
+        )
+
+        ax.legend()
+        ax.set_title(titles[i])
+
+        ax.set_xlabel("Scale factor")
+        secax = ax.secondary_xaxis("top", functions=(a2z, z2a))
+        secax.set_xlabel("Redshift")
+
+        ax.yaxis.get_offset_text().set_position((-0.05, 1))
+
+        if analytic_exponent[i] == 0.0:
+            ax.set_ylim(plot_param[i][0] * 0.95, plot_param[i][0] * 1.05)
+        if i == 0:
+            units = plot_param[i][0].units.latex_representation()
+            ax.set_ylabel(f"{ylabel} [${units}$]")
+
+    plt.tight_layout()
+    plt.savefig(figname)
+    plt.close()
+
+
+if __name__ in ("__main__"):
+    # Get command line args
+    args = parse_args()
+    domain = args.redshift.lower()
+    if domain in ("low", "l", "low_redshift", "low redshift", "low-redshift"):
+        redshift_domain = "low_redshift"
+    elif domain in (
+        "medium",
+        "m",
+        "medium_redshift",
+        "medium redshift",
+        "medium-redshift",
+    ):
+        redshift_domain = "medium_redshift"
+    elif domain in ("high", "h", "high_redshift", "high redshift", "high-redshift"):
+        redshift_domain = "high_redshift"
+    else:
+        print("Redshift domain not recognised!")
+        sys.exit(1)
+
+    snaplist = get_snapshot_list(snapshot_base + f"_{redshift_domain}")
+
+    if len(snaplist) < 1:
+        print("No snapshots found!")
+        exit(1)
+
+    for param in ["energy density", "volume", "total energy", "mass"]:
+        plot_param_over_time(snaplist, param, redshift_domain)
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..f4177d6804a5865a530e5085f8f7032144ea4b98
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_high_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: cosmo_RT_uniform_box-3D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2   # The maximal time-step size of the simulation (in internal units).
+  time_begin: 0
+  time_end:   1.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_high_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_high_redshift
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.00990099
+  delta_time:          1.06   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./uniform_3D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0                       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99                             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [1.]                          # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                 # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [0.]           # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                   # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                        # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17  # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1        # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                         # Skip thermochemsitry.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.00990099  # z~100
+  a_end:          0.01408451  # z~70
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6bac120ae6d9d4e8dcf0239a1e656073d52b97fa
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_low_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: cosmo_RT_uniform_box-3D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2   # The maximal time-step size of the simulation (in internal units).
+  time_begin: 0
+  time_end:   1.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_low_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_low_redshift
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.93
+  delta_time:          1.005   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./uniform_3D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0                       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99                             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [1.]                          # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                 # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [0.]           # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                   # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                        # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17  # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1        # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                         # Skip thermochemsitry.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.93 
+  a_end:          1.     
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9189a2ca73dbfa70d48a3e9b1765b82c8ab8e864
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/rt_uniform3D_medium_redshift.yml
@@ -0,0 +1,69 @@
+MetaData:
+  run_name: cosmo_RT_uniform_box-3D
+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841586e+33 # 1 M_Sol
+  UnitLength_in_cgs:   3.08567758e21 # Mpc in cm
+  UnitVelocity_in_cgs: 1.e5 # km/s
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1. # K
+
+# Parameters governing the time integration
+TimeIntegration:
+  max_nr_rt_subcycles: 1
+  dt_min:     1.e-17  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1.e-2   # The maximal time-step size of the simulation (in internal units).
+  time_begin: 0
+  time_end:   1.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_medium_redshift # Common part of the name of output files
+  output_list_on:      1
+  output_list:         snapshot_times_medium_redshift
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.0909
+  delta_time:          1.02   # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.6      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./uniform_3D.hdf5  # The file to read
+  periodic:   1                     # periodic ICs
+
+Restarts:
+  delta_hours:        72        # (Optional) decimal hours between dumps of restart files.
+
+GEARRT:
+  f_reduce_c: 1.                                  # reduce the speed of light for the RT solver by multiplying c with this factor
+  f_limit_cooling_time: 0.0                       # (Optional) multiply the cooling time by this factor when estimating maximal next time step. Set to 0.0 to turn computation of cooling time off.
+  CFL_condition: 0.99                             # CFL condition for RT, independent of hydro
+  photon_groups_Hz: [1.]                          # Lower photon frequency group bin edges in Hz. Needs to have exactly N elements, where N is the configured number of bins --with-RT=GEAR_N
+  stellar_luminosity_model: const                 # Which model to use to determine the stellar luminosities.
+  const_stellar_luminosities_LSol: [0.]           # (Conditional) constant star luminosities for each photon frequency group to use if stellar_luminosity_model:const is set, in units of Solar Luminosity.
+  hydrogen_mass_fraction:  0.76                   # total hydrogen (H + H+) mass fraction in the metal-free portion of the gas
+  stellar_spectrum_type: 0                        # Which radiation spectrum to use. 0: constant from 0 until some max frequency set by stellar_spectrum_const_max_frequency_Hz. 1: blackbody spectrum.
+  stellar_spectrum_const_max_frequency_Hz: 1.e17  # (Conditional) if stellar_spectrum_type=0, use this maximal frequency for the constant spectrum. 
+  set_initial_ionization_mass_fractions: 1        # (Optional) manually overwrite initial mass fraction of each species (using the values you set below)
+  mass_fraction_HI: 0.76                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HI mass fractions to this value
+  mass_fraction_HII: 0.                           # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HII mass fractions to this value
+  mass_fraction_HeI: 0.24                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeI mass fractions to this value
+  mass_fraction_HeII: 0.                          # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeII mass fractions to this value
+  mass_fraction_HeIII: 0.                         # (Conditional) If overwrite_initial_ionization_fractions=1, needed to set initial HeIII mass fractions to this value
+  skip_thermochemistry: 1                         # Skip thermochemsitry.
+
+Cosmology:        # Planck13 (EAGLE flavour)
+  a_begin:        0.0909    # z=10
+  a_end:          0.1428571 # z=6
+  h:              0.6777
+  Omega_cdm:      0.2587481
+  Omega_lambda:   0.693
+  Omega_b:        0.0482519
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/run.sh b/examples/RadiativeTransferTests/CosmoUniformBox_3D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..42c200af6c300fa2913a907cab764b22c0fadb4a
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/run.sh
@@ -0,0 +1,52 @@
+#!/bin/bash
+
+# exit if anything fails
+set -e
+set -o pipefail
+
+# Generate the initial conditions if they are not present.
+if [ ! -e uniform_3D.hdf5 ]
+then
+    echo "Generating initial conditions for the 3D uniform box example..."
+    python3 makeIC.py
+fi
+
+# Default run
+ymlfile=rt_uniform3D_high_redshift.yml
+zdomain="h"
+
+# Do we have a cmdline argument provided?
+if [ $# -gt 0 ]; then
+    case "$1" in
+    -l | -low | --l | --low | l | ./rt_uniform3D_low_redshift.yml | rt_uniform3D_low_redshift | rt_uniform3D_low_redshift.yml )
+        ymlfile=rt_uniform3D_low_redshift.yml
+	zdomain="l"
+        ;;
+    -m | -mid | --m | --mid | m | ./rt_uniform3D_medium_redshift.yml | rt_uniform3D_medium_redshift | rt_uniform3D_medium_redshift.yml )
+        ymlfile=rt_uniform3D_medium_redshift.yml
+	zdomain="m"
+        ;;
+    -h | -high | -hi | --h | --hi | --high | h | ./rt_uniform3D_high_redshift.yml | rt_uniform3D_high_redshift | rt_uniform3D_high_redshift.yml )
+        ymlfile=rt_uniform3D_high_redshift.yml
+	zdomain="h"
+        ;;
+    *)
+        echo unknown cmdline param, running default $ymlfile
+        ;;
+    esac
+fi
+
+
+# Run SWIFT with RT and cosmology
+../../../swift \
+    --hydro \
+    --cosmology \
+    --threads=4 \
+    --verbose=0  \
+    --radiation \
+    --stars \
+    --feedback \
+    --external-gravity \
+    $ymlfile 2>&1 | tee output.log
+
+python3 ./plotSolution.py -z $zdomain
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_high_redshift b/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_high_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e504014177837a7b2c1d134d1c4a0d5a785d0044
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_high_redshift
@@ -0,0 +1,12 @@
+# Redshift
+99.98460572382712
+96.47031888378449
+93.23367859034735
+90.25299848922855
+87.49735694548497
+84.9407463291136
+82.56112036003019
+80.33965676782498
+78.26018022924342
+76.30870609169932
+74.47307621234665
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_low_redshift b/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_low_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..e58febac09ba12e5677c6fc8abafe94cc31824ad
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_low_redshift
@@ -0,0 +1,12 @@
+# Redshift
+0.0751932785953775
+0.06869649828543634
+0.06225945012746337
+0.055881161353223074
+0.04956068264358926
+0.04329708741463545
+0.03708947113706218
+0.030936950674366415
+0.02483866364616727
+0.018793767817153695
+0.01280144050017884
diff --git a/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_medium_redshift b/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_medium_redshift
new file mode 100644
index 0000000000000000000000000000000000000000..958165894af540df676443cfdf3fbf20ce46dd1c
--- /dev/null
+++ b/examples/RadiativeTransferTests/CosmoUniformBox_3D/snapshot_times_medium_redshift
@@ -0,0 +1,12 @@
+# Redshift
+9.983090485639307
+9.499283920730043
+8.92297287774699
+8.419654113369706
+7.975570397565402
+7.580294600085116
+7.225767198868688
+6.905651913777566
+6.614892571630509
+6.349401127100822
+6.1058334302500645
diff --git a/src/cosmology.c b/src/cosmology.c
index 23915f925378b2f5dd76e7b23cd858ef85c40d44..269c808e07d46b38e5b7364823973e7b4aaf5f79 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -813,7 +813,9 @@ void cosmology_init_tables(struct cosmology *c) {
 
     /* Choose value of expansion factor for check */
     const double dloga = (c->log_a_end - c->log_a_begin) / (n - 1);
-    const double a = exp(c->log_a_begin + dloga * i);
+    double a = exp(c->log_a_begin + dloga * i);
+    a = fmax(a, c->a_begin);
+    a = fmin(a, c->a_end);
 
     /* Verify that converting expansion factor to time and back recovers the
      * original value */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 171484205512e86718f4c5de60d826b0c77c7985..3d2e3cbded1866cf6c94c53aa86d8e8b46a76c2a 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -207,7 +207,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
       "co-moving positions of the particles");
 
   list[2] =
-      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 1.f, parts,
+      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts,
                            conserved.mass, "Co-moving masses of the particles");
 
   list[3] = io_make_output_field(
diff --git a/src/rt/GEAR/rt.h b/src/rt/GEAR/rt.h
index 72a9e25b221c162063dbffbad0e8d84abe6ed37e..f2c2f31e7f6f709c591fad8fdb62494b7776e6dc 100644
--- a/src/rt/GEAR/rt.h
+++ b/src/rt/GEAR/rt.h
@@ -402,8 +402,7 @@ __attribute__((always_inline)) INLINE static double rt_part_dt(
     const double time_base, const int with_cosmology,
     const struct cosmology* cosmo) {
   if (with_cosmology) {
-    error("GEAR RT with cosmology not implemented yet! :(");
-    return 0.f;
+    return cosmology_get_delta_time(cosmo, ti_beg, ti_end);
   } else {
     return (ti_end - ti_beg) * time_base;
   }
@@ -483,12 +482,17 @@ __attribute__((always_inline)) INLINE static void rt_finalise_transport(
 
   for (int g = 0; g < RT_NGROUPS; g++) {
     const float e_old = rtd->radiation[g].energy_density;
+
     /* Note: in this scheme, we're updating d/dt (U * V) + sum F * A * dt = 0.
      * So we'll need the division by the volume here. */
     rtd->radiation[g].energy_density += rtd->flux[g].energy * Vinv;
+
     rtd->radiation[g].flux[0] += rtd->flux[g].flux[0] * Vinv;
+
     rtd->radiation[g].flux[1] += rtd->flux[g].flux[1] * Vinv;
+
     rtd->radiation[g].flux[2] += rtd->flux[g].flux[2] * Vinv;
+
     rt_check_unphysical_state(&rtd->radiation[g].energy_density,
                               rtd->radiation[g].flux, e_old, /*callloc=*/4);
   }
diff --git a/src/timestep.h b/src/timestep.h
index c09895bac00e2793b621d68697d7d51f33839b7d..5e5af8074632fcfaa486e26c23f2fc1e170c3748 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -259,12 +259,11 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_rt_timestep(
                           e->physical_constants, e->internal_units);
 
   if ((e->policy & engine_policy_cosmology))
-    error("Cosmology factor in get_part_rt_timestep not implemented yet");
-  /* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/
-  /* new_dt = min(new_dt, e->dt_max_RMS_displacement); */
+    /* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/
+    new_dt = min(new_dt, e->dt_max_RMS_displacement);
 
   /* Apply cosmology correction (This is 1 if non-cosmological) */
-  /* new_dt *= e->cosmology->time_step_factor; */
+  new_dt *= e->cosmology->time_step_factor;
 
   /* Limit timestep within the allowed range */
   new_dt = min(new_dt, e->dt_max);
diff --git a/swift.c b/swift.c
index 038f078a68031c585e17f93df9cbda7f318db60f..34760a34e7606fa8077bc65cef59baae0e34a6e9 100644
--- a/swift.c
+++ b/swift.c
@@ -670,9 +670,6 @@ int main(int argc, char *argv[]) {
   if (with_rt && with_cooling) {
     error("Error: Cannot use radiative transfer and cooling simultaneously.");
   }
-  if (with_rt && with_cosmology) {
-    error("Error: Cannot use run radiative transfer with cosmology (yet).");
-  }
 #endif /* idfef RT_NONE */
 
 #ifdef SINK_NONE