Add EPOCH Workshop 2026 scripts + datasets

.gitignore

@@ -339,3 +339,4 @@ src/sdf_xarray/_version.py
 
 # Downloaded doc tutorial datasets
 docs/tutorial_*
+docs/epoch_workshop_2026/datasets

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)

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)

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)

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")

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")

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)

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)

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
+```

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)

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)

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)