diff --git a/src/struphy/io/options.py b/src/struphy/io/options.py
index 6ecfff7c7..ece2ea02f 100644
--- a/src/struphy/io/options.py
+++ b/src/struphy/io/options.py
@@ -52,6 +52,7 @@ class LiteralOptions:
 
     # fields background
     BackgroundTypes = Literal["LogicalConst", "FluidEquilibrium"]
+    KineticDimensionsToPlot = Literal["e1", "e2", "e3", "v1", "v2", "v3"]
 
     # models
     ModelTypes = Literal["Toy", "Kinetic", "Fluid", "Hybrid"]
diff --git a/src/struphy/kinetic_background/base.py b/src/struphy/kinetic_background/base.py
index a223ab7d6..fc09444d1 100644
--- a/src/struphy/kinetic_background/base.py
+++ b/src/struphy/kinetic_background/base.py
@@ -4,9 +4,14 @@
 from typing import Callable
 
 import cunumpy as xp
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from matplotlib.colors import Normalize
 
 from struphy.fields_background.base import FluidEquilibriumWithB
+from struphy.geometry.base import Domain
 from struphy.initial.base import Perturbation
+from struphy.io.options import LiteralOptions
 
 
 class KineticBackground(metaclass=ABCMeta):
@@ -109,6 +114,248 @@ def __call__(self, *args):
         """
         pass
 
+    def plot_density_profile(
+        self,
+        dim_1: LiteralOptions.KineticDimensionsToPlot = "e1",
+        dim_2: LiteralOptions.KineticDimensionsToPlot | None = None,
+        v_lim: float = 5.0,
+        resol: int = 100,
+        integrate_resol: int = 50,
+        logical_coord: tuple[float] = (0.5, 0.5, 0.5),
+        in_physical: bool = False,
+        domain: Domain | None = None,
+        proj_axis: tuple[float,] = (0, 1),
+        plot_3D: bool = False,
+        title: str | None = None,
+        use_mu=False,
+        equil=None,
+    ):
+        """
+        Plots the density profile (slices) of the phase space background distribution. The slice can either be 1D or 2D, in logical or in Cartesian coordinates.
+
+        Parameters
+        ----------
+        dim_1, dim_2 : LiteralOptions.KineticDimensionsToPlot = ["e1","e2","e3","v1","v2","v3"]
+            The axes used in the projection, they refere to logical space axes. If dim_2 is not defined the projection is 1D, it is 2D if dim_2 is attributed.
+
+        v_lim : float = 5.0
+            Limit value of the velocity axes.
+
+        resol : int = 100
+            Resolution along each axes
+
+        integrate_resol : int = 10
+            Resolution along not used velocity axes. The density is reduced (with a maximum function) over these axes before being plotted. High values (>50) may require much memory.
+
+        logical_coord : tuple[float] = (0.5, 0.5, 0.5)
+            Refere to the default coordinate (in logical space) attributed to each axe which is not used in the projection.
+
+        in_physical : bool = False
+            Specify if the result is plotted in logical coordinates or in Cartesian coordinates, has a effect in 2D plotting. If True, you must specify a domain.
+
+        domain : Domain | None = None
+            Domain used to plot the density if in_physical=True.
+
+        proj_axis : tuple[float] = (0,1)
+            Axes of the Cartesian coordinates used to plot the density: 0->x, 1->y, 2->z.
+            I you do not see the density profile in 2D, you may change these axes.
+
+        plot_3D : bool = False
+            Plots a surface in a 3D environment. Only for physical projection.
+        """
+        integrate_resol = min(integrate_resol, 100)  # to avoid memory issues
+        if in_physical:
+            if not (dim_1 in ["e1", "e2", "e3"] and dim_2 in ["e1", "e2", "e3"]):
+                AssertionError(
+                    'To perform a plot in physical space you must use two space axes (dim_1, dim_2 in ["e1","e2","e3"]).'
+                )
+        if plot_3D and not in_physical:
+            AssertionError("To perform a 3D plot you must plot in physical space (activate in_physical).")
+        assert 0 <= proj_axis[0] < proj_axis[1] < 3
+        if dim_2 is None:
+            if dim_1 == "e1":
+                axe_to_plot = 0
+            elif dim_1 == "e2":
+                axe_to_plot = 1
+            elif dim_1 == "e3":
+                axe_to_plot = 2
+            elif dim_1 == "v1":
+                axe_to_plot = 3
+            elif dim_1 == "v2":
+                axe_to_plot = 4
+            elif dim_1 == "v3":
+                axe_to_plot = 5
+            else:
+                AssertionError("dim_1argument must match an exiting dimension")
+            if axe_to_plot - 3 > self.vdim:
+                AssertionError("Coordinate " + dim_1 + " does not exist with this background")
+            linspace_space = xp.array([0.0])
+            integrate_linspace_vel = xp.linspace(0.0, v_lim, integrate_resol)
+            if axe_to_plot < 3:
+                plot_linspace = xp.linspace(0.0, 1.0, resol)
+            else:
+                plot_linspace = xp.linspace(-v_lim, v_lim, resol)
+            tabs = 3 * [linspace_space] + self.vdim * [integrate_linspace_vel]
+            for i in range(3):
+                tabs[i][0] = logical_coord[i]
+            tabs[axe_to_plot] = plot_linspace
+            etas = xp.abs(xp.meshgrid(*tabs, indexing="ij"))
+            total_density = self(*etas)
+            if use_mu and axe_to_plot == 4:
+                B_norm_tab = equil.absB0(
+                    etas[0][tuple([slice(None), slice(None), slice(None)] + self.vdim * [0])],
+                    etas[1][tuple([slice(None), slice(None), slice(None)] + self.vdim * [0])],
+                    etas[2][tuple([slice(None), slice(None), slice(None)] + self.vdim * [0])],
+                )
+                B_norm_tab = xp.broadcast_to(
+                    B_norm_tab[
+                        tuple(
+                            [
+                                ...,
+                            ]
+                            + self.vdim * [None]
+                        )
+                    ],
+                    etas[0].shape,
+                )
+                total_density *= B_norm_tab / etas[4]
+                plot_linspace = plot_linspace**2
+            axes_to_integrate = [i for i in range(3 + self.vdim)]
+            axes_to_integrate.remove(axe_to_plot)
+            total_density = xp.mean(total_density, tuple(axes_to_integrate))
+            fig, ax = plt.subplots(1, 1)
+            ax.plot(plot_linspace, total_density)
+            ax.set_xlabel(dim_1)
+            ax.set_ylabel("density")
+            ax.set_title("background profile")
+            plt.show(block=True)
+        else:
+            if dim_1 == "e1":
+                axe_to_plot1 = 0
+            elif dim_1 == "e2":
+                axe_to_plot1 = 1
+            elif dim_1 == "e3":
+                axe_to_plot1 = 2
+            elif dim_1 == "v1":
+                axe_to_plot1 = 3
+            elif dim_1 == "v2":
+                axe_to_plot1 = 4
+            elif dim_1 == "v3":
+                axe_to_plot1 = 5
+            else:
+                AssertionError("dim_1 argument must match an existing dimension")
+            if dim_2 == "e1":
+                axe_to_plot2 = 0
+            elif dim_2 == "e2":
+                axe_to_plot2 = 1
+            elif dim_2 == "e3":
+                axe_to_plot2 = 2
+            elif dim_2 == "v1":
+                axe_to_plot2 = 3
+            elif dim_2 == "v2":
+                axe_to_plot2 = 4
+            elif dim_2 == "v3":
+                axe_to_plot2 = 5
+            else:
+                AssertionError("dim_2 argument must match an existing dimension")
+            if axe_to_plot2 == axe_to_plot1:
+                AssertionError("You must specify different dimensions for dim_1 and dim_2")
+            integrate_linspace_vel = xp.linspace(0.0, v_lim, integrate_resol)
+            tabs = [xp.array([logical_coord[i]]) for i in range(3)] + self.vdim * [integrate_linspace_vel]
+            if axe_to_plot1 < 3:
+                plot_linspace1 = xp.linspace(0.0, 1.0, resol)
+            elif axe_to_plot1 != 4 or (not use_mu):
+                plot_linspace1 = xp.linspace(-v_lim, v_lim, resol)
+            else:
+                plot_linspace1 = xp.linspace(0.0, xp.sqrt(v_lim), resol)
+            if axe_to_plot2 < 3:
+                plot_linspace2 = xp.linspace(0.0, 1.0, resol)
+            elif axe_to_plot2 != 4 or (not use_mu):
+                plot_linspace2 = xp.linspace(-v_lim, v_lim, resol)
+            else:
+                plot_linspace2 = xp.linspace(0.0, xp.sqrt(v_lim), resol)
+            tabs[axe_to_plot1] = plot_linspace1
+            tabs[axe_to_plot2] = plot_linspace2
+            etas = xp.meshgrid(*tabs, indexing="ij")
+            if in_physical:
+                physical_coords = domain(*tabs[:3])
+                physical_coords = list(physical_coords)
+                for i in range(3):
+                    physical_coords[i] = physical_coords[i][tuple([...] + self.vdim * [None])]
+                    physical_coords[i] = xp.broadcast_to(physical_coords[i], etas[0].shape)
+                for i in range(self.vdim):
+                    physical_coords.append(etas[i + 3])
+                total_density = self(*xp.abs(etas))
+            else:
+                physical_coords = list(etas)
+                total_density = self(*xp.abs(etas))
+            if use_mu:
+                if axe_to_plot1 == 4 or axe_to_plot2 == 4:
+                    B_norm_tab = equil.absB0(
+                        etas[0][tuple([slice(None), slice(None), slice(None)] + self.vdim * [0])],
+                        etas[1][tuple([slice(None), slice(None), slice(None)] + self.vdim * [0])],
+                        etas[2][tuple([slice(None), slice(None), slice(None)] + self.vdim * [0])],
+                    )
+                    B_norm_tab = xp.broadcast_to(
+                        B_norm_tab[
+                            tuple(
+                                [
+                                    ...,
+                                ]
+                                + self.vdim * [None]
+                            )
+                        ],
+                        etas[0].shape,
+                    )
+                    total_density *= B_norm_tab / etas[4]
+                    physical_coords[4] = physical_coords[4] ** 2
+            axes_to_integrate = [i for i in range(3 + self.vdim)]
+            axes_to_integrate.remove(axe_to_plot1)
+            axes_to_integrate.remove(axe_to_plot2)
+            total_density = xp.mean(total_density, tuple(axes_to_integrate))
+            id_dim = [0] * len(etas)
+            id_dim[axe_to_plot1] = slice(None)
+            id_dim[axe_to_plot2] = slice(None)
+            if not plot_3D:
+                if in_physical:
+                    X = physical_coords[proj_axis[0]][tuple(id_dim)]
+                    Y = physical_coords[proj_axis[1]][tuple(id_dim)]
+                else:
+                    X = physical_coords[axe_to_plot1][tuple(id_dim)]
+                    Y = physical_coords[axe_to_plot2][tuple(id_dim)]
+                fig, ax = plt.subplots()
+                for_color = ax.pcolor(X, Y, total_density)
+                if in_physical:
+                    ax.set_xlabel(["x", "y", "z"][proj_axis[0]])
+                    ax.set_ylabel(["x", "y", "z"][proj_axis[1]])
+                else:
+                    if use_mu:
+                        if dim_1 == "v2":
+                            dim_1 = "mu"
+                        if dim_2 == "v2":
+                            dim_2 = "mu"
+                    ax.set_xlabel(dim_1)
+                    ax.set_ylabel(dim_2)
+            else:
+                norm = Normalize(total_density.min(), total_density.max() + 0.01)
+                colors = cm.viridis(norm(total_density))
+                fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+                for_color = ax.plot_surface(
+                    X=physical_coords[0][tuple(id_dim)],
+                    Y=physical_coords[1][tuple(id_dim)],
+                    Z=physical_coords[2][tuple(id_dim)],
+                    facecolors=colors,
+                )
+                ax.set_xlabel("x")
+                ax.set_ylabel("y")
+                ax.set_zlabel("z")
+            if title is None:
+                ax.set_title(f"Density in ({dim_1}, {dim_2}) space")
+            else:
+                ax.set_title(title)
+            fig.colorbar(for_color)
+            plt.show(block=True)
+
     def __add__(self, other_f0):
         return SumKineticBackground(self, other_f0)
 
diff --git a/src/struphy/kinetic_background/maxwellians.py b/src/struphy/kinetic_background/maxwellians.py
index d9e34dcab..ed19be5c6 100644
--- a/src/struphy/kinetic_background/maxwellians.py
+++ b/src/struphy/kinetic_background/maxwellians.py
@@ -6,7 +6,9 @@
 
 from struphy.fields_background.base import FluidEquilibriumWithB
 from struphy.fields_background.equils import set_defaults
+from struphy.geometry.base import Domain
 from struphy.initial.base import Perturbation
+from struphy.io.options import LiteralOptions
 from struphy.kinetic_background.base import Maxwellian
 
 
@@ -300,6 +302,40 @@ def vth(self, eta1, eta2, eta3):
         out += [self._evaluate_moment(eta1, eta2, eta3, name="vth_perp")]
         return [ou * mom_fac for ou, mom_fac in zip(out, self.moment_factors["vth"])]
 
+    def plot_density_profile(
+        self,
+        dim_1: LiteralOptions.KineticDimensionsToPlot = "e1",
+        dim_2: LiteralOptions.KineticDimensionsToPlot | None = None,
+        v_lim: float = 5.0,
+        resol: int = 100,
+        integrate_resol: int = 10,
+        logical_coord: tuple[float] = (0.5, 0.5, 0.5),
+        in_physical: bool = False,
+        domain: Domain | None = None,
+        proj_axis: tuple[float,] = (0, 1),
+        plot_3D: bool = False,
+        title: str | None = None,
+        use_mu: bool = False,
+        equil: FluidEquilibriumWithB | None = None,
+    ):
+        if equil is None:
+            equil = self.equil
+        super().plot_density_profile(
+            dim_1,
+            dim_2,
+            v_lim,
+            resol,
+            integrate_resol,
+            logical_coord,
+            in_physical,
+            domain,
+            proj_axis,
+            plot_3D,
+            title,
+            use_mu=use_mu,
+            equil=equil,
+        )
+
 
 class CanonicalMaxwellian:
     r"""canonical Maxwellian distribution function.
