diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml
index ba1ddca..655558d 100644
--- a/.github/workflows/docs.yaml
+++ b/.github/workflows/docs.yaml
@@ -3,6 +3,9 @@ on:
   push:
     branches:
       - main
+  pull_request:
+    branches:
+      - main
 
 permissions:
   contents: write
@@ -30,15 +33,22 @@ jobs:
         run: |
           for f in examples/*.py; do 
             echo "Processing $f..."
-            python "$f" -j 2
+            python "$f" -j 2 || exit 1
           done
 
       - name: Build Sphinx
         run: |
-          sphinx-build -b html docs/source docs/_build
+          sphinx-build -W -b html docs/source docs/_build/html
+
+      - name: Upload Documentation Artifact
+        uses: actions/upload-artifact@v4
+        with:
+          name: docs-preview
+          path: ./docs/_build/html
 
       - name: Deploy to GitHub Pages
+        if: github.event_name == 'push' && github.ref == 'refs/heads/main'
         uses: peaceiris/actions-gh-pages@v3
         with:
           github_token: ${{ secrets.GITHUB_TOKEN }}
-          publish_dir: ./docs/_build
\ No newline at end of file
+          publish_dir: ./docs/_build/html
\ No newline at end of file
diff --git a/README.md b/README.md
index 353d54d..e70c7a6 100644
--- a/README.md
+++ b/README.md
@@ -9,22 +9,22 @@ The module is built around 4 components:
 +  **`RunContext`**: Handles IO, directory creation, and data discovery. Detects simulation geometry, dimensions, and available fields.
 +  **`PhysicsProcessor`**: Performs mathematical transformations. Handles grid conversions (e.g., converting internal coordinates to Cartesian $x, z$ for plotting), applies zooms, and computes derived quantities.
 +  **`SliceRenderer`**: Matplotlib engine. Manages multi-panel layouts, colorbars (Log, Linear, TwoSlope), streamlines, contours.
-+  **`Pipeline`**: The coordinator:
-    *   Pre-calculates global bounds across all frames.
-    *   Distributes rendering tasks across multiple CPU cores.
-    *   Manages the temporal evolution of SpaceTime heatmaps and particle data.
++  **`Pipeline`**: Coordinates the detection, processing, and rendering of the simulation data. Computes the bounds and distributes with `multiprocessing`.
 
 ## Supported quantity types
 
 Users define what to plot by passing lists of "Quantity" objects to the Pipeline:
 
-*   **`MapMovie2D`**: For $f(x, z, t)$ fields, 2D pcolormesh plots generated for every frame. Supports:
+*   **`MapMovie2D`**: For $f(x, z, t)$ fields that will be rendered as a heatmap (pcolormesh) animation. Supports:
     *   `streamlines`: Vector overlays (e.g., velocity fields).
-    *   `contours`: Scalar overlays (e.g., density levels).
+    *   `contours`: draw contour lines over the pcolormesh. (e.g., density levels).
     *   `compute`: Custom functions to calculate new fields on the fly.
-*   **`LineMovie1D`**: For $f(x,t)$ fields, 1D plot generated for every frame (e.g., radial profile over time).
-*   **`SpaceTimeHeatmap`**: For $f(x,t)$ fields, generates a space-time heapmap. Supports `ref_function` to overlay analytical trajectories.
-*   **`PartQuantity`**: Tracks Lagrangian particle properties (like `PART_X1`) over time. Supports `ref_function` to overlay analytical functions.
+*   **`LineMovie1D`**: For $f(x,t)$ fields, renders a line plot that updates every frame. (e.g., radial profile over time).
+*   **`SpaceTimeHeatmap`**: For $f(x,t)$ fields, renders a space-time heatmap. Supports `ref_function` to overlay analytical trajectories.
+*   **`PartQuantity`**: Tracks Lagrangian particle properties (like `PART_X1`) over time. 
+
+`SpaceTimeHeatmap` and `PartQuantity` supports `ref_function` to overlay analytical functions.
+One can also plot a particle quantity over a spacetime heatmap with `trace_over`
 
 ## Usage
 
diff --git a/docs/source/_static/.gitkeep b/docs/source/_static/.gitkeep
new file mode 100644
index 0000000..e69de29
diff --git a/docs/source/api.rst b/docs/source/api.rst
index ee86958..cd98240 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -7,8 +7,6 @@ RunContext
 ----------
 .. autoclass:: RunContext
    :members:
-   :undoc-members:
-   :show-inheritance:
 
 Pipeline
 --------
@@ -19,6 +17,10 @@ Pipeline
 Quantity Types
 --------------
 .. autoclass:: MapMovie2D
-   :members:
+   :show-inheritance:
+.. autoclass:: LineMovie1D
+   :show-inheritance:
 .. autoclass:: SpaceTimeHeatmap
-   :members:
+   :show-inheritance:
+.. autoclass:: PartQuantity
+   :show-inheritance:
diff --git a/docs/source/examples.rst b/docs/source/examples.rst
index dc2a7ea..09fb191 100644
--- a/docs/source/examples.rst
+++ b/docs/source/examples.rst
@@ -1,7 +1,7 @@
 Examples Gallery
 ================
 
-Select an example below to explore the capabilities of the pipeline.
+Here are three most basic examples to use the Pipeline, depending on the nature of the quantity: 1D, 2D, or particles.
 
 .. grid:: 1 2 2 3
    :gutter: 3
@@ -32,4 +32,17 @@ Select an example below to explore the capabilities of the pipeline.
 
    examples/run_1D
    examples/run_2D
-   examples/run_particles
\ No newline at end of file
+   examples/run_particles
+   examples/run_particles_over
+
+Here are more complex examples
+
+.. grid:: 1 2 2 3
+   :gutter: 3
+
+   .. grid-item-card:: Particle over a heatmap
+      :link: examples/run_particles_over
+      :link-type: doc
+      :img-top: ../../examples/data_test/frames/particles_over_test/particles_over_test_spacetimeheatmap.png
+
+      Superimpose a particle trajectory over a :math:`(x,t)` heatmap
\ No newline at end of file
diff --git a/docs/source/examples/run_1D.rst b/docs/source/examples/run_1D.rst
index 1081347..b77215f 100644
--- a/docs/source/examples/run_1D.rst
+++ b/docs/source/examples/run_1D.rst
@@ -41,6 +41,8 @@ On the heatmap, one can add a reference line (e.g., an analytical trajectory) by
    :language: python
    :linenos:
 
-.. image:: ../../../examples/data_test/frames/1D_test_withref/1D_test_spacetimeheatmap.png
+.. figure:: ../../../examples/data_test/frames/1D_test_withref/1D_test_spacetimeheatmap.png
    :width: 80%
-   :align: center
\ No newline at end of file
+   :align: center
+
+   1D_test_spacetimeheatmap.png
\ No newline at end of file
diff --git a/docs/source/examples/run_particles.rst b/docs/source/examples/run_particles.rst
index 8b184a6..08a70a1 100644
--- a/docs/source/examples/run_particles.rst
+++ b/docs/source/examples/run_particles.rst
@@ -5,7 +5,9 @@ Particle quantities
    :language: python
    :linenos:
 
-.. image:: ../../../examples/data_test/frames/particles_test/particles_test_timeseries.png
+.. figure:: ../../../examples/data_test/frames/particles_test/particles_test_timeseries.png
     :width: 80%
     :align: center
 
+    particles_test_timeseries.png
+
diff --git a/docs/source/examples/run_particles_over.rst b/docs/source/examples/run_particles_over.rst
new file mode 100644
index 0000000..370d2e4
--- /dev/null
+++ b/docs/source/examples/run_particles_over.rst
@@ -0,0 +1,13 @@
+Particle quantities over heatmap
+================================
+
+.. literalinclude:: ../../../examples/run_particles_over.py
+   :language: python
+   :linenos:
+
+.. figure:: ../../../examples/data_test/frames/particles_over_test/particles_over_test_spacetimeheatmap.png
+    :width: 80%
+    :align: center
+    
+    particles_over_test_spacetimeheatmap.png
+
diff --git a/examples/data_test/inputs/particles_over_test.ini b/examples/data_test/inputs/particles_over_test.ini
new file mode 100644
index 0000000..f5c166e
--- /dev/null
+++ b/examples/data_test/inputs/particles_over_test.ini
@@ -0,0 +1,55 @@
+##
+[Grid]
+X1-grid    1  1.99      1  u  2.01
+X2-grid    1  -0.0125      1  u  0.0125
+X3-grid    1  -0.3  512    u  0.3
+
+[TimeIntegrator]
+CFL         0.5
+tstop       20
+first_dt    1.0e-5
+nstages     2
+
+[Hydro]
+solver    hllc
+csiso     userdef
+# viscosity    explicit  userdef
+
+[Dust]
+nSpecies         3
+drag             userdef  1   0.2   0.04    # St=1, 0.2, 0.04
+drag_feedback    no
+
+[Particles]
+count            per_proc  1
+stopping_time    constant  1.0
+ParticleMass     3e-3
+DustToGas        3e-3
+
+[Gravity]
+potential    central 
+Mcentral     1.0
+
+[Boundary]
+X1-beg    userdef
+X1-end    userdef
+X2-beg    periodic
+X2-end    periodic
+X3-beg    userdef
+X3-end    userdef
+
+[Setup]
+sigma0        0.125
+sigmaSlope    -0.5
+CsSlope       -0.5
+h0            0.05
+alpha         1.0e-4
+
+[Output]
+vtk    0.1
+dmp    1000.0
+log        1000
+dmp_dir    /home/dp316/dp316/dc-fang1/IdefixRuns/VerticalSettling/outputs/SettlingL_Tau
+vtk_dir    /home/dp316/dp316/dc-fang1/IdefixRuns/VerticalSettling/outputs/SettlingL_Tau/vtks
+# File produced automatically by a Perl script
+# Do not edit
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0000.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0000.vtk
new file mode 100644
index 0000000..73cc640
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0000.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0001.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0001.vtk
new file mode 100644
index 0000000..1c9d8aa
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0001.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0002.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0002.vtk
new file mode 100644
index 0000000..ca73628
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0002.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0003.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0003.vtk
new file mode 100644
index 0000000..8495991
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0003.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0004.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0004.vtk
new file mode 100644
index 0000000..35209a0
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0004.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0005.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0005.vtk
new file mode 100644
index 0000000..ad1f8ab
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0005.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0006.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0006.vtk
new file mode 100644
index 0000000..15aa7e4
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0006.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0007.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0007.vtk
new file mode 100644
index 0000000..77d9f24
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0007.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0008.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0008.vtk
new file mode 100644
index 0000000..44f1d61
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0008.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0009.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0009.vtk
new file mode 100644
index 0000000..26e5563
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0009.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0010.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0010.vtk
new file mode 100644
index 0000000..2a67e17
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0010.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0011.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0011.vtk
new file mode 100644
index 0000000..0e676a8
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0011.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0012.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0012.vtk
new file mode 100644
index 0000000..2ced5dc
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0012.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0013.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0013.vtk
new file mode 100644
index 0000000..eb82e53
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0013.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0014.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0014.vtk
new file mode 100644
index 0000000..56eb8a3
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0014.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0015.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0015.vtk
new file mode 100644
index 0000000..690ebd5
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0015.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0016.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0016.vtk
new file mode 100644
index 0000000..4590df4
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0016.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0017.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0017.vtk
new file mode 100644
index 0000000..c2f5d0c
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0017.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0018.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0018.vtk
new file mode 100644
index 0000000..7f0104d
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0018.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0019.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0019.vtk
new file mode 100644
index 0000000..6ba0d79
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0019.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/data.0020.vtk b/examples/data_test/outputs/particles_over_test/vtks/data.0020.vtk
new file mode 100644
index 0000000..4137cb8
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/data.0020.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0000.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0000.vtk
new file mode 100644
index 0000000..c082d56
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0000.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0001.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0001.vtk
new file mode 100644
index 0000000..7ef27c4
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0001.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0002.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0002.vtk
new file mode 100644
index 0000000..70182e9
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0002.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0003.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0003.vtk
new file mode 100644
index 0000000..58ccffa
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0003.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0004.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0004.vtk
new file mode 100644
index 0000000..b39ef78
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0004.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0005.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0005.vtk
new file mode 100644
index 0000000..8d182ec
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0005.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0006.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0006.vtk
new file mode 100644
index 0000000..30a06e6
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0006.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0007.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0007.vtk
new file mode 100644
index 0000000..595fc65
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0007.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0008.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0008.vtk
new file mode 100644
index 0000000..ddce872
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0008.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0009.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0009.vtk
new file mode 100644
index 0000000..9b84b3e
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0009.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0010.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0010.vtk
new file mode 100644
index 0000000..5d63b31
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0010.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0011.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0011.vtk
new file mode 100644
index 0000000..507251c
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0011.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0012.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0012.vtk
new file mode 100644
index 0000000..b30b881
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0012.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0013.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0013.vtk
new file mode 100644
index 0000000..d67d9c0
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0013.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0014.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0014.vtk
new file mode 100644
index 0000000..334793c
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0014.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0015.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0015.vtk
new file mode 100644
index 0000000..0869152
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0015.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0016.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0016.vtk
new file mode 100644
index 0000000..8a08a93
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0016.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0017.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0017.vtk
new file mode 100644
index 0000000..1511289
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0017.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0018.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0018.vtk
new file mode 100644
index 0000000..c7a1968
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0018.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0019.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0019.vtk
new file mode 100644
index 0000000..f4c7537
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0019.vtk differ
diff --git a/examples/data_test/outputs/particles_over_test/vtks/part.0020.vtk b/examples/data_test/outputs/particles_over_test/vtks/part.0020.vtk
new file mode 100644
index 0000000..a8ea58f
Binary files /dev/null and b/examples/data_test/outputs/particles_over_test/vtks/part.0020.vtk differ
diff --git a/examples/run_particles_over.py b/examples/run_particles_over.py
new file mode 100644
index 0000000..705b40a
--- /dev/null
+++ b/examples/run_particles_over.py
@@ -0,0 +1,46 @@
+from idefix2python import RunContext, Pipeline, PartQuantity, SpaceTimeHeatmap
+import utilities
+from pathlib import Path
+
+projectPath = Path(__file__).parent / "data_test"
+task = "particles_over_test"
+# By default the vtks are expected to be in {projetPath}/{task}/outputs/vtks/
+# In this example, the vtks/ folder contains both part*.vtk and data*.vtk
+
+
+def analytical_trajectory(t):
+    z0 = 0.1
+    fluid = utilities.Fluid(0.05, -0.5, 0.125, -0.5, Stokes0=1, z0=z0)
+    return utilities.solve_2nd_order_ode(fluid.azSettling, z0, 0, t)
+
+
+analytical_trajectory.plot_kwargs = {"ls": "--", "color": "cyan", "lw": 2}
+
+z_part = PartQuantity(
+    "PART_X3",
+    r"$z^\mathrm{part}$",
+    plot_coords=[0, 0],
+)
+
+custom_partQuantities = [z_part]
+SpaceTimeHeatmaps = [
+    SpaceTimeHeatmap(
+        "Dust0_RHO",
+        r"$\rho^\mathrm{dust}$",
+        plot_coords=[0, 0],
+        trace_over=[z_part],
+        ref_function=analytical_trajectory,
+    )
+]
+runContext = RunContext(
+    task,
+    projectPath,
+    active_directions=[2],  # currently necessary for lagrangian particles.)
+)
+pipeline = Pipeline(
+    runContext,
+    partQuantities=custom_partQuantities,
+    spaceTimeHeatmaps=SpaceTimeHeatmaps,
+)
+
+pipeline.run()
diff --git a/examples/utilities.py b/examples/utilities.py
index 23dd93d..d51ad3d 100644
--- a/examples/utilities.py
+++ b/examples/utilities.py
@@ -25,7 +25,7 @@ def vK(r):
 
 
 class Fluid:
-    def __init__(self, cs0, csSlope, sigma0, sigmaSlope, Stokes0):
+    def __init__(self, cs0, csSlope, sigma0, sigmaSlope, Stokes0, z0=0.1):
         self.cs0 = cs0
         self.csSlope = csSlope
         self.sigma0 = sigma0
@@ -39,6 +39,8 @@ def __init__(self, cs0, csSlope, sigma0, sigmaSlope, Stokes0):
         # self.rhoSlope = -1.5
         self.Stokes0 = Stokes0
 
+        self.z0 = z0
+
     def eta(self, r):
         cs0 = self.cs0
         csSlope = self.csSlope
@@ -58,6 +60,26 @@ def vrDrift(self, r):
         st = self.Stokes(r)
         return self.eta(r) * vK(r) / (st + 1 / st)
 
+    def vzSettling_approx(self, z):
+        r = 2.0
+        st = self.Stokes(r)
+        OmegaK = r ** (-1.5)
+        return -OmegaK * st * z
+
+    def azSettling(self, z, vz, t):
+        r = 2.0
+        tstop = self.Stokes0
+        OmegaK = r ** (-1.5)
+        return -vz / tstop - z * OmegaK**2
+
+    # def z_drift(self, t):
+    #     r = 2.0
+    #     tstop = self.Stokes0
+    #     OmegaK = r ** (-1.5)
+    #     return self.z0 * np.exp(OmegaK*st*z**2/2)
+    #     z = self.z0 * np.exp(-t / (2 * tstop))
+    #     st = self.Stokes(r)
+
     # def eta(self, r):
     #     cs0 = self.cs0
     #     csSlope = self.csSlope
@@ -72,3 +94,19 @@ def vrDrift(self, r):
 
     # def vrDrift(self, r):
     #     return -self.eta(r) / np.sqrt(r) / (self.Stokes(r) + 1 / self.Stokes(r))
+
+
+def solve_2nd_order_ode(f, u0, du0, times):
+    """
+    u" = f(u',u,t)
+    """
+    u = [u0]
+    du = [du0]
+    dt = np.diff(times)
+    for i, t in enumerate(times[:-1]):
+        u_prev = u[-1]
+        du_prev = du[-1]
+
+        du += [du_prev + f(u_prev, du_prev, t) * dt[i]]
+        u += [u_prev + du_prev * dt[i]]
+    return u
diff --git a/idefix2python/core.py b/idefix2python/core.py
index 9ee8fa4..4325ca5 100644
--- a/idefix2python/core.py
+++ b/idefix2python/core.py
@@ -85,13 +85,29 @@ def set_times(self, times):
 
 
 class Data:
-    # Types of data that I can imagine:
-    # 2D Field
-    # 1D Field time graph
-    # Scalar time graph
+    """
+    Base class for all data quantities in the pipeline.
+
+    :param key: Unique identifier for the field.
+    :type key: str
+    :param symbol: Symbol for labels (e.g., r"$\rho$").
+    :type symbol: str
+    :param plot_coords: [row, col] position in the subplot grid, defaults to [0, 0].
+    :type plot_coords: list[int], optional
+    :param vmin: Minimum value for manual scaling, defaults to None.
+    :type vmin: float, optional
+    :param vmax: Maximum value for manual scaling, defaults to None.
+    :type vmax: float, optional
+    :param \**kwargs:
+        * **title** (str): Custom title for the plot. Defaults to `symbol`.
+        * **id** (str): Unique ID to distinguish instances of the same field nature.
+        * **scale** (str): Scaling type, e.g., 'linear' or 'log'.
+        * **ref_function** (callable): Analytical function for comparison.
+    """
+
     timeline_instances = count(1)
 
-    def __init__(self, key, symbol, plot_coords, vmin=None, vmax=None, **kwargs):
+    def __init__(self, key, symbol, plot_coords=[0, 0], vmin=None, vmax=None, **kwargs):
         self.key = key
         self.symbol = symbol
         self.plot_coords = plot_coords
@@ -121,16 +137,39 @@ def set_data(self, points, values):
 
 
 class MapMovie2D(Data):
+    r"""
+    2D spatial field :math:`f(x, z, t)` rendered as a heatmap (pcolormesh) animation.
+    """
+
     def __init__(
         self,
         key,
         symbol,
-        plot_coords,
+        plot_coords=[0, 0],
         cmap=DEFAULT_CMAP,
         norm="linear",
         streamlines=None,
         **kwargs,
     ):
+        r"""
+        Initializes a 2D movie field.
+
+         (Refer to :class:`Data` for base parameters)
+        :param cmap: Matplotlib colormap name, defaults to DEFAULT_CMAP.
+        :type cmap: str, optional
+        :param norm: Colorbar scaling. Options usually include 'linear', 'log', or 'TwoSlopeNorm'.
+                     Defaults to "linear".
+        :type norm: str, optional
+        :param streamlines: A list of two Idefix field keys used to show vector streamlines,
+                            e.g., ``["VX1", "VX2"]``. Defaults to None.
+        :type streamlines: list[str], optional
+        :param \**kwargs: Additional rendering options.
+            :keyword streamline_color (str): Color of streamline arrows. Defaults to "w".
+            :keyword compute (callable): Custom function to calculate new fields on the fly.
+            :keyword contours (str): Field key used to draw contour lines over the pcolormesh.
+            :keyword contour_color (str): Color of the contour lines. Defaults to "green".
+        """
+
         # streamlines should be a list like ["VX1", "VX2"]
 
         super().__init__(key, symbol, plot_coords, **kwargs)
@@ -149,31 +188,56 @@ def set_cmap(self, cmap):
         self.cmap = cmap
 
     def set_XYgrid(self, X, Y):
+        """
+        Assign the spatial cartesian grid used for rendering the 2D pcolormesh.
+
+        :param X: 2D array of horizontal coordinates.
+        :type X: numpy.ndarray
+        :param Y: 2D array of vertical coordinates.
+        :type Y: numpy.ndarray
+        """
         self.X, self.Y = X, Y
 
 
 class Field1D(Data):
+    """
+    Base class for 1D fields :math:`f(x, t)`.
+    Increments a global counter for indexing in results arrays.
+    """
+
     def __init__(self, *args, **kwargs):
         self.index = next(Data.timeline_instances)
         super().__init__(*args, **kwargs)
 
 
 class LineMovie1D(Field1D):
-    """1D field meant to be rendered as an animated line plot over time."""
+    """
+    For :math:`f(x, t)` fields, renders as a line plot :math:`f(x, t)` that updates every frame.
+    """
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
 
 class SpaceTimeHeatmap(Field1D):
-    """1D field meant to be rendered as a 2D heatmap (x-t)."""
+    """
+    For :math:`f(x, t)` fields, renders a space-time heatmap.
+
+    :keyword cmap: Colormap for the heatmap.
+    :keyword trace_over: List of :class:`PartQuantity` objects to overlay as trajectories.
+    """
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.cmap = kwargs.get("cmap", DEFAULT_CMAP)
+        self.trace_over = kwargs.get("trace_over", [])
 
 
 class PartQuantity(Data):
+    """
+    Tracks Lagrangian particle properties over time.
+    """
+
     partQuantities_instances = count(1)
 
     def __init__(self, *args, **kwargs):
@@ -341,7 +405,23 @@ def __init__(
         streamLines=None,
     ):
         """
-        If there are n dumps, and end=0.5, only 0.5n of the dumps will be read.
+        Coordinates the detection, processing, and rendering of the simulation data.
+        :param Context: The RunContext object containing simulation metadata.
+        :type Context: RunContext
+        :param spaceTimeHeatmaps: Objects defining (x, t) heatmap plots.
+        :type spaceTimeHeatmaps: list[SpaceTimeHeatmap], optional
+        :param movies1D: Objects defining 1D line plot animations.
+        :type movies1D: list[LineMovie1D], optional
+        :param movies2D: Objects defining 2D heatmaps animations.
+        :type movies2D: list[MapMovie2D], optional
+        :param partQuantities: Objects defining particles quantities.
+        :type partQuantities: list[PartQuantity], optional
+        :param zoom: Zoom level for the rendering view (for 2D only currently).
+        :type zoom: float, optional
+        :param end: Fraction of the simulation dumps to process (0 to 1) (deprecated).
+        :type end: float, optional
+        :param streamLines: Configuration for streamlines overlays.
+        :type streamLines: StreamlineConfig, optional
         """
         self.context = Context
         self.userArgs = self.context.args
@@ -364,6 +444,14 @@ def _to_dict(obj_input):
         self.spaceTimeHeatmaps = _to_dict(spaceTimeHeatmaps)
         self.movies1D = _to_dict(movies1D)
         self.movies2D = _to_dict(movies2D)
+
+        original_part_quantity_keys = set(_to_dict(partQuantities).keys())
+        for heatmap in self.spaceTimeHeatmaps.values():
+            for traceover in heatmap.trace_over:
+                if isinstance(traceover, PartQuantity):
+                    if traceover.key not in original_part_quantity_keys:
+                        traceover.is_trace_over = True
+                        partQuantities.append(traceover)
         self.partQuantities = _to_dict(partQuantities)
 
         combined_1D = {**self.movies1D, **self.spaceTimeHeatmaps}
@@ -373,6 +461,9 @@ def _to_dict(obj_input):
         self._apply_config()
 
     def run(self):
+        """
+        Pray.
+        """
         partInfo = self.context.outputTypes_info["particles"]
 
         if self.userArgs.onlyMovie:
@@ -469,7 +560,7 @@ def run(self):
             if self.userArgs.doOnlyFrames:
                 render_list = [render_list[i] for i in self.userArgs.doOnlyFrames]
             with Pool(self.userArgs.jobs) as pool:
-                pool.map(self.process_and_render_frame, render_list)
+                pool.map(self._process_and_render_frame, render_list)
 
             if self.doMovie:
                 tools.movie(
@@ -477,7 +568,7 @@ def run(self):
                     movie_path=self.framesPaths.slice1_video_path,
                 )
 
-    def process_and_render_frame(self, vtkPath):
+    def _process_and_render_frame(self, vtkPath):
         """1. Read VTK Data, 2. Process Physics Math, 3. Render Requested Frame"""
         V = readVTK(vtkPath)
         self.processor.process(V)
@@ -518,7 +609,7 @@ def _apply_config(self):
 
         if len(fields_tobound) > 0:
             bound_list = self.slice1_list if len(self.slice1_list) > 0 else self.vtkList
-            all_bounds = self.get_bounds(
+            all_bounds = self._get_bounds(
                 bound_list[min(len(bound_list), 5) :],
                 fields_tobound,
             )
@@ -546,7 +637,7 @@ def _apply_config(self):
         for qt in self.movies2D:
             LOG(qt, self.movies2D[qt].bounds)
 
-    def get_bounds(self, vtkList, fields):
+    def _get_bounds(self, vtkList, fields):
         """
         Get the bounds (min, max) of all given fields. I recommend not passing the entire vtkList but rather vtkList[1:] to discard the first output(s ?).
 
@@ -562,7 +653,7 @@ def get_bounds(self, vtkList, fields):
 
         with Pool(self.userArgs.jobs) as pool:
             all_bounds = pool.map(
-                self.get_bounds_indiv,
+                self._get_bounds_indiv,
                 [[vtk, mapfields_indexes] for vtk in vtkList],
             )
         all_bounds = np.array(all_bounds)
@@ -577,7 +668,7 @@ def get_bounds(self, vtkList, fields):
             )
         return bounds
 
-    def get_bounds_indiv(self, args):
+    def _get_bounds_indiv(self, args):
         """
         args (list[2]) must have two components:
             first:   vtk_path (str)
diff --git a/idefix2python/renderer.py b/idefix2python/renderer.py
index 4746e06..e26deec 100644
--- a/idefix2python/renderer.py
+++ b/idefix2python/renderer.py
@@ -249,6 +249,7 @@ def render_1D(self, V, vtkPath):
 
             ax.plot(points, V.data[key])
 
+            # To revove?
             if len(field1D.pointsRef) > 0:
                 ax.plot(
                     field1D.points,
@@ -291,14 +292,36 @@ def render_SpaceTimeHeatmap(self):
                 Points,
                 np.transpose(field1D.values),
                 shading="nearest",
-                # TODO More flexible norm
                 cmap="inferno",
+                # TODO More flexible norm and cmap
             )
             cbar = fig.colorbar(cmesh, ax=ax)
 
+            has_legend_items = False
             if len(field1D.pointsRef) > 0:
-                ax.plot(field1D.pointsRef, field1D.valuesRef, label="Predicted")
+                plot_kwargs = {}
+                if hasattr(field1D.ref_function, "plot_kwargs"):
+                    plot_kwargs = field1D.ref_function.plot_kwargs
+                if "label" not in plot_kwargs:
+                    plot_kwargs["label"] = "Predicted"
+                ax.plot(
+                    field1D.pointsRef,
+                    field1D.valuesRef,
+                    **plot_kwargs,
+                )
+                has_legend_items = True
+
+            for trace_over in field1D.trace_over:
+                ax.plot(
+                    trace_over.points,
+                    trace_over.values,
+                    label=trace_over.symbol,
+                    color="lime",  # TODO customize this in later PR
+                )
+                has_legend_items = True
+            if has_legend_items:
                 ax.legend()
+
             ax.set_ylim(np.min(field1D.points), np.max(field1D.points))
             cbar.ax.set_title(field1D.symbol)
             ax.set_xlabel(r"$t$", fontsize=LABEL_FONTSIZE)
@@ -315,6 +338,8 @@ def render_timeSeries(self):
         fig, axs = self._setup_figure(self.partQuantities)
 
         for key, qty in self.partQuantities.items():
+            if getattr(qty, "is_trace_over", False):
+                continue
             ax = axs[*qty.plot_coords]
             T = np.asarray(self.context.outputTypes_info["particles"].times)
             ax.plot(T, qty.values, lw=2)