diff --git a/examples/GravityTests/Plummer_Selfgravitating/README b/examples/GravityTests/Plummer_Selfgravitating/README
new file mode 100644
index 0000000000000000000000000000000000000000..b6a061b5d45c1b20778613d7c63c9aa5ca3e7cbe
--- /dev/null
+++ b/examples/GravityTests/Plummer_Selfgravitating/README
@@ -0,0 +1,50 @@
+This example generates a self consistent realization of an anisotropic Plummer model, by sampling
+the distribution function derived by (Dejonghe 1987).
+
+The parameters that control the ICs are essentially:
+*  a: 			The Plummer softening length
+*  M:			The total mass of the system
+*  q:			A parameter that controls the radial anisotropy of the system,
+			ranging from -inf to 2, where -inf corresponds to the "Einstein Sphere",
+			a model with only circular orbits and 2 to a model with only radial orbits.
+			q is related to the traditional beta parameter (which measures anisotropy) via
+			beta(r) = q/2 * r^2/(1+r^2)
+*  N:			The Number of particles in the system
+
+Some additional parameters are:
+*  bound:		A given cut radius above which particles are not included in the IC file
+			(this may reduce the effective number of particles to N_eff < N). N_eff is
+			printed out when the script is executed.
+*  use_parallel:	Use python's multiprocessing library for sampling
+*  nb_threads:		Number of threads to use for the above
+*  pickle_ics:		Also output .pkl files of the pos, vel and mass arrays
+
+All other variables should be self-explanatory.
+
+The radii of the particles are renerated through inverse transform sampling, exploiting
+the analiticity of the inverse of M(<r) for the Plummer model. (Masses are drawn uniformly
+between 0 and M, and their corresponding radii are computed via r(M))
+
+The radial and tangential components of the velocities are generated trough rejection
+sampling on the DF of the model. For each radius sampled above, a max on the DF is computed
+numerically to achieve this.
+
+The masses of all the particles are simply given by M/N.
+
+The script will also print out a guiding estimate of a good softening length to use, given by the radius
+of a sphere which contains on average one particle at the number density at the softening length a.
+
+NOTE: the positions of the particles produced by this script are centered at the origin, and must
+be shifted by some distance in x,y,z to be accepted by SWIFT. This will depend on the box size and
+the bound parameter described above.
+
+## USAGE ##
+Set the parameters described above as desired in plummerIC.py
+Run plummerIC.py to produce .hdf5 initial conditions
+Set suitable parameters in params.yml
+Run SWIFT with --self-gravity on
+
+To display results, another script (plotdensity.py) is provided. The basic usage is
+python3 plotdensity.py [path to snapshot files to plot (max 20)] -a [softening value] -M [mass value] -shift [shift applied in params.yml]
+
+A run with default values, 10k particles and subsequent density plotting can be automated with the provided run.sh script.
diff --git a/examples/GravityTests/Plummer_Selfgravitating/params.yml b/examples/GravityTests/Plummer_Selfgravitating/params.yml
new file mode 100755
index 0000000000000000000000000000000000000000..ae86dde1f25aaf521e3411f84a99696f2d650448
--- /dev/null
+++ b/examples/GravityTests/Plummer_Selfgravitating/params.yml
@@ -0,0 +1,38 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.988e+43 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.086e+21 # kpc
+  UnitVelocity_in_cgs: 1e5       # km / s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            2.0     # The end time of the simulation (in internal units).
+  dt_min:              1e-11   # The minimal time-step size of the simulation (in internal units).
+  dt_max:              1e-3    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                           0.002      # Constant dimensionless multiplier for time integration.
+  MAC:                           geometric  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+  epsilon_fmm:                   0.001      # Tolerance parameter for the adaptive multipole acceptance criterion.
+  theta_cr:                      0.7        # Opening angle for the purely gemoetric criterion.
+  max_physical_DM_softening:     0.005      # Physical softening length (in internal units).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snap/snapshot  # Common part of the name of output files
+  time_first:          0.      # Time of the first output (in internal units)
+  delta_time:          2e-1    # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          plummer.hdf5 # The file to read
+  shift:              [2.,2.,2.]
+  periodic:           0
diff --git a/examples/GravityTests/Plummer_Selfgravitating/plotdensity.py b/examples/GravityTests/Plummer_Selfgravitating/plotdensity.py
new file mode 100644
index 0000000000000000000000000000000000000000..27903fd1412aa9ff6838a90b38c276d87b7207ef
--- /dev/null
+++ b/examples/GravityTests/Plummer_Selfgravitating/plotdensity.py
@@ -0,0 +1,107 @@
+import numpy as np
+from matplotlib import pyplot as plt
+import argparse
+from scipy.optimize import curve_fit
+import h5py
+
+# Parse user input
+parser = argparse.ArgumentParser(
+    description="Plot multiple density profiles against theoretical prediction"
+)
+parser.add_argument("files", nargs="+", help="snapshot files to be imaged")
+parser.add_argument("--notex", action="store_true", help="Flag to not use LaTeX markup")
+parser.add_argument(
+    "-a", type=float, default=0.1, help="Softening length of theoretical model"
+)
+parser.add_argument(
+    "-M", type=float, default=1.0e-5, help="Total Mass of theoretical model"
+)
+parser.add_argument("-Rmin", type=float, default=0.01, help="Min Radius")
+parser.add_argument("-Rmax", type=float, default=2.1, help="Max Radius")
+parser.add_argument(
+    "-nbsamples", type=int, default=64, help="Number of radii to sample (bins)"
+)
+parser.add_argument(
+    "-shift", type=float, default=2.0, help="Shift applied to particles in params.yml"
+)
+
+args = parser.parse_args()
+fnames = args.files
+
+# Limit number of snapshots to plot
+if len(fnames) > 20:
+    raise ValueError(
+        "Too many ({:d}) files provided (cannot plot more than 20).".format(len(fnames))
+    )
+
+# Set parameters
+tex = not args.notex
+if tex:
+    plt.rcParams.update({"text.usetex": tex, "font.size": 14, "font.family": "serif"})
+else:
+    plt.rcParams.update({"font.size": 12})
+figsize = 7
+
+# Model Parameters (Mestel surface density)
+G = 4.299581e04
+rsp = np.logspace(np.log10(args.Rmin), np.log10(args.Rmax), args.nbsamples)
+
+
+def plummer_analytical(r):
+    return (
+        3.0
+        * args.M
+        / (4.0 * np.pi * args.a ** 3)
+        * (1.0 + r ** 2 / args.a ** 2) ** (-2.5)
+    )
+
+
+# Plot densities
+fig, ax = plt.subplots(figsize=(1.2 * figsize, figsize))
+for fname in fnames:
+    print(fname)
+    f = h5py.File(fname, "r")
+    pos = np.array(f["DMParticles"]["Coordinates"]) - args.shift
+    time = (
+        f["Header"].attrs["Time"][0]
+        * f["Units"].attrs["Unit time in cgs (U_t)"][0]
+        / 31557600.0e9
+    )
+    mass = np.array(f["DMParticles"]["Masses"])
+    x = pos[:, 0]
+    y = pos[:, 1]
+    r = np.sqrt(np.sum(pos ** 2, 1))
+
+    # Methods to compute density profile
+    def mass_ins(R):
+        return ((r < R) * mass).sum()
+
+    mass_ins_vect = np.vectorize(mass_ins)
+
+    def density(R):
+        return np.diff(mass_ins_vect(R)) / np.diff(R) / (4.0 * np.pi * R[1:] ** 2)
+
+    dens = density(rsp)
+    rs = rsp[1:]
+
+    # remove empty bins
+    c = dens > 0
+    dens = np.compress(c, dens)
+    rs = np.compress(c, rs)
+
+    # Plot
+    # ax.loglog(rsp[1:], density(rsp), "o", ms=1.7, label=r"$t=$ {:.3f} Gyr".format(time))
+    ax.plot(rs, dens, label=r"$t=$ {:.3f} Gyr".format(time))
+
+ax.plot(rsp, plummer_analytical(rsp), c="black", label="Analytical Prediction")
+ax.set_xlabel("r [kpc]")
+ax.legend()
+ax.loglog()
+ax.set_title(
+    r"Plummer Density Profile: $a = {:.1e}$ kpc, $M = {:.1e}$ M$_{{\odot}}$".format(
+        args.a, args.M * 1e10
+    )
+)
+plt.tight_layout()
+ax.set_ylabel(r"$\rho(r)$ [$M_{\odot}$ kpc$^{-3}$]")
+plt.savefig("density_profile.png")
diff --git a/examples/GravityTests/Plummer_Selfgravitating/plummerIC.py b/examples/GravityTests/Plummer_Selfgravitating/plummerIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..173caa7aaff0aa816509d74810e7ef3a1c42e799
--- /dev/null
+++ b/examples/GravityTests/Plummer_Selfgravitating/plummerIC.py
@@ -0,0 +1,297 @@
+################################################################################
+# Copyright (c) 2022 Patrick Hirling (patrick.hirling@epfl.ch)
+#
+# 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 numpy as np
+import scipy.special as sci
+from scipy.optimize import minimize
+import time
+import argparse
+from swiftsimio import Writer
+import unyt
+
+parser = argparse.ArgumentParser(
+    description="Generate discrete realization of Anisotropic Plummer Model"
+)
+
+# Parse Main Arguments
+parser.add_argument("-a", type=float, default=0.05, help="Softening Length (in kpc)")
+parser.add_argument(
+    "-M", type=float, default=1.0e-5, help="Total Mass of Model (in 10^10 solar masses)"
+)
+parser.add_argument(
+    "-q", type=float, default=1.0, help="Anisotropy parameter (between -inf and +2)"
+)
+parser.add_argument(
+    "-N", type=int, default=1e4, help="Number of particles in the simulation"
+)
+parser.add_argument(
+    "-bound", type=float, default=2.0, help="Upper bound of radius sampling"
+)
+parser.add_argument("-boxsize", type=float, default=4.0, help="Box Size of simulation")
+parser.add_argument(
+    "--noparallel",
+    action="store_false",
+    help="Do not use multiprocessing library (default: false)",
+)
+parser.add_argument(
+    "-nbthreads",
+    type=int,
+    default=6,
+    help="Number of threads to use for multiprocessing (default: 6)",
+)
+
+args = parser.parse_args()
+if args.q > 2.0:
+    raise ValueError(
+        "Anisotropy parameter must be smaller than 2.0 (q = {:.2f} given)".format(
+            args.q
+        )
+    )
+
+#### Parameters
+# Plummer Model
+G = 4.299581e04  # Gravitational constant [kpc / 10^10 M_s * (kms)^2]
+a = args.a  # Plummer softening length [kpc]
+M = args.M  # Total Mass [10^10 M_s]
+q = args.q  # Anisotropy Parameter (-inf,2]
+
+# IC File
+N = int(args.N)  # Number of Particles
+fname = "plummer.hdf5"  # Name of the ic file (dont forget .hdf5)
+pickle_ics = 0  # Optional: Pickle ics (pos,vel,mass) for future use
+periodic_bdry = 0
+
+# Parallelism for velocity sampling
+use_parallel = args.noparallel
+nb_threads = int(args.nbthreads)  # Set to None to use os.cpucount()
+
+### Print parameters
+print(
+    "\n############################################################################################# \n"
+)
+print("Generating Plummer ICs with the following parameters::\n")
+print(
+    "Softening length:                                  a = " + "{:.3e} kpc".format(a)
+)
+print(
+    "Total Mass:                                        M = "
+    + "{:.3e} Solar Masses".format(M * 1e10)
+)
+print("Anisotropy parameter:                              q = " + str(q))
+print("Output file:                                       " + fname)
+print("Number of particles:                               N = " + "{:.1e}".format(N))
+print(
+    "\n############################################################################################# \n"
+)
+
+# Estimate good softening length (density at origin, sphere s.t. contains on avg 1 particle)
+epsilon0 = a * N ** (-1.0 / 3.0)
+print(
+    "Recommended Softening length (times conventional factor): "
+    + "{:.4e}".format(epsilon0)
+    + " kpc"
+)
+
+### Generate Positions
+# Inverse transform sampling from analytical inverse of M(<r)
+def r_of_m(m, a, M):
+    return a * ((M / m) ** (2.0 / 3.0) - 1) ** (-1.0 / 2.0)
+
+
+m_rand = M * np.random.uniform(0.0, 1.0, N)
+r_rand = r_of_m(m_rand, a, M)
+phi_rand = np.random.uniform(0.0, 2 * np.pi, N)
+theta_rand = np.arccos(np.random.uniform(-1.0, 1.0, N))
+
+x = r_rand * np.sin(theta_rand) * np.cos(phi_rand)
+y = r_rand * np.sin(theta_rand) * np.sin(phi_rand)
+z = r_rand * np.cos(theta_rand)
+
+X = np.array([x, y, z]).transpose()
+
+### Generate Velocities
+# Sample the DF of the Plummer model at the radii generated above (give vr,vt)
+# Dejonghe (1987) Anisotropic Distribution Functionl
+normalization = 3.0 * sci.gamma(6.0 - q) / (2.0 * (2.0 * np.pi) ** (2.5))
+units = G ** (q - 5.0) * M ** (q - 4.0) * a ** (2.0 - q)
+
+
+def H(E, L):
+    x = L ** 2 / (2.0 * E * a ** 2)
+    if x <= 1.0:
+        return 1.0 / (sci.gamma(4.5 - q)) * sci.hyp2f1(0.5 * q, q - 3.5, 1.0, x)
+    else:
+        return (
+            1.0
+            / (sci.gamma(1.0 - 0.5 * q) * sci.gamma(4.5 - 0.5 * q) * (x ** (0.5 * q)))
+            * sci.hyp2f1(0.5 * q, 0.5 * q, 4.5 - 0.5 * q, 1.0 / x)
+        )
+
+
+def Fq(E, L):
+    if E < 0.0:
+        return 0.0
+    else:
+        return normalization * units * E ** (3.5 - q) * H(E, L)
+
+
+# Helper Functions
+# Total energy: E = phi - 1/2v^2. relative potential: psi = -phi. Relative energy: -E = psi - 1/2v^2
+def relative_potential(r):  # = - phi
+    return G * M / np.sqrt(r ** 2 + a ** 2)
+
+
+def relative_energy(r, vr, vt):
+    return relative_potential(r) - 0.5 * (vr ** 2 + vt ** 2)
+
+
+# N.B: Angular momentum: L = r * vt
+
+# Convenience Function for scipy.minimize negative of Fq*vt, indep. vars passed as array
+def Fq_tomin(v, r):
+    return -Fq(relative_energy(r, v[0], v[1]), r * v[1]) * v[1]
+
+
+# Find max of DF at given radius
+def fmax(r, vmax):
+    # Constraint function (defined locally, dep on r)
+    def vel_constr2(v):
+        return vmax ** 2 - v[0] ** 2 - v[1] ** 2
+
+    # Initial Guess
+    v0 = [0.1 * vmax, 0.2 * vmax]
+
+    # Constraint Object
+    cons = {"type": "ineq", "fun": vel_constr2}
+
+    # Args
+    args = (r,)
+
+    # Minimize through scipy.optimize.minimize
+    # fm = minimize(Fq_tomin,v0,constraints=cons,method = 'COBYLA',args=args)
+    fm = minimize(
+        Fq_tomin,
+        v0,
+        constraints=cons,
+        method="SLSQP",
+        args=args,
+        bounds=[(0, vmax), (0, vmax)],
+    )
+
+    # Min of negative df == max of df
+    return -fm.fun
+
+
+# Sample vr,vt from DF at given Radius through rejection sampling
+def sample_vel(r):
+    # Compute max velocity (equiv. condition for E>=0)
+    vmax = np.sqrt(2.0 * relative_potential(r))
+    # Compute max of DF at this radius
+    fm = 1.1 * fmax(r, vmax)  # 1.1 factor to be sure to include max
+    while True:
+        # Select candidates for vr,vt based on max full velocity
+        while True:
+            vr = np.random.uniform(0.0, vmax)
+            vt = np.random.uniform(0.0, vmax)
+            if vr ** 2 + vt ** 2 <= vmax ** 2:
+                break
+        # Rejection Sampling on Fq
+        f = np.random.uniform(0.0, fm)
+        if Fq(relative_energy(r, vr, vt), r * vt) * vt >= f:
+            return vr, vt
+
+
+print("Sampling velocities...")
+ti = time.time()
+if use_parallel:
+    from multiprocessing import Pool
+
+    with Pool(nb_threads) as p:
+        vels = np.array(p.map(sample_vel, r_rand)).transpose()
+else:
+    from tqdm import tqdm
+
+    vels = np.empty((2, N))
+    for j, r in enumerate(tqdm(r_rand)):
+        vels[:, j] = sample_vel(r)
+tf = time.time()
+print("Sampling took " + "{:.2f}".format(tf - ti) + " seconds.")
+
+# Convert to Cartesian
+# First: project vt on e_theta, e_phi with random orientation
+alph = np.random.uniform(0, 2 * np.pi, N)
+sgn = np.random.choice([-1, 1], size=N)
+vphi = vels[1] * np.cos(alph)
+vtheta = vels[1] * np.sin(alph)
+# project vr on e_r (random sign)
+vr = sgn * vels[0]
+
+# Convert Spherical to cartesian coordinates
+v_x = (
+    np.sin(theta_rand) * np.cos(phi_rand) * vr
+    + np.cos(theta_rand) * np.cos(phi_rand) * vtheta
+    - np.sin(phi_rand) * vphi
+)
+v_y = (
+    np.sin(theta_rand) * np.sin(phi_rand) * vr
+    + np.cos(theta_rand) * np.sin(phi_rand) * vtheta
+    + np.cos(phi_rand) * vphi
+)
+v_z = np.cos(theta_rand) * vr - np.sin(theta_rand) * vtheta
+
+# Create velocity array
+V = np.array([v_x, v_y, v_z]).transpose()
+
+### Generate masses
+
+m = M / N * np.ones(N)
+
+### Exclude extreme outliers
+idx = r_rand < args.bound
+X = X[idx]
+V = V[idx]
+m = m[idx]
+new_N = len(m)
+
+### Write to hdf5
+galactic_units = unyt.UnitSystem(
+    "galactic",
+    unyt.kpc,
+    unyt.unyt_quantity(1e10, units=unyt.Solar_Mass),
+    unyt.unyt_quantity(1.0, units=unyt.s * unyt.Mpc / unyt.km).to(unyt.Gyr),
+)
+wrt = Writer(galactic_units, args.boxsize * unyt.kpc)
+wrt.dark_matter.coordinates = X * unyt.kpc
+wrt.dark_matter.velocities = V * (unyt.km / unyt.s)
+wrt.dark_matter.masses = m * 1e10 * unyt.msun
+wrt.dark_matter.particle_ids = np.arange(new_N)
+wrt.write(fname)
+print("Writing IC file...")
+
+### Optional: Pickle ic arrays for future use
+if pickle_ics:
+    import pickle as pkl
+
+    with open("X.pkl", "wb") as f:
+        pkl.dump(X[:], f)
+
+    with open("V.pkl", "wb") as f:
+        pkl.dump(V[:], f)
+
+    with open("M.pkl", "wb") as f:
+        pkl.dump(m[:], f)
diff --git a/examples/GravityTests/Plummer_Selfgravitating/run.sh b/examples/GravityTests/Plummer_Selfgravitating/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..cabcfd013b8d627caefd7eb1d7e1820199899659
--- /dev/null
+++ b/examples/GravityTests/Plummer_Selfgravitating/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+if [ ! -e plummer.hdf5 ]
+then
+    echo "Generating initial conditions for Plummer example..."
+    python3 plummerIC.py -a 0.1 -N 65536
+fi
+
+# create output directory
+if [ ! -e snap ]
+then
+  mkdir snap
+fi
+
+../../swift --self-gravity --threads=8 params.yml 2>&1 | tee output.log
+
+echo "Plotting results..."
+# If params.yml is left at default values, should produce 10 snapshots
+python3 plotdensity.py snap/snapshot_*.hdf5