diff --git a/src/struphy/kinetic_background/tests/test_base.py b/src/struphy/kinetic_background/tests/test_base.py
index df32b44d2..340d5740a 100644
--- a/src/struphy/kinetic_background/tests/test_base.py
+++ b/src/struphy/kinetic_background/tests/test_base.py
@@ -84,5 +84,60 @@ def test_kinetic_background_magics(show_plot=False):
         plt.show()
 
 
+def test_plotting_function():
+
+    import cunumpy as xp
+
+    from struphy import domains, equils, maxwellians
+
+    equil = equils.HomogenSlab(B0x=0.0, B0y=0.0, B0z=1.0)
+    equil.domain = domains.Cuboid()
+
+    # definition of test functions
+    l, m, n = 3, 4, 5
+
+    def n_init(*etas):
+        if len(etas) == 1:
+            e1, e2, e3 = etas[0][:, 0], etas[0][:, 1], etas[0][:, 1]
+        else:
+            assert len(etas) == 3
+            e1, e2, e3 = etas[0], etas[1], etas[2]
+        return 1 + 0.5 * xp.cos(2 * xp.pi * e1 * l) * xp.cos(2 * xp.pi * e2 * m) * xp.cos(2 * xp.pi * e3 * n)
+
+    def vth(*etas):
+        if len(etas) == 1:
+            e1, e2, e3 = etas[0][:, 0], etas[0][:, 1], etas[0][:, 1]
+        else:
+            assert len(etas) == 3
+            e1, e2, e3 = etas[0], etas[1], etas[2]
+        return 1 + 0.2 * xp.cos(2 * xp.pi * e1 * l) * xp.cos(2 * xp.pi * e2 * m) * xp.cos(2 * xp.pi * e3 * n)
+
+    # Testing with GyroMaxwellian2D:
+    background = maxwellians.GyroMaxwellian2D(n=(n_init, None), vth_para=(vth, None), vth_perp=(vth, None), equil=equil)
+    background.plot_density_profile("e1")
+    background.plot_density_profile("e2")
+    background.plot_density_profile("e3")
+    background.plot_density_profile("v1")
+    background.plot_density_profile("v2")
+    background.plot_density_profile("e1", "e2")
+    background.plot_density_profile("e1", "e2", domain=domains.HollowCylinder(), proj_axis=(0, 1), in_physical=True)
+    background.plot_density_profile("e1", "e2", domain=domains.HollowTorus(), proj_axis=(1, 2), in_physical=True)
+    background.plot_density_profile(
+        "e1", "e2", domain=domains.HollowTorus(), proj_axis=(0, 2), in_physical=True, plot_3D=True
+    )
+    background.plot_density_profile(
+        "e2", "e3", domain=domains.HollowTorus(), proj_axis=(1, 2), in_physical=True, plot_3D=True
+    )
+    background.plot_density_profile("v1", "v2")
+    background.plot_density_profile("v1", "v2", use_mu=True)
+
+    # Testing with Maxwellian3D:
+    background = maxwellians.Maxwellian3D(n=(n_init, None), vth1=(vth, None), vth2=(vth, None), vth3=(vth, None))
+    background.plot_density_profile("v1", "v2")
+    background.plot_density_profile("v1", "v3")
+    background.plot_density_profile("e1", "v3")
+
+
 if __name__ == "__main__":
