diff --git a/.gitignore b/.gitignore
index 53014d8..09919c8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -339,3 +339,4 @@ src/sdf_xarray/_version.py
 
 # Downloaded doc tutorial datasets
 docs/tutorial_*
+docs/epoch_workshop_2026/datasets
diff --git a/docs/conf.py b/docs/conf.py
index a7e4d87..3ac20d5 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -153,3 +153,23 @@
         continue
     else:
         fetch_dataset(dataset, save_path=cwd)
+
+epoch_workshop_datasets = [
+    "1_1_drifting_bunch",
+    "2_1_two_stream_instability",
+    "3_3_Gaussian_1d_laser",
+    "3_5_Gaussian_beam",
+    "4_2_self_heating",
+    "4_3_basic_target",
+    "4_4_momentum_distribution",
+    "5_1_probe",
+    "5_2_subsets",
+]
+
+epoch_workshop_dataset_folder = cwd / "epoch_workshop_2026" / "datasets"
+for dataset in epoch_workshop_datasets:
+    # If the dataset already exists then don't download it again
+    if (epoch_workshop_dataset_folder / dataset).exists():
+        continue
+    else:
+        fetch_dataset(dataset, save_path=epoch_workshop_dataset_folder)
diff --git a/docs/epoch_workshop_2026/animating/animate_1D_density.py b/docs/epoch_workshop_2026/animating/animate_1D_density.py
new file mode 100644
index 0000000..073aa3e
--- /dev/null
+++ b/docs/epoch_workshop_2026/animating/animate_1D_density.py
@@ -0,0 +1,19 @@
+from pathlib import Path
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/1_1_drifting_bunch")
+ds = sdfxr.open_mfdataset(input_dir)
+
+# Convert the time to femtoseconds
+ds = ds.epoch.rescale_coords(1e15, "fs", "time")
+# Convert the x and y coords to microns
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid"])
+
+anim = ds["Derived_Number_Density"].epoch.animate()
+
+# Visualise it in a Jupyter notebook
+anim.show()
+
+# Or save the animation
+# anim.save(input_dir / "number_density.gif", fps=5)
diff --git a/docs/epoch_workshop_2026/animating/animate_1D_poynting_flux.py b/docs/epoch_workshop_2026/animating/animate_1D_poynting_flux.py
new file mode 100644
index 0000000..dc8bd0d
--- /dev/null
+++ b/docs/epoch_workshop_2026/animating/animate_1D_poynting_flux.py
@@ -0,0 +1,33 @@
+from pathlib import Path
+
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/3_3_Gaussian_1d_laser")
+ds = sdfxr.open_mfdataset(input_dir)
+
+# Convert the time to femtoseconds
+ds = ds.epoch.rescale_coords(1e15, "fs", "time")
+# Convert the x and y coords to microns
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid"])
+
+# Calculate Poynting flux magnitude
+flux_magnitude = np.sqrt(
+    ds["Derived_Poynting_Flux_x"] ** 2
+    + ds["Derived_Poynting_Flux_y"] ** 2
+    + ds["Derived_Poynting_Flux_z"] ** 2
+)
+
+# convert to W/cm^2
+I_Wcm2 = flux_magnitude * 1e-4
+I_Wcm2.attrs["long_name"] = "Poynting Flux Magnitude"
+I_Wcm2.attrs["units"] = "W/cm$^2$"
+
+anim = I_Wcm2.epoch.animate()
+
+# Visualise it in a Jupyter notebook
+anim.show()
+
+# Or save the animation
+# anim.save(input_dir / "laser.gif", fps=10)
diff --git a/docs/epoch_workshop_2026/animating/animate_2D_dist_fn_multispecies.py b/docs/epoch_workshop_2026/animating/animate_2D_dist_fn_multispecies.py
new file mode 100644
index 0000000..a5bb913
--- /dev/null
+++ b/docs/epoch_workshop_2026/animating/animate_2D_dist_fn_multispecies.py
@@ -0,0 +1,25 @@
+from pathlib import Path
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/2_1_two_stream_instability")
+ds = sdfxr.open_mfdataset(
+    input_dir, data_vars=["dist_fn_x_px_Left", "dist_fn_x_px_Right"]
+)
+
+# Rescale coords to account for kilometers
+ds = ds.epoch.rescale_coords(1e-3, "km", ["X_x_px_Left"])
+
+# Sum phase-space of species "Left" and "Right" in "x_px" distribution function
+# NOTE: We only use the values from the right distribution function as if we inherit
+# the coords we from the right we end up with 4 coords instead of 2
+total_phase_space = ds["dist_fn_x_px_Left"] + ds["dist_fn_x_px_Right"].values
+total_phase_space.attrs["long_name"] = "Phase Space Distribution"
+total_phase_space.attrs["units"] = "kg.m/s"
+
+anim = total_phase_space.epoch.animate()
+# Visualise it in a Jupyter notebook
+anim.show()
+
+# Or save the animation
+# anim.save(input_dir / "phase_space.gif")
diff --git a/docs/epoch_workshop_2026/animating/animate_2D_dist_fn_multispecies_alternative.py b/docs/epoch_workshop_2026/animating/animate_2D_dist_fn_multispecies_alternative.py
new file mode 100644
index 0000000..e8fb86a
--- /dev/null
+++ b/docs/epoch_workshop_2026/animating/animate_2D_dist_fn_multispecies_alternative.py
@@ -0,0 +1,56 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/2_1_two_stream_instability")
+ds = sdfxr.open_mfdataset(
+    input_dir, data_vars=["dist_fn_x_px_Left", "dist_fn_x_px_Right"]
+)
+
+# Rescale coords to account for kilometers
+ds = ds.epoch.rescale_coords(1e-3, "km", ["X_x_px_Left"])
+
+# Sum phase-space of species "Left" and "Right" in "x_px" distribution function
+# NOTE: We only use the values from the right distribution function as if we inherit
+# the coords we from the right we end up with 4 coords instead of 2
+total_phase_space = ds["dist_fn_x_px_Left"] + ds["dist_fn_x_px_Right"].values
+total_phase_space.attrs["long_name"] = "Phase Space Distribution"
+total_phase_space.attrs["units"] = "kg.m/s"
+
+fig, ax = plt.subplots()
+
+# construct first frame
+# Data needs to be transposed so that axes match
+total_phase_space = total_phase_space.T
+col_min = total_phase_space.values.min()
+col_max = total_phase_space.values.max()
+gif_frame = ax.pcolormesh(
+    total_phase_space["X_x_px_Left"],
+    total_phase_space["Px_x_px_Left"],
+    total_phase_space.isel(time=0),
+    vmin=col_min,
+    vmax=col_max,
+)
+ax.set_xlabel(
+    f"{total_phase_space['X_x_px_Left'].attrs['long_name']} [{total_phase_space['X_x_px_Left'].attrs['units']}]"
+)
+ax.set_ylabel(
+    f"{total_phase_space.attrs['long_name']} [{total_phase_space.attrs['units']}]"
+)
+ax.set_title(f"t = {ds['time'].values[0]:.3f} [{ds['time'].attrs['units']}]")
+cbar = fig.colorbar(gif_frame, ax=ax)
+cbar.set_label("Summed particle weight per bin")
+
+
+# Make GIF
+def update(i):
+    gif_frame.set_array(total_phase_space.isel(time=i))
+    ax.set_title(f"t = {ds['time'].values[i]:.3f} [{ds['time'].attrs['units']}]")
+    return (gif_frame,)
+
+
+anim = FuncAnimation(fig, update, frames=ds.sizes["time"])
+anim.save(input_dir / "phase_space.gif")
diff --git a/docs/epoch_workshop_2026/animating/animate_2D_poynting_flux.py b/docs/epoch_workshop_2026/animating/animate_2D_poynting_flux.py
new file mode 100644
index 0000000..0f09b09
--- /dev/null
+++ b/docs/epoch_workshop_2026/animating/animate_2D_poynting_flux.py
@@ -0,0 +1,33 @@
+from pathlib import Path
+
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/3_5_Gaussian_beam")
+ds = sdfxr.open_mfdataset(input_dir)
+
+# Convert the time to femtoseconds
+ds = ds.epoch.rescale_coords(1e15, "fs", "time")
+# Convert the x and y coords to microns
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid", "Y_Grid_mid"])
+
+# Calculate Poynting flux magnitude
+flux_magnitude = np.sqrt(
+    ds["Derived_Poynting_Flux_x"] ** 2
+    + ds["Derived_Poynting_Flux_y"] ** 2
+    + ds["Derived_Poynting_Flux_z"] ** 2
+)
+
+# convert to W/cm^2
+I_Wcm2 = flux_magnitude * 1e-4
+I_Wcm2.attrs["long_name"] = "Poynting Flux Magnitude"
+I_Wcm2.attrs["units"] = "W/cm$^2$"
+
+anim = I_Wcm2.epoch.animate()
+
+# Visualise it in a Jupyter notebook
+anim.show()
+
+# Or save the animation
+# anim.save(input_dir / "laser.gif", fps=10)
diff --git a/docs/epoch_workshop_2026/animations.md b/docs/epoch_workshop_2026/animations.md
new file mode 100644
index 0000000..1876dc1
--- /dev/null
+++ b/docs/epoch_workshop_2026/animations.md
@@ -0,0 +1,48 @@
+---
+file_format: mystnb
+kernelspec:
+  name: python3
+---
+
+# Animations
+
+```{note}
+Unfortunately, these animations are not interactive like the rest of the documentation as building them on readthedocs takes too long so we display the GIFs from the datasets instead.
+```
+
+## Animate Density 1D
+
+- **input deck:** <path:datasets/1_1_drifting_bunch/input.deck>
+- **Python File:** <path:animating/animate_1D_density.py>
+
+```{literalinclude} ./animating/animate_1D_density.py
+```
+![datasets/1_1_drifting_bunch/number_density.gif](datasets/1_1_drifting_bunch/number_density.gif)
+
+## Animate Poynting Flux 1D
+
+- **input deck:** <path:datasets/3_3_Gaussian_1d_laser/input.deck>
+- **Python File:** <path:animating/animate_1D_poynting_flux.py>
+
+```{literalinclude} ./animating/animate_1D_poynting_flux.py
+```
+![datasets/3_3_Gaussian_1d_laser/laser.gif](datasets/3_3_Gaussian_1d_laser/laser.gif)
+
+## Animate Distribution Function Multi-Species 2D
+
+- **input deck:** <path:datasets/2_1_two_stream_instability/input.deck>
+- **Python File:** <path:animating/animate_2D_dist_fn_multispecies.py>
+- **Python File (alternative):** <path:animating/animate_2D_dist_fn_multispecies_alternative.py>
+
+```{literalinclude} ./animating/animate_2D_dist_fn_multispecies.py
+```
+![datasets/2_1_two_stream_instability/phase_space.gif](datasets/2_1_two_stream_instability/phase_space.gif)
+
+## Animate Poynting Flux 2D
+
+- **input deck:** <path:datasets/3_5_Gaussian_beam/input.deck>
+- **Python File:** <path:animating/animate_2D_poynting_flux.py>
+
+```{literalinclude} ./animating/animate_2D_poynting_flux.py
+```
+![datasets/3_5_Gaussian_beam/laser.gif](datasets/3_5_Gaussian_beam/laser.gif)
diff --git a/docs/epoch_workshop_2026/histograms.md b/docs/epoch_workshop_2026/histograms.md
new file mode 100644
index 0000000..cd6f913
--- /dev/null
+++ b/docs/epoch_workshop_2026/histograms.md
@@ -0,0 +1,34 @@
+---
+file_format: mystnb
+kernelspec:
+  name: python3
+---
+
+# Histograms
+
+## Plot Histogram Kinetic Energy Probes
+
+- **input deck:** <path:datasets/5_1_probe/input.deck>
+- **Python File:** <path:histograms/plot_hist_ke_probe.py>
+
+```{code-cell} ipython3
+:load: ./histograms/plot_hist_ke_probe.py
+```
+
+## Plot Histogram Kinetic Energy
+
+- **input deck:** <path:datasets/4_3_basic_target/input.deck>
+- **Python File:** <path:histograms/plot_hist_ke.py>
+
+```{code-cell} ipython3
+:load: ./histograms/plot_hist_ke.py
+```
+
+## Plot Histogram Electron Momentum
+
+- **input deck:** <path:datasets/4_4_momentum_distribution/input.deck>
+- **Python File:** <path:histograms/plot_hist_px.py>
+
+```{code-cell} ipython3
+:load: ./histograms/plot_hist_px.py
+```
\ No newline at end of file
diff --git a/docs/epoch_workshop_2026/histograms/plot_hist_ke.py b/docs/epoch_workshop_2026/histograms/plot_hist_ke.py
new file mode 100644
index 0000000..3275d82
--- /dev/null
+++ b/docs/epoch_workshop_2026/histograms/plot_hist_ke.py
@@ -0,0 +1,38 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+# Constants
+mass = 9.1093837e-31
+c = 299792458
+q0 = 1.60217663e-19
+
+input_dir = Path("datasets/4_3_basic_target")
+ds = sdfxr.open_dataset(input_dir / "0000.sdf", keep_particles=True)
+particles_magnitude = (
+    ds["Particles_Px_Electron"] ** 2
+    + ds["Particles_Py_Electron"] ** 2
+    + ds["Particles_Pz_Electron"] ** 2
+)
+ke_MeV = (np.sqrt(particles_magnitude * c**2 + mass**2 * c**4) - mass * c**2) / (
+    1.0e6 * q0
+)
+
+bin_edges = np.linspace(ke_MeV.min().values, ke_MeV.max().values, 101)
+bin_N, _ = np.histogram(
+    ke_MeV.values, bins=bin_edges, weights=ds["Particles_Weight_Electron"]
+)
+dKE = np.diff(bin_edges)
+dN_dKE = bin_N / dKE
+bin_centres = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+
+plt.plot(bin_centres, dN_dKE)
+plt.xlabel("Kinetic energy [MeV]")
+plt.ylabel("dN/dKE [1/MeV]")
+
+plt.savefig(input_dir / "ke_spectrum.png", dpi=300)
+np.savetxt(input_dir / "ke_spectrum_vals.txt", ke_MeV.values)
+np.savetxt(input_dir / "ke_spectrum_ke.txt", bin_centres)
diff --git a/docs/epoch_workshop_2026/histograms/plot_hist_ke_probe.py b/docs/epoch_workshop_2026/histograms/plot_hist_ke_probe.py
new file mode 100644
index 0000000..ea046c8
--- /dev/null
+++ b/docs/epoch_workshop_2026/histograms/plot_hist_ke_probe.py
@@ -0,0 +1,47 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+# Constants
+mass = 9.1093837e-31
+c = 299792458
+q0 = 1.60217663e-19
+
+input_dir = Path("datasets/5_1_probe")
+# Since we're loading particle and probes we need to specify their name
+# from the input.deck
+ds = sdfxr.open_mfdataset(
+    input_dir, keep_particles=True, probe_names=["Electron_probe"]
+)
+
+px = ds["Electron_probe_Px"].values.flatten()
+py = ds["Electron_probe_Py"].values.flatten()
+pz = ds["Electron_probe_Pz"].values.flatten()
+probe_weights_raw = ds["Electron_probe_weight"].values.flatten()
+
+# Create a mask to remove NaNs
+mask = ~np.isnan(probe_weights_raw)
+probe_weights = probe_weights_raw[mask]
+probe_px = px[mask]
+probe_py = py[mask]
+probe_pz = pz[mask]
+
+probe_magnitude = probe_px**2 + probe_py**2 + probe_pz**2
+ke_MeV = (np.sqrt(probe_magnitude * c**2 + mass**2 * c**4) - mass * c**2) / (1.0e6 * q0)
+
+bin_edges = np.linspace(ke_MeV.min(), ke_MeV.max(), 101)
+bin_N, _ = np.histogram(ke_MeV, bins=bin_edges, weights=probe_weights)
+dKE = np.diff(bin_edges)
+dN_dKE = bin_N / dKE
+bin_centres = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+
+plt.plot(bin_centres, dN_dKE)
+plt.xlabel("Kinetic energy [MeV]")
+plt.ylabel("dN/dKE [1/MeV]")
+
+plt.savefig(input_dir / "ke_spectrum.png", dpi=300)
+np.savetxt(input_dir / "ke_spectrum_vals.txt", ke_MeV)
+np.savetxt(input_dir / "ke_spectrum_ke.txt", bin_centres)
diff --git a/docs/epoch_workshop_2026/histograms/plot_hist_px.py b/docs/epoch_workshop_2026/histograms/plot_hist_px.py
new file mode 100644
index 0000000..715ec61
--- /dev/null
+++ b/docs/epoch_workshop_2026/histograms/plot_hist_px.py
@@ -0,0 +1,26 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/4_4_momentum_distribution")
+ds = sdfxr.open_dataset(input_dir / "0000.sdf", keep_particles=True)
+px = ds["Particles_Px_Electron_user"].values
+
+bin_edges = np.linspace(px.min(), px.max(), 101)
+bin_N, _ = np.histogram(
+    px, bins=bin_edges, weights=ds["Particles_Weight_Electron_user"]
+)
+dpx = np.diff(bin_edges)
+dN_dpx = bin_N / dpx
+bin_centres = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+
+plt.plot(bin_centres, dN_dpx)
+plt.xlabel("Px [kg.m/s]")
+plt.ylabel("dN/dp$_x$ [s/(kg.m)]")
+
+plt.savefig(input_dir / "px_spectrum.png", dpi=300)
+np.savetxt(input_dir / "px_spectrum_vals.txt", dN_dpx)
+np.savetxt(input_dir / "ke_spectrum_px.txt", bin_centres)
diff --git a/docs/epoch_workshop_2026/plots.md b/docs/epoch_workshop_2026/plots.md
new file mode 100644
index 0000000..cb4be12
--- /dev/null
+++ b/docs/epoch_workshop_2026/plots.md
@@ -0,0 +1,70 @@
+---
+file_format: mystnb
+kernelspec:
+  name: python3
+---
+
+# Plots
+
+## Plot 1D Density
+
+- **input deck:** <path:datasets/1_1_drifting_bunch/input.deck>
+- **Python File:** <path:plotting/plot_1D_density.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_1D_density.py
+```
+
+## Plot 1D Poynting Flux
+
+- **input deck:** <path:datasets/3_3_Gaussian_1d_laser/input.deck>
+- **Python File:** <path:plotting/plot_1D_poynting_flux.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_1D_poynting_flux.py
+```
+
+## Plot 2D Density
+
+- **input deck:** <path:datasets/4_3_basic_target/input.deck>
+- **Python File:** <path:plotting/plot_2D_density.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_2D_density.py
+```
+
+## Plot 2D Distribution Function Multispecies
+
+- **input deck:** <path:datasets/2_1_two_stream_instability/input.deck>
+- **Python File:** <path:plotting/plot_2D_dist_fn_multispecies.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_2D_dist_fn_multispecies.py
+```
+
+## Plot 2D Poynting Flux
+
+- **input deck:** <path:datasets/3_5_Gaussian_beam/input.deck>
+- **Python File:** <path:plotting/plot_2D_poynting_flux.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_2D_poynting_flux.py
+```
+
+## Plot Time vs Temperature
+
+- **input deck:** <path:datasets/4_2_self_heating/input.deck>
+- **Python File:** <path:plotting/plot_time_temperature.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_time_temperature.py
+```
+
+## Plot X vs Px Scatter
+
+- **input deck:** <path:datasets/5_2_subsets/input.deck>
+- **Python File:** <path:plotting/plot_x_px_scatter.py>
+
+```{code-cell} ipython3
+:load: ./plotting/plot_x_px_scatter.py
+```
diff --git a/docs/epoch_workshop_2026/plotting/plot_1D_density.py b/docs/epoch_workshop_2026/plotting/plot_1D_density.py
new file mode 100644
index 0000000..7220d40
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_1D_density.py
@@ -0,0 +1,18 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/1_1_drifting_bunch")
+
+ds = sdfxr.open_dataset(input_dir / "0000.sdf")
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid"])
+
+ds["Derived_Number_Density"].epoch.plot()
+plt.tight_layout()
+
+plt.savefig(input_dir / "number_density.png", dpi=300)
+np.savetxt(input_dir / "number_density_vals.txt", ds["Derived_Number_Density"].values)
+np.savetxt(input_dir / "number_density_x.txt", ds["X_Grid_mid"].values)
diff --git a/docs/epoch_workshop_2026/plotting/plot_1D_poynting_flux.py b/docs/epoch_workshop_2026/plotting/plot_1D_poynting_flux.py
new file mode 100644
index 0000000..71e99f6
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_1D_poynting_flux.py
@@ -0,0 +1,32 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/3_3_Gaussian_1d_laser")
+
+ds = sdfxr.open_dataset(input_dir / "0020.sdf")
+
+# Convert the x and y coords to microns
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid"])
+
+# Calculate Poynting flux magnitude
+flux_magnitude = np.sqrt(
+    ds["Derived_Poynting_Flux_x"] ** 2
+    + ds["Derived_Poynting_Flux_y"] ** 2
+    + ds["Derived_Poynting_Flux_z"] ** 2
+)
+
+# convert to W/cm^2
+I_Wcm2 = flux_magnitude * 1e-4
+I_Wcm2.attrs["long_name"] = "Poynting Flux Magnitude"
+I_Wcm2.attrs["units"] = "W/cm$^2$"
+
+I_Wcm2.epoch.plot()
+plt.tight_layout()
+
+plt.savefig(input_dir / "intensity.png", dpi=300)
+np.savetxt(input_dir / "intensity_vals.txt", I_Wcm2.values)
+np.savetxt(input_dir / "intensity_x.txt", ds["X_Grid_mid"].values)
diff --git a/docs/epoch_workshop_2026/plotting/plot_2D_density.py b/docs/epoch_workshop_2026/plotting/plot_2D_density.py
new file mode 100644
index 0000000..811924f
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_2D_density.py
@@ -0,0 +1,18 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/4_3_basic_target")
+
+ds = sdfxr.open_dataset(input_dir / "0000.sdf")
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid", "Y_Grid_mid"])
+
+ds["Derived_Number_Density"].epoch.plot()
+plt.tight_layout()
+
+plt.savefig(input_dir / "number_density.png", dpi=300)
+np.savetxt(input_dir / "number_density_vals.txt", ds["Derived_Number_Density"].values)
+np.savetxt(input_dir / "number_density_x.txt", ds["X_Grid_mid"].values)
diff --git a/docs/epoch_workshop_2026/plotting/plot_2D_dist_fn_multispecies.py b/docs/epoch_workshop_2026/plotting/plot_2D_dist_fn_multispecies.py
new file mode 100644
index 0000000..a56cea8
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_2D_dist_fn_multispecies.py
@@ -0,0 +1,26 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/2_1_two_stream_instability")
+ds = sdfxr.open_dataset(input_dir / "0000.sdf")
+
+ds = ds.epoch.rescale_coords(1e-3, "km", ["X_x_px_Left"])
+
+# Sum phase-space of species "Left" and "Right" in "x_px" distribution function
+# NOTE: We only use the values from the right distribution function as if we inherit
+# the coords we from the right we end up with 4 coords instead of 2
+total_phase_space = ds["dist_fn_x_px_Left"] + ds["dist_fn_x_px_Right"].values
+total_phase_space.attrs["long_name"] = "Summed Particle Weight Per Bin"
+total_phase_space.attrs["units"] = ""
+
+total_phase_space.epoch.plot()
+plt.tight_layout()
+
+plt.savefig(input_dir / "phase_space.png", dpi=300)
+np.savetxt(input_dir / "phase_space_vals.txt", total_phase_space.values)
+np.savetxt(input_dir / "phase_space_x.txt", ds["X_x_px_Left"].values)
+np.savetxt(input_dir / "phase_space_y.txt", ds["Px_x_px_Left"].values)
diff --git a/docs/epoch_workshop_2026/plotting/plot_2D_poynting_flux.py b/docs/epoch_workshop_2026/plotting/plot_2D_poynting_flux.py
new file mode 100644
index 0000000..99a6bb3
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_2D_poynting_flux.py
@@ -0,0 +1,33 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/3_5_Gaussian_beam")
+
+ds = sdfxr.open_dataset(input_dir / "0001.sdf")
+
+# Convert the x and y coords to microns
+ds = ds.epoch.rescale_coords(1e6, "µm", ["X_Grid_mid", "Y_Grid_mid"])
+
+# Calculate Poynting flux magnitude
+flux_magnitude = np.sqrt(
+    ds["Derived_Poynting_Flux_x"] ** 2
+    + ds["Derived_Poynting_Flux_y"] ** 2
+    + ds["Derived_Poynting_Flux_z"] ** 2
+)
+
+# convert to W/cm^2
+I_Wcm2 = flux_magnitude * 1e-4
+I_Wcm2.attrs["long_name"] = "Poynting Flux Magnitude"
+I_Wcm2.attrs["units"] = "W/cm$^2$"
+
+I_Wcm2.epoch.plot()
+plt.tight_layout()
+
+plt.savefig(input_dir / "intensity.png", dpi=300)
+np.savetxt(input_dir / "intensity_vals.txt", I_Wcm2.values)
+np.savetxt(input_dir / "intensity_x.txt", ds["X_Grid_mid"].values)
+np.savetxt(input_dir / "intensity_y.txt", ds["Y_Grid_mid"].values)
diff --git a/docs/epoch_workshop_2026/plotting/plot_time_temperature.py b/docs/epoch_workshop_2026/plotting/plot_time_temperature.py
new file mode 100644
index 0000000..cdab64c
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_time_temperature.py
@@ -0,0 +1,29 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+# Constants
+kB = 1.380649e-23
+q0 = 1.60217663e-19
+
+input_dir = Path("datasets/4_2_self_heating")
+ds = sdfxr.open_mfdataset(input_dir)
+
+# Convert the time to femtoseconds
+ds = ds.epoch.rescale_coords(1e15, "fs", "time")
+
+# Averate temperature over all spatial cells at each time-step
+kB = 1.380649e-23
+q0 = 1.60217663e-19
+
+temp_ev = ds["Derived_Temperature"].mean(dim=("X_Grid_mid", "Y_Grid_mid")) * kB / q0
+
+temp_ev.epoch.plot()
+plt.tight_layout()
+
+plt.savefig(input_dir / "temperature_time.png", dpi=300)
+np.savetxt(input_dir / "temperature_time_time.txt", temp_ev["time"].values)
+np.savetxt(input_dir / "temperature_time_temperature.txt", temp_ev.values)
diff --git a/docs/epoch_workshop_2026/plotting/plot_x_px_scatter.py b/docs/epoch_workshop_2026/plotting/plot_x_px_scatter.py
new file mode 100644
index 0000000..19df9ba
--- /dev/null
+++ b/docs/epoch_workshop_2026/plotting/plot_x_px_scatter.py
@@ -0,0 +1,24 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import sdf_xarray as sdfxr
+
+input_dir = Path("datasets/5_2_subsets")
+ds = sdfxr.open_dataset(input_dir / "0000.sdf", keep_particles=True)
+
+ds = ds.epoch.rescale_coords(
+    1e6, "µm", ["X_Particles_subset_Refluxers_Electron", "Y_Target_mid"]
+)
+
+x = ds["X_Particles_subset_Refluxers_Electron"]
+px = ds["Particles_Px_subset_Refluxers_Electron"]
+
+plt.scatter(x, px)
+plt.xlabel("x [µm]")
+plt.ylabel("Px [kg.m/s]")
+
+plt.savefig(input_dir / "x_px.png", dpi=300)
+np.savetxt(input_dir / "x_px_vals.txt", px.values)
+np.savetxt(input_dir / "x_px_x.txt", x.values)
diff --git a/docs/index.md b/docs/index.md
index 89903c5..c824c38 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -24,6 +24,15 @@ unit_conversion.md
 animation.md
 ```
 
+```{toctree}
+:caption: 'EPOCH Workshop 2026'
+
+epoch_workshop_2026/plots.md
+epoch_workshop_2026/animations.md
+epoch_workshop_2026/histograms.md
+```
+
+
 ```{toctree}
 :caption: 'Development'
 :maxdepth: 1
