diff --git a/examples/MHDTests/MagnetisedCloudCollapse/makeIC.py b/examples/MHDTests/MagnetisedCloudCollapse/makeIC.py
index 870c2fa64ce57925f38a777aec9a557ab2839f00..3b681dcfe51328f8b8374be596ed9dc5df7690f6 100644
--- a/examples/MHDTests/MagnetisedCloudCollapse/makeIC.py
+++ b/examples/MHDTests/MagnetisedCloudCollapse/makeIC.py
@@ -40,7 +40,7 @@ Omega = 2 * pi / (T * 3.1536e7) # initial angular frequency of cloud
mu = (
10
) # mass to magnetic field flux through sphere, normalised to a critical value for collapse. Refer to e.g. Henebelle & Fromang 2008 for details.
-Bini = 3.0 / c1 * sqrt(mu0 * G / 5.0) * M / (Rcloud * Rcloud) * 1 / mu
+Bini = (3.0 / (2.0 * c1)) * sqrt(mu0 * G / (5.0 * pi)) * M / (Rcloud * Rcloud) * 1 / mu
# Barotropic EoS parameters
cs0 = 2e4
diff --git a/examples/MHDTests/MagnetisedCloudCollapse/makePrimitiveCubicLattice.py b/examples/MHDTests/MagnetisedCloudCollapse/makePrimitiveCubicLattice.py
new file mode 100644
index 0000000000000000000000000000000000000000..e8685092f18e03e16731e9780b08894919cd2281
--- /dev/null
+++ b/examples/MHDTests/MagnetisedCloudCollapse/makePrimitiveCubicLattice.py
@@ -0,0 +1,72 @@
+"""
+Create a cubic lattice.
+"""
+
+import numpy as np
+import h5py
+
+from unyt import cm, g, s, erg
+from unyt.unit_systems import cgs_unit_system
+
+from swiftsimio import Writer
+
+
+def generate_cube(num_on_side, side_length=1.0):
+ """
+ Generates a cube
+ """
+
+ values, step = np.linspace(0.0, side_length, num_on_side, endpoint=False, retstep=True)
+ values += 0.5*step
+
+ positions = np.empty((num_on_side ** 3, 3), dtype=float)
+
+ for x in range(num_on_side):
+ for y in range(num_on_side):
+ for z in range(num_on_side):
+ index = x * num_on_side + y * num_on_side ** 2 + z
+
+ positions[index, 0] = values[x]
+ positions[index, 1] = values[y]
+ positions[index, 2] = values[z]
+
+ return positions
+
+
+def write_out_glass(filename, cube, side_length=1.0):
+ x = Writer(cgs_unit_system, side_length * cm)
+
+ x.gas.coordinates = cube * cm
+
+ x.gas.velocities = np.zeros_like(cube) * cm / s
+
+ x.gas.masses = np.ones(cube.shape[0], dtype=float) * g
+
+ x.gas.internal_energy = np.ones(cube.shape[0], dtype=float) * erg / g
+
+ x.gas.generate_smoothing_lengths(boxsize=side_length * cm, dimension=3)
+
+ x.write(filename)
+
+ return
+
+
+if __name__ == "__main__":
+ import argparse as ap
+
+ parser = ap.ArgumentParser(description="Generate a BCC lattice")
+
+ parser.add_argument(
+ "-n",
+ "--numparts",
+ help="Number of particles on a side. Default: 32",
+ default=32,
+ type=int,
+ )
+
+ args = parser.parse_args()
+
+ output = "CubicPrimitive_{}.hdf5".format(args.numparts)
+
+ glass_cube = generate_cube(args.numparts)
+ write_out_glass(output, glass_cube)
diff --git a/examples/MHDTests/MagnetisedCloudCollapse/plotSolutionReference.py b/examples/MHDTests/MagnetisedCloudCollapse/plotSolutionReference.py
new file mode 100644
index 0000000000000000000000000000000000000000..09c95044e8c46daf57a09035f79aedc290d32c3f
--- /dev/null
+++ b/examples/MHDTests/MagnetisedCloudCollapse/plotSolutionReference.py
@@ -0,0 +1,204 @@
+# colormaps for the script were chosen to match 1909.09650v2
+
+import numpy as np
+import h5py
+import sys
+import unyt
+import matplotlib.pyplot as plt
+
+from swiftsimio import load
+from swiftsimio.visualisation.rotation import rotation_matrix_from_vector
+from swiftsimio.visualisation.slice import slice_gas
+
+filename = sys.argv[1]
+data = load(filename)
+center = 0.5 * data.metadata.boxsize
+
+# Constants
+G = 6.67430e-11 * unyt.m **3 / unyt.kg / unyt.s **2 #6.67430e-8
+G = G.to(unyt.cm **3 / unyt.g / unyt.s **2)
+mu0 = 1.25663706127e-1 * unyt.g * unyt.cm / (unyt.s**2 * unyt.A**2)
+
+# Parameters (taken from Hopkins 2016)
+R0 = 4.628516371e16 * unyt.cm
+M = 1.99e33 * unyt.g # total mass of the sphere
+rhocloud0 = M/(4/3 * np.pi * R0**3 )
+rhouniform = rhocloud0/360
+t_ff = np.sqrt(3/(2*np.pi*G*rhocloud0)) * unyt.s
+tsim = data.metadata.time
+toplot = tsim / t_ff
+print("Showing results at %3f free fall times" % toplot)
+
+# First create a mass-weighted temperature dataset
+r = data.gas.coordinates - center
+v = data.gas.velocities
+norm_r = np.sqrt(r[:, 0] ** 2 + r[:, 1] ** 2 + r[:, 2] ** 2)
+vr = np.sum(v * r, axis=1) / norm_r
+B = data.gas.magnetic_flux_densities
+J = data.gas.magnetic_flux_curl
+normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
+normJ = np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2)
+divB = data.gas.magnetic_divergences
+h = data.gas.smoothing_lengths
+rho = data.gas.densities
+
+
+data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities
+data.gas.mass_weighted_error = data.gas.masses * np.log10(
+ np.maximum(h * abs(divB) / normB, 1e-6)
+)
+data.gas.mass_weighted_vr = data.gas.masses * vr
+data.gas.mass_weighted_Bz = data.gas.masses * abs(B[:, 2])
+data.gas.mass_weighted_J = data.gas.masses * np.sqrt(J[:, 0] ** 2 + J[:, 1] ** 2 + J[:, 2] ** 2) / mu0
+
+"""
+rotation_matrix = [[1,0,0],
+ [0,1,0],
+ [0,0,1]]
+"""
+rotation_matrix = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]
+
+Lslice_AU = 1000
+Lslice = Lslice_AU * 1.49597871 * 10**13 * unyt.cm
+
+visualise_region_xy = [
+ center[0] - Lslice,
+ center[0] + Lslice,
+ center[1] - Lslice,
+ center[1] + Lslice,
+]
+
+visualise_region = [
+ center[0] - Lslice,
+ center[0] + Lslice,
+ center[2] - Lslice,
+ center[2] + Lslice,
+]
+
+
+common_arguments_xy = dict(
+ data=data,
+ resolution=512,
+ parallel=True,
+ region=visualise_region_xy,
+ z_slice=center[2], #0.0 * unyt.cm,
+)
+
+common_arguments = dict(
+ data=data,
+ resolution=512,
+ parallel=True,
+ rotation_center=unyt.unyt_array(center),
+ rotation_matrix=rotation_matrix,
+ region=visualise_region,
+ z_slice=0.0 * unyt.cm,
+)
+
+
+# Map in msun / mpc^3
+
+mass_map_xy = slice_gas(**common_arguments_xy, project="masses")
+mass_weighted_density_map_xy = slice_gas(
+ **common_arguments_xy, project="mass_weighted_densities"
+)
+
+mass_map = slice_gas(**common_arguments, project="masses")
+mass_weighted_density_map = slice_gas(
+ **common_arguments, project="mass_weighted_densities"
+)
+mass_weighted_error_map = slice_gas(**common_arguments, project="mass_weighted_error")
+mass_weighted_vr_map = slice_gas(**common_arguments, project="mass_weighted_vr")
+mass_weighted_Bz_map = slice_gas(**common_arguments, project="mass_weighted_Bz")
+mass_weighted_J_map = slice_gas(**common_arguments, project="mass_weighted_J")
+
+density_map_xy = mass_weighted_density_map_xy / mass_map_xy
+density_map = mass_weighted_density_map / mass_map
+error_map = mass_weighted_error_map / mass_map
+vr_map = mass_weighted_vr_map / mass_map
+Bz_map = mass_weighted_Bz_map / mass_map
+J_map = mass_weighted_J_map / mass_map
+
+density_map_xy.convert_to_units(unyt.g * unyt.cm ** (-3))
+density_map.convert_to_units(unyt.g * unyt.cm ** (-3))
+vr_map.convert_to_units(unyt.km / unyt.s)
+Bz_map.convert_to_units(1e-7 * unyt.g / (unyt.A * unyt.s * unyt.s))
+J_map.convert_to_units(unyt.A / (unyt.m**2))
+
+fig, ax = plt.subplots(3, 2, sharey=True, figsize=(10, 15))
+
+fig.suptitle("Plotting at %.3f free fall times" % toplot)
+
+a00 = ax[0, 0].contourf(
+ np.log10(density_map.value),
+ cmap="RdBu",
+ extend="both",
+ levels=np.linspace(-17.0, -11.0, 100),
+)
+
+a01 = ax[0, 1].contourf(
+ np.log10(density_map_xy.value),
+ cmap="RdBu",
+ extend="both",
+ levels=np.linspace(-17.0, -11.0, 100),
+)
+
+a10 = ax[1, 0].contourf(
+ vr_map.value, cmap="CMRmap", extend="both", levels=np.linspace(-1.0, 1.0, 100)
+)
+a11 = ax[1, 1].contourf(
+ error_map.value, cmap="jet", extend="both", levels=np.arange(-6.0, 3.0,1.0)
+)
+a20 = ax[2, 0].contourf(
+ np.log10(np.maximum(J_map.value, -15)),
+ cmap='nipy_spectral',
+ extend="both",
+ levels=np.linspace(-15, -11, 100),
+)
+a21 = ax[2, 1].contourf(
+ np.log10(np.maximum(Bz_map.value, -6.9)),
+ cmap="nipy_spectral_r",
+ extend="both",
+ levels=np.linspace(2.0, 5.0, 100),
+)
+
+locs = [512/4,512/2,3*512/4]
+labels = [-Lslice_AU/2,0,Lslice_AU/2]
+
+for ii in range(3):
+ ax[ii, 0].set_ylabel(r"$z$ [A.U.]")
+ ax[ii, 0].set_yticks(locs, labels)
+
+for ii in range(3):
+ for jj in range(2):
+ ax[ii, jj].set_xlabel(r"$x$ [A.U.]")
+ ax[ii, jj].set_xticks(locs, labels)
+ ax[ii, jj].set_aspect("equal")
+
+cbar1 = plt.colorbar(
+ a00, ax=ax[0, 0], fraction=0.046, pad=0.04, ticks=np.linspace(-17, -11, 7)
+)
+cbar1.set_label(r"$\mathrm{log}_{10} \rho \; [g {cm}^{-3}]$")
+cbar2 = plt.colorbar(
+ a01, ax=ax[0, 1], fraction=0.046, pad=0.04, ticks=np.linspace(-17, -11, 7)
+)
+ax[0, 1].set_ylabel(r"$y$ [A.U.]")
+cbar2.set_label(r"$\mathrm{log}_{10} \rho \; [g {cm}^{-3}]$")
+cbar3 = plt.colorbar(
+ a10, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=np.linspace(-1.0, 1.0, 9)
+)
+cbar3.set_label(r"$\mathrm{log}_{10} v_r \; [km s^{-1}]$")
+
+cbar4 = plt.colorbar(a11, ax=ax[1, 1], fraction=0.046, pad=0.04)
+cbar4.set_label(r"$\mathrm{log}_{10} \frac{h \: \nabla \cdot B}{|B|}$")
+cbar5 = plt.colorbar(
+ a20, ax=ax[2, 0], fraction=0.046, pad=0.04, ticks=np.linspace(-15, -11, 5)
+)
+cbar5.set_label(r"$\mathrm{log}_{10} |J| \; [A/m^2]$")
+cbar6 = plt.colorbar(
+ a21, ax=ax[2, 1], fraction=0.046, pad=0.04, ticks=np.linspace(2, 5, 4)
+)
+cbar6.set_label(r"$\mathrm{log}_{10} B_z \; [\mu G]$")
+
+plt.tight_layout()
+
+plt.savefig(sys.argv[2])