-    test_kinetic_background_magics(show_plot=True)
+    # test_kinetic_background_magics(show_plot=True)
+    test_plotting_function()
diff --git a/tutorials/tutorial_03_test_particles.ipynb b/tutorials/tutorial_03_test_particles.ipynb
index ed3ab9c89..7b9ae89dc 100644
--- a/tutorials/tutorial_03_test_particles.ipynb
+++ b/tutorials/tutorial_03_test_particles.ipynb
@@ -443,6 +443,83 @@
    "cell_type": "markdown",
    "id": "29",
    "metadata": {},
+   "source": [
+    "Sometimes you need to load particles more uniformly. It is usefull to reach an homogeneous description of the plasma with a limited number of particles.\n",
+    "\n",
+    "It is performed using `loading=\"sobol_standard\"` or `loading=\"sobol_antithetic\"` in the `LoadingParameters` function:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "30",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "env_3 = EnvironmentOptions(sim_folder=\"sim_3\")\n",
+    "model_3 = Vlasov()\n",
+    "sim_3 = Simulation.spawn_sister(sim, env=env_3, model=model_3)\n",
+    "\n",
+    "# species parameters\n",
+    "loading_params_3 = LoadingParameters(Np=1000, spatial=\"disc\", loading=\"sobol_standard\")\n",
+    "\n",
+    "model_3.kinetic_ions.set_markers(\n",
+    "    loading_params=loading_params_3, weights_params=weights_params, boundary_params=boundary_params\n",
+    ")\n",
+    "model_3.kinetic_ions.set_sorting_boxes()\n",
+    "\n",
+    "model_3.kinetic_ions.set_save_data(n_markers=1.0)\n",
+    "\n",
+    "# propagator options\n",
+    "model_3.propagators.push_vxb.options = model_3.propagators.push_vxb.Options()\n",
+    "model_3.propagators.push_eta.options = model_3.propagators.push_eta.Options()\n",
+    "\n",
+    "# initial conditions (background + perturbation)\n",
+    "model_3.kinetic_ions.var.add_background(background)\n",
+    "\n",
+    "sim_3.run()\n",
+    "sim_3.pproc()\n",
+    "sim_3.load_plotting_data()\n",
+    "\n",
+    "fig = plt.figure(figsize=(15, 6))\n",
+    "\n",
+    "orbits_standard = sim_3.orbits.kinetic_ions\n",
+    "\n",
+    "plt.subplot(1, 3, 1)\n",
+    "plt.scatter(orbits[0, :, 0], orbits[0, :, 1], s=2.0)\n",
+    "circle1 = plt.Circle((0, 0), a2, color=\"k\", fill=False)\n",
+    "ax = plt.gca()\n",
+    "ax.add_patch(circle1)\n",
+    "ax.set_aspect(\"equal\")\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.title(\"sim_1: draw uniform in logical space\")\n",
+    "\n",
+    "plt.subplot(1, 3, 2)\n",
+    "plt.scatter(orbits_uni[0, :, 0], orbits_uni[0, :, 1], s=2.0)\n",
+    "circle2 = plt.Circle((0, 0), a2, color=\"k\", fill=False)\n",
+    "ax = plt.gca()\n",
+    "ax.add_patch(circle2)\n",
+    "ax.set_aspect(\"equal\")\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.title(\"sim_2: draw uniform on disc\")\n",
+    "\n",
+    "plt.subplot(1, 3, 3)\n",
+    "plt.scatter(orbits_standard[0, :, 0], orbits_standard[0, :, 1], s=2.0)\n",
+    "circle3 = plt.Circle((0, 0), a2, color=\"k\", fill=False)\n",
+    "ax = plt.gca()\n",
+    "ax.add_patch(circle3)\n",
+    "ax.set_aspect(\"equal\")\n",
+    "plt.xlabel(\"x\")\n",
+    "plt.ylabel(\"y\")\n",
+    "plt.title(\"sim_3: draw using sobol_standard\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "31",
+   "metadata": {},
    "source": [
     "## Reflecting boundary conditions \n",
     "\n",
@@ -453,7 +530,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "30",
+   "id": "32",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -465,7 +542,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "31",
+   "id": "33",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -484,7 +561,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "32",
+   "id": "34",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -497,7 +574,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "33",
+   "id": "35",
    "metadata": {},
    "source": [
     "We still have to set the propagator options and the initial conditions:"
@@ -506,7 +583,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "34",
+   "id": "36",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -523,7 +600,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "35",
+   "id": "37",
    "metadata": {},
    "source": [
     "We can now run the simulation, then post-process the data and plot the resulting orbits:"
@@ -532,7 +609,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "36",
+   "id": "38",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -542,7 +619,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "37",
+   "id": "39",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -552,7 +629,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "38",
+   "id": "40",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -561,7 +638,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "39",
+   "id": "41",
    "metadata": {},
    "source": [
     "Under `sim.orbits[<species_name>]` one finds a three-dimensional numpy array; the first index refers to the time step, the second index to the particle and the third index to the particle attribute. The first three attributes are the particle positions, followed by the velocities and the (initial and time-dependent) weights."
@@ -570,7 +647,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "40",
+   "id": "42",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -584,7 +661,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "41",
+   "id": "43",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -614,7 +691,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "42",
+   "id": "44",
    "metadata": {},
    "source": [
     "## Particles in a cylinder with a magnetic field\n",
@@ -638,7 +715,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "43",
+   "id": "45",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -650,7 +727,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "44",
+   "id": "46",
    "metadata": {},
    "source": [
     "In order to project the equilibrium on the spline basis for fast evaluation in the particle kernels, we need a Derham complex:"
@@ -659,7 +736,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "45",
+   "id": "47",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -669,7 +746,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "46",
+   "id": "48",
    "metadata": {},
    "source": [
     "Now we create the light-weight instance of the model and set the species options. We shall `remove` particles that hit the boundary in $\\eta_1$ (radial) direction:"
@@ -678,7 +755,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "47",
+   "id": "49",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -693,7 +770,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "48",
+   "id": "50",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -718,7 +795,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "49",
+   "id": "51",
    "metadata": {},
    "source": [
     "Now the usual procedure: run, post-process, load data and finally plot the orbits:"
@@ -727,7 +804,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "50",
+   "id": "52",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -737,7 +814,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "51",
+   "id": "53",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -747,7 +824,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "52",
+   "id": "54",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -757,7 +834,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "53",
+   "id": "55",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -770,7 +847,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "54",
+   "id": "56",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -798,7 +875,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "55",
+   "id": "57",
    "metadata": {},
    "source": [
     "## Particles in a Tokamak equilibrium\n",
@@ -809,7 +886,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "56",
+   "id": "58",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -821,7 +898,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "57",
+   "id": "59",
    "metadata": {},
    "source": [
     "Since `EQDSKequilibrium` is an `AxisymmMHDequilibrium`, which is a subclass of `CartesianMHDequilibrium`, we are free to choose any mapping for the simulation (e.g. a Cuboid for Cartesian coordinates). In order to be conforming to the boundary of the equilibrium, we shall choose the `Tokamak` mapping:"
@@ -830,7 +907,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "58",
+   "id": "60",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -843,7 +920,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "59",
+   "id": "61",
    "metadata": {},
    "source": [
     "The `Tokamak` domain is a `PoloidalSplineTorus`, hence\n",
@@ -894,7 +971,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "60",
+   "id": "62",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -912,7 +989,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "61",
+   "id": "63",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -927,7 +1004,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "62",
+   "id": "64",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -944,7 +1021,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "63",
+   "id": "65",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1001,7 +1078,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "64",
+   "id": "66",
    "metadata": {},
    "source": [
     "We need a Derham complex for the projection of the equilibrium onto the spline basis:"
@@ -1010,7 +1087,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "65",
+   "id": "67",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1024,7 +1101,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "66",
+   "id": "68",
    "metadata": {},
    "source": [
     "We aim to simulate 15000 time steps with a second-order splitting algorithm:"
@@ -1033,7 +1110,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "67",
+   "id": "69",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1042,7 +1119,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "68",
+   "id": "70",
    "metadata": {},
    "source": [
     "We now set up a simulation of 4 specific particle orbits in this equilibrium:"
@@ -1051,7 +1128,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "69",
+   "id": "71",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1070,7 +1147,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "70",
+   "id": "72",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1104,7 +1181,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "71",
+   "id": "73",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1114,7 +1191,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "72",
+   "id": "74",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1124,7 +1201,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "73",
+   "id": "75",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1134,7 +1211,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "74",
+   "id": "76",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1147,7 +1224,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "75",
+   "id": "77",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1172,7 +1249,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "76",
+   "id": "78",
    "metadata": {},
    "source": [
     "## Guiding-centers in a Tokamak equilibrium\n",
@@ -1183,7 +1260,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "77",
+   "id": "79",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1196,7 +1273,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "78",
+   "id": "80",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1208,7 +1285,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "79",
+   "id": "81",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1242,7 +1319,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "80",
+   "id": "82",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1300,7 +1377,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "81",
+   "id": "83",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1310,7 +1387,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "82",
+   "id": "84",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1320,7 +1397,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "83",
+   "id": "85",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1330,7 +1407,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "84",
+   "id": "86",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1343,7 +1420,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "85",
+   "id": "87",
    "metadata": {},
    "outputs": [],
    "source": [