diff --git a/src/sdf_xarray/download.py b/src/sdf_xarray/download.py
index 0a67b7d..4faef35 100644
--- a/src/sdf_xarray/download.py
+++ b/src/sdf_xarray/download.py
@@ -1,9 +1,6 @@
 from pathlib import Path
 from shutil import move
-from typing import TYPE_CHECKING, Literal, TypeAlias
-
-if TYPE_CHECKING:
-    import pooch  # noqa: F401
+from typing import Literal, TypeAlias
 
 DatasetName: TypeAlias = Literal[
     "test_array_no_grids",
@@ -17,6 +14,15 @@
     "tutorial_dataset_2d",
     "tutorial_dataset_2d_moving_window",
     "tutorial_dataset_3d",
+    "1_1_drifting_bunch",
+    "2_1_two_stream_instability",
+    "3_3_Gaussian_1d_laser",
+    "3_5_Gaussian_beam",
+    "4_2_self_heating",
+    "4_3_basic_target",
+    "4_4_momentum_distribution",
+    "5_1_probe",
+    "5_2_subsets",
 ]
 
 
@@ -54,22 +60,11 @@ def fetch_dataset(
     logger = pooch.get_logger()
     datasets = pooch.create(
         path=pooch.os_cache("sdf_datasets"),
-        base_url="https://zenodo.org/records/17991042/files",
-        registry={
-            "test_array_no_grids.zip": "md5:583c85ed8c31d0e34e7766b6d9f2d6da",
-            "test_dist_fn.zip": "md5:a582ff5e8c59bad62fe4897f65fc7a11",
-            "test_files_1D.zip": "md5:42e53b229556c174c538c5481c4d596a",
-            "test_files_2D_moving_window.zip": "md5:3744483bbf416936ad6df8847c54dad1",
-            "test_files_3D.zip": "md5:a679e71281bab1d373dc4980e6da1a7c",
-            "test_mismatched_files.zip": "md5:710fdc94666edf7777523e8fc9dd1bd4",
-            "test_two_probes_2D.zip": "md5:0f2a4fefe84a15292d066b3320d4d533",
-            "tutorial_dataset_1d.zip": "md5:7fad744d8b8b2b84bba5c0e705fdef7b",
-            "tutorial_dataset_2d.zip": "md5:b7f35c05703a48eb5128049cdd106ffa",
-            "tutorial_dataset_2d_moving_window.zip": "md5:a795f40d18df69263842055de4559501",
-            "tutorial_dataset_3d.zip": "md5:d9254648867016292440fdb028f717f7",
-        },
+        base_url="doi:10.5281/zenodo.17618509",
+        registry=None,
         retry_if_failed=10,
     )
+    datasets.load_registry_from_doi()
 
     datasets.fetch(
         f"{dataset_name}.zip", processor=pooch.Unzip(extract_dir="."), progressbar=True