From bbc00915b467a12059708dcf87b7cccd09e717ba Mon Sep 17 00:00:00 2001
From: Max <max.lindqvist@ipp.mpg.de>
Date: Mon, 28 Sep 2026 08:07:07 +0200
Subject: [PATCH 1/4] Use the improvements branch

---
 submodules/struphy-plots | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/submodules/struphy-plots b/submodules/struphy-plots
index 3bdfea5..bf7ba94 160000
--- a/submodules/struphy-plots
+++ b/submodules/struphy-plots
@@ -1 +1 @@
-Subproject commit 3bdfea5caa28dd3720e2246a90331232f2f3dac7
+Subproject commit bf7ba94fdb2332bb0d06e6a04c134b16b1740843

From bd7c2b23879e6da40b3204d6601a59b95adb7e89 Mon Sep 17 00:00:00 2001
From: Max <max.lindqvist@ipp.mpg.de>
Date: Mon, 28 Sep 2026 08:33:32 +0200
Subject: [PATCH 2/4] Update submodules/struphy-plots

---
 submodules/struphy-plots | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/submodules/struphy-plots b/submodules/struphy-plots
index bf7ba94..e6dfc8d 160000
--- a/submodules/struphy-plots
+++ b/submodules/struphy-plots
@@ -1 +1 @@
-Subproject commit bf7ba94fdb2332bb0d06e6a04c134b16b1740843
+Subproject commit e6dfc8d1920c45f2c7b39903eb51d6db6cf4a81c

From df6b1d2ca85b85bad351c8d2e4f236e765c6b22c Mon Sep 17 00:00:00 2001
From: Max <max.lindqvist@ipp.mpg.de>
Date: Mon, 28 Sep 2026 08:33:46 +0200
Subject: [PATCH 3/4] Updated example, use the new save helper

---
 docs/src/examples/README.md                   | 10 ++---
 docs/src/examples/acoustic-pulse.py           | 24 ++--------
 docs/src/examples/alfven-standing-wave.py     | 26 ++---------
 docs/src/examples/beltrami-sph.py             | 32 +++-----------
 docs/src/examples/bump-on-tail.py             | 28 ++----------
 docs/src/examples/coaxial-waveguide.py        | 28 ++----------
 docs/src/examples/cold-plasma-oscillation.py  | 26 ++---------
 docs/src/examples/cold-plasma-wave-packet.py  | 26 ++---------
 docs/src/examples/cold-plasma-waves.py        | 26 ++---------
 docs/src/examples/dam-break.py                | 26 ++---------
 docs/src/examples/damped-alfven-wave.py       | 26 ++---------
 docs/src/examples/diffusion-methods.py        | 22 +---------
 docs/src/examples/diocotron-instability.py    | 28 ++----------
 docs/src/examples/faraday-rotation.py         | 26 ++---------
 docs/src/examples/gas-expansion.py            | 30 +++----------
 docs/src/examples/grad-b-drift.py             | 26 ++---------
 docs/src/examples/guiding-center-orbits.py    | 30 +++----------
 docs/src/examples/gvec-equilibrium.py         | 30 +++----------
 docs/src/examples/gyromotion.py               | 26 ++---------
 docs/src/examples/hall-mhd-waves.py           | 26 ++---------
 docs/src/examples/hasegawa-wakatani.py        | 26 ++---------
 .../examples/hybrid-alfven-ion-coupling.py    | 24 +---------
 docs/src/examples/hybrid-current-coupling.py  | 28 ++----------
 .../incompressible-shear-relaxation.py        | 22 +---------
 docs/src/examples/itg-drift-wave.py           | 36 ++++-----------
 docs/src/examples/langmuir-wave-dispersion.py | 26 ++---------
 .../linear-dissipative-alfven-wave.py         | 26 ++---------
 .../src/examples/maxwell-cavity-resonances.py | 26 ++---------
 docs/src/examples/maxwell-curved-mesh.py      | 26 ++---------
 .../maxwell-structure-preservation.py         | 26 ++---------
 docs/src/examples/maxwell-wave.py             | 26 ++---------
 docs/src/examples/mhd-slab-waves.py           | 24 +---------
 docs/src/examples/ordinary-mode-dispersion.py | 26 ++---------
 docs/src/examples/orszag-tang-vortex.py       | 28 ++----------
 docs/src/examples/poisson-convergence.py      | 26 ++---------
 docs/src/examples/poisson-source.py           | 24 +---------
 docs/src/examples/pressureless-transport.py   | 26 ++---------
 docs/src/examples/resistive-diffusion.py      | 26 ++---------
 docs/src/examples/resistive-x-point.py        | 28 ++----------
 docs/src/examples/shear-alfven-wave.py        | 26 ++---------
 docs/src/examples/sph-velocity-diffusion.py   | 24 +---------
 docs/src/examples/strong-landau-damping.py    | 24 +---------
 docs/src/examples/toroidal-shear-alfven.py    | 44 +++++--------------
 docs/src/examples/two-stream-instability.py   | 28 ++----------
 docs/src/examples/vlasov-tokamak.py           | 26 ++---------
 docs/src/examples/vortex-merger.py            | 26 ++---------
 docs/src/examples/weak-landau-damping.py      | 24 ++--------
 docs/src/examples/weibel-instability.py       | 28 ++----------
 docs/src/examples/zeldovich-caustic.py        | 22 +---------
 49 files changed, 173 insertions(+), 1121 deletions(-)

diff --git a/docs/src/examples/README.md b/docs/src/examples/README.md
index 2652295..db19c4d 100644
--- a/docs/src/examples/README.md
+++ b/docs/src/examples/README.md
@@ -86,10 +86,10 @@ filename.
 - `pproc(sim, show=False)` contains only the physics. Draw the figures with struphy-plots and
   `backend="plotly"`, e.g. `e_x.struphy.plot.slice(x="z", y="t", symmetric=True, backend="plotly")`
   or `spectrum.struphy.plot.dispersion(kmin=0, branches=..., fits=..., backend="plotly")` (or with
-  `plotly.graph_objects` for anything struphy-plots does not draw), and save each with the script's
-  own `save(figure, name, show=show)`, copied from `maxwell-wave.py`, which writes `name.html`,
-  `.png` and `.plotly.json` on MPI rank 0. Name the page's main figure `<script-stem>` and the
-  others `<script-stem>-<key>`.
+  `plotly.graph_objects` for anything struphy-plots does not draw), and save each with
+  `struphy_plots.save_figure(figure, name, show=show)`, which writes `name.html`, `.png` and
+  `.plotly.json` on MPI rank 0 (`frame=` or `still=` choose the image of an animation). Name the
+  page's main figure `<script-stem>` and the others `<script-stem>-<key>`.
 - The website's part is `run_example.py`, which `python cli.py run` and CI call: it runs the script
   in `docs/public/examples/`, profiles its simulation, exports the profiling data into
   `<script-stem>.metadata.json` and copies the PNGs to `docs/public/images/examples/`.
@@ -147,7 +147,7 @@ GitHub-hosted runners execute each example with `run_example.py` in a single pro
 For optional local MPI runs: `python cli.py run <example> --mpi 4`, or, from the repository root,
 `mpirun -n 4 python run_example.py <script-stem>`.
 
-- The simulation, `output.pproc(...)` and the analysis run on every rank; the scripts' `save`
+- The simulation, `output.pproc(...)` and the analysis run on every rank; `save_figure`
   writes on rank 0 only (struphy-plots draws nothing on the other ranks), as does `run_example.py`.
 - The grid must split over the ranks: with four ranks, keep at least a few cells per rank and direction
   (Orszag–Tang uses 32 × 32 × 1, i.e. 16 × 16 cells per rank).
diff --git a/docs/src/examples/acoustic-pulse.py b/docs/src/examples/acoustic-pulse.py
index 01e4d02..d079dfc 100644
--- a/docs/src/examples/acoustic-pulse.py
+++ b/docs/src/examples/acoustic-pulse.py
@@ -30,6 +30,7 @@
 )
 from struphy.linear_algebra.solver import NonlinearSolverParameters
 from struphy.models import VariationalCompressibleFluid
+from struphy_plots import save_figure
 
 length = 2 * np.pi
 width = 0.5  # of the Gaussian pulse
@@ -104,25 +105,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -176,7 +158,7 @@ def profile_traces(index):
     figure.update_yaxes(title_text="density", range=[1 - 0.3 * amplitude, 1 + 1.1 * amplitude], row=1, col=1)
     figure.update_xaxes(title_text="t", row=2, col=1)
     figure.update_yaxes(title_text="energy change / largest kinetic energy", row=2, col=1)
-    save(figure, "acoustic-pulse", width=900, height=850, show=show)
+    save_figure(figure, "acoustic-pulse", width=900, height=850, show=show)
 
     density_change = rho.copy(data=density - 1.0)
     space_time = density_change.assign_coords(eta1=x).struphy.plot.slice(
@@ -190,7 +172,7 @@ def profile_traces(index):
         colorbar_label="ρ − 1",
         backend="plotly",
     )
-    save(space_time, "acoustic-pulse-space-time", show=show)
+    save_figure(space_time, "acoustic-pulse-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/alfven-standing-wave.py b/docs/src/examples/alfven-standing-wave.py
index cbafd55..ad3cf04 100644
--- a/docs/src/examples/alfven-standing-wave.py
+++ b/docs/src/examples/alfven-standing-wave.py
@@ -22,6 +22,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import ShearAlfven
+from struphy_plots import save_figure
 
 # A uniform background with B0 along z and n0 = 1, so the Alfvén speed is 1.
 B0z, n0, beta = 1.0, 1.0, 0.1
@@ -70,27 +71,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc(physical=True)
@@ -132,7 +112,7 @@ def pproc(sim: Simulation, show: bool = False):
         margin={"l": 75, "r": 30, "t": 80, "b": 100},
     )
     figure.update_yaxes(range=[-0.05, 1.15])
-    save(figure, "alfven-standing-wave", show=show)
+    save_figure(figure, "alfven-standing-wave", show=show)
 
     # The velocity along z over time: a standing wave keeps its nodes, so the stripes are vertical,
     # unlike the diagonal stripes of the travelling waves in `shear-alfven-wave`.
@@ -151,7 +131,7 @@ def pproc(sim: Simulation, show: bool = False):
         colorbar_label="u₁ (logical component)",
         backend="plotly",
     )
-    save(space_time, "alfven-standing-wave-space-time", show=show)
+    save_figure(space_time, "alfven-standing-wave-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/beltrami-sph.py b/docs/src/examples/beltrami-sph.py
index a78a29b..15ee7a5 100644
--- a/docs/src/examples/beltrami-sph.py
+++ b/docs/src/examples/beltrami-sph.py
@@ -41,6 +41,7 @@
 )
 from struphy.models import PressureLessSPH
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 box_min = -0.5
 box_max = 0.5
@@ -139,27 +140,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc()
@@ -359,7 +339,7 @@ def error_trace(index, values, name, color):
     # The frames name only the markers and errors (the contour stays), so build the final still.
     still = go.Figure(data=final_data, layout=figure.layout)
     still.layout.sliders[0].active = len(figure.frames) - 1
-    save(figure, "beltrami-sph", still=still, width=1100, height=680, show=show)
+    save_figure(figure, "beltrami-sph", still=still, width=1100, height=680, show=show)
 
     # The kernel reconstruction shows the simulated mass density independently of the marker view.
     # The exact divergence-free Beltrami transport preserves the initially uniform rho = 1.
@@ -594,10 +574,10 @@ def signed_cell_area(grid_x, grid_y):
         f"maximum marker-cell area deviation from 1 {max_area_deviation:.3e}"
     )
 
-    save(density_figure, "beltrami-sph-density", frame=-1, show=show)
-    save(compression_figure, "beltrami-sph-compression", frame=-1, show=show)
-    save(area_figure, "beltrami-sph-area-deformation", frame=-1, show=show)
-    save(trajectory_figure, "beltrami-sph-trajectories", show=show)
+    save_figure(density_figure, "beltrami-sph-density", frame=-1, show=show)
+    save_figure(compression_figure, "beltrami-sph-compression", frame=-1, show=show)
+    save_figure(area_figure, "beltrami-sph-area-deformation", frame=-1, show=show)
+    save_figure(trajectory_figure, "beltrami-sph-trajectories", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/bump-on-tail.py b/docs/src/examples/bump-on-tail.py
index 2a45460..c56c5d2 100644
--- a/docs/src/examples/bump-on-tail.py
+++ b/docs/src/examples/bump-on-tail.py
@@ -31,6 +31,7 @@
     perturbations,
 )
 from struphy.models import VlasovAmpereOneSpecies
+from struphy_plots import save_figure
 
 # A 90% bulk Maxwellian plus a 10% "bump" population drifting at u1 = -4.5.
 perturbation_amplitude = 0.05
@@ -91,27 +92,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     domain = sim.domain
 
@@ -125,7 +105,7 @@ def pproc(sim: Simulation, show: bool = False):
 
     figure = field_energy.struphy.plot.timeseries(logy=True, title="Bump-on-tail instability: electric field energy", backend="plotly")
 
-    save(figure, "bump-on-tail", show=show)
+    save_figure(figure, "bump-on-tail", show=show)
 
     # Evaluate the saved products on their grids, then look at them in more than one way.
     output.pproc()
@@ -162,8 +142,8 @@ def pproc(sim: Simulation, show: bool = False):
         backend="plotly",
     )
 
-    save(phase_space, "bump-on-tail-phasespace", frame=len(phase_space.fig.frames) // 2, show=show)
-    save(velocity_time, "bump-on-tail-velocity-time", show=show)
+    save_figure(phase_space, "bump-on-tail-phasespace", frame=len(phase_space.fig.frames) // 2, show=show)
+    save_figure(velocity_time, "bump-on-tail-velocity-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/coaxial-waveguide.py b/docs/src/examples/coaxial-waveguide.py
index febbdb2..0670fab 100644
--- a/docs/src/examples/coaxial-waveguide.py
+++ b/docs/src/examples/coaxial-waveguide.py
@@ -23,6 +23,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import Maxwell
+from struphy_plots import save_figure
 
 # The inner and outer radii are two zeros of the derivative of the radial profile
 # J_m(r) - 0.28 Y_m(r), so that the tangential electric field vanishes on both conducting walls.
@@ -75,27 +76,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from struphy_plots.analysis import evaluate_on
 
@@ -243,7 +223,7 @@ def frame_traces(index):
 
     # The still image and thumbnail show the mode a little after the start of the run.
     still_position = len(frame_indices) // 4
-    save(figure, "coaxial-waveguide", height=650, frame=still_position, show=show)
+    save_figure(figure, "coaxial-waveguide", height=650, frame=still_position, show=show)
 
     # The probe signal against the exact mode.
     probe_figure = signal.struphy.plot.timeseries(
@@ -291,8 +271,8 @@ def frame_traces(index):
     energy_figure.update_xaxes(title_text="t [a.u.]", row=2, col=1)
     energy_figure.update_layout(template="plotly_white", autosize=True, margin={"l": 80, "r": 30, "t": 80, "b": 60})
 
-    save(probe_figure, "coaxial-waveguide-frequency", show=show)
-    save(energy_figure, "coaxial-waveguide-energy", show=show)
+    save_figure(probe_figure, "coaxial-waveguide-frequency", show=show)
+    save_figure(energy_figure, "coaxial-waveguide-energy", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/cold-plasma-oscillation.py b/docs/src/examples/cold-plasma-oscillation.py
index 2176fa4..fea1371 100644
--- a/docs/src/examples/cold-plasma-oscillation.py
+++ b/docs/src/examples/cold-plasma-oscillation.py
@@ -19,6 +19,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import ColdPlasma
+from struphy_plots import save_figure
 
 # Time in units of the inverse cyclotron frequency and alpha = omega_p / omega_c = 1: at n0 = 1 the plasma
 # frequency is 1 (as in the cold-plasma-waves example).
@@ -74,27 +75,6 @@ def oscillation_frequency(times, values):
     return float(np.pi / np.mean(np.diff(roots)))
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -160,7 +140,7 @@ def profile_traces(index):
     figure.update_yaxes(title_text="E_z", range=[-1.2 * amplitude, 1.2 * amplitude], row=1, col=1)
     figure.update_xaxes(title_text="t", row=2, col=1)
     figure.update_yaxes(title_text="energy / initial energy", row=2, col=1)
-    save(figure, "cold-plasma-oscillation", width=900, height=850, show=show)
+    save_figure(figure, "cold-plasma-oscillation", width=900, height=850, show=show)
 
     # The frequency against the density, on the line omega = omega_p (drawn from n0 = 0).
     n_line = np.linspace(0.0, max(densities) * 1.05, 100)
@@ -170,7 +150,7 @@ def profile_traces(index):
         xlabel="density n₀", ylabel="angular frequency ω", title="Oscillation frequency against density",
         backend="plotly",
     )
-    save(scan, "cold-plasma-oscillation-frequency-scan", show=show)
+    save_figure(scan, "cold-plasma-oscillation-frequency-scan", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/cold-plasma-wave-packet.py b/docs/src/examples/cold-plasma-wave-packet.py
index 316e09b..eb40b9d 100644
--- a/docs/src/examples/cold-plasma-wave-packet.py
+++ b/docs/src/examples/cold-plasma-wave-packet.py
@@ -21,6 +21,7 @@
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids
 from struphy.initial.base import GenericPerturbation
 from struphy.models import ColdPlasma
+from struphy_plots import save_figure
 
 # Plasma frequency equal to the cyclotron frequency, time in units of the inverse cyclotron frequency, c = 1.
 alpha, epsilon, n0, B0z = 1.0, 1.0, 1.0, 1.0
@@ -80,27 +81,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from struphy_plots.theory.waves import Species, cold_plasma_waves, group_velocity
 
@@ -169,7 +149,7 @@ def miss(side, speed):
     )
     figure.update_xaxes(range=[0.0, length])
     figure.update_yaxes(range=[0.0, time_opts.Tend])
-    save(figure, "cold-plasma-wave-packet", width=1100, height=750, show=show)
+    save_figure(figure, "cold-plasma-wave-packet", width=1100, height=750, show=show)
 
     centroids = go.Figure()
     colors = {"whistler": "#168aad", "R wave": "#2a9d8f", "L wave": "#d62828"}
@@ -186,7 +166,7 @@ def miss(side, speed):
     )
     centroids.update_xaxes(range=[0.0, time_opts.Tend])
     centroids.update_yaxes(range=[0.0, center])
-    save(centroids, "cold-plasma-wave-packet-centroid", show=show)
+    save_figure(centroids, "cold-plasma-wave-packet-centroid", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/cold-plasma-waves.py b/docs/src/examples/cold-plasma-waves.py
index a7f42d4..d61b40a 100644
--- a/docs/src/examples/cold-plasma-waves.py
+++ b/docs/src/examples/cold-plasma-waves.py
@@ -17,6 +17,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import ColdPlasma
+from struphy_plots import save_figure
 
 # Plasma frequency equal to the cyclotron frequency (alpha = 1), and time in units of the inverse
 # cyclotron frequency (epsilon = 1). The R cutoff is then at (1 + sqrt 5)/2 and the L cutoff at (sqrt 5 - 1)/2.
@@ -83,27 +84,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc(physical=True)
@@ -156,7 +136,7 @@ def pproc(sim: Simulation, show: bool = False):
         title="Cold-plasma waves along B₀: power spectrum of E_x",
         backend="plotly",
     )
-    save(figure, "cold-plasma-waves", height=700, show=show)
+    save_figure(figure, "cold-plasma-waves", height=700, show=show)
 
     # Energy channels: the noise starts purely electric, then shares its energy with the magnetic field
     # and the electron current, while the sum stays constant.
@@ -166,7 +146,7 @@ def pproc(sim: Simulation, show: bool = False):
     energy_figure = energies[0].struphy.plot.timeseries(
         *energies[1:], logy=False, title="Energy channels of the cold plasma", backend="plotly"
     )
-    save(energy_figure, "cold-plasma-waves-energy", show=show)
+    save_figure(energy_figure, "cold-plasma-waves-energy", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/dam-break.py b/docs/src/examples/dam-break.py
index 15df174..d0ace5b 100644
--- a/docs/src/examples/dam-break.py
+++ b/docs/src/examples/dam-break.py
@@ -34,6 +34,7 @@
 )
 from struphy.models import ViscousEulerSPH
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 # Weakly compressible SPH: the isothermal pressure p = kappa * rho, with kappa = c_s^2 small
 # enough that the flow stays roughly subsonic, and a little viscosity for stability.
@@ -107,27 +108,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc()
@@ -243,7 +223,7 @@ def frame_traces(index, webgl=True):
     still_index = int(np.argmin(abs(times - 0.5)))
     still = go.Figure(data=frame_traces(still_index, webgl=False), layout=figure.layout)
     still.layout.sliders[0].active = still_index
-    save(figure, "dam-break", height=750, still=still, show=show)
+    save_figure(figure, "dam-break", height=750, still=still, show=show)
 
     trajectory = go.Figure()
     trajectory.add_scatter(
@@ -269,7 +249,7 @@ def frame_traces(index, webgl=True):
         legend={"x": 0.98, "y": 0.5, "xanchor": "right", "bgcolor": "rgba(255,255,255,0.82)"},
         margin={"l": 70, "r": 30, "t": 80, "b": 60},
     )
-    save(trajectory, "dam-break-front", show=show)
+    save_figure(trajectory, "dam-break-front", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/damped-alfven-wave.py b/docs/src/examples/damped-alfven-wave.py
index b45191b..45c1aed 100644
--- a/docs/src/examples/damped-alfven-wave.py
+++ b/docs/src/examples/damped-alfven-wave.py
@@ -29,6 +29,7 @@
 )
 from struphy.linear_algebra.solver import NonlinearSolverParameters
 from struphy.models import ViscoResistiveMHD
+from struphy_plots import save_figure
 
 length = 2 * np.pi
 wavenumber = 2 * np.pi / length  # one wavelength in the box
@@ -100,27 +101,6 @@ def mode_amplitude(run):
     return u_y.struphy.analysis.project_mode(dim="eta1", number=1)
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -184,7 +164,7 @@ def profile_traces(index):
     figure.update_yaxes(title_text="u_y", range=[-1.2 * amplitude, 1.2 * amplitude], row=1, col=1)
     figure.update_xaxes(title_text="t", row=2, col=1)
     figure.update_yaxes(title_text="amplitude of sin(kx)", row=2, col=1)
-    save(figure, "damped-alfven-wave", width=900, height=850, show=show)
+    save_figure(figure, "damped-alfven-wave", width=900, height=850, show=show)
 
     colors = {0.05: "#168aad", 0.1: "#d62828", 0.2: "#f77f00"}
     decay = go.Figure()
@@ -199,7 +179,7 @@ def profile_traces(index):
         xaxis_title="t", yaxis_title="amplitude of the sin(kx) mode", yaxis_type="log",
         margin={"l": 75, "r": 30, "t": 80, "b": 60},
     )
-    save(decay, "damped-alfven-wave-decay", show=show)
+    save_figure(decay, "damped-alfven-wave-decay", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/diffusion-methods.py b/docs/src/examples/diffusion-methods.py
index 06e5c4a..e81ed81 100644
--- a/docs/src/examples/diffusion-methods.py
+++ b/docs/src/examples/diffusion-methods.py
@@ -39,6 +39,7 @@
 )
 from struphy.models import DeterministicParticleDiffusion, RandomParticleDiffusion
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 diffusion = 0.05  # D
 amplitude = 0.5  # of the initial density perturbation, relative to the uniform density 1
@@ -101,25 +102,6 @@ def create_simulation(method="Random walk") -> Simulation:
     )
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -188,7 +170,7 @@ def traces(index):
     figure.update_yaxes(title_text="density", range=[0.4, 1.6], row=1, col=1)
     figure.update_xaxes(title_text="t", row=2, col=1)
     figure.update_yaxes(title_text="amplitude", type="log", row=2, col=1)
-    save(figure, "diffusion-methods", width=900, height=850, show=show)
+    save_figure(figure, "diffusion-methods", width=900, height=850, show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/diocotron-instability.py b/docs/src/examples/diocotron-instability.py
index e534a4f..5834f86 100644
--- a/docs/src/examples/diocotron-instability.py
+++ b/docs/src/examples/diocotron-instability.py
@@ -40,6 +40,7 @@
     perturbations,
 )
 from struphy.models import ToyDrift
+from struphy_plots import save_figure
 
 # A uniform-density ring between r = 4 and r = 5, seeded with a tiny m = 4 azimuthal mode.
 r_minus, r_plus, mode_number = 4.0, 5.0, 4
@@ -107,27 +108,6 @@ def ring_density(etas, r_minus=r_minus, r_plus=r_plus):
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     a1, a2 = sim.domain.params["a1"], sim.domain.params["a2"]
 
@@ -205,7 +185,7 @@ def pproc(sim: Simulation, show: bool = False):
         ],
     )
 
-    save(figure, "diocotron-instability", height=800, show=show)
+    save_figure(figure, "diocotron-instability", height=800, show=show)
 
     mode_figure = go.Figure()
     for m, amplitude in mode_amplitudes.items():
@@ -274,8 +254,8 @@ def curve(values, name, color):
         sliders=[{"steps": [{"args": [[frame.name], {"frame": {"duration": 0, "redraw": True}, "mode": "immediate"}], "label": frame.name, "method": "animate"} for frame in interface_frames], "active": len(interface_frames) - 1, "x": 0.12, "len": 0.88, "y": -0.18, "currentvalue": {"prefix": "t = "}}],
     )
 
-    save(mode_figure, "diocotron-instability-mode-growth", show=show)
-    save(interface_figure, "diocotron-instability-ring-interfaces", show=show)
+    save_figure(mode_figure, "diocotron-instability-mode-growth", show=show)
+    save_figure(interface_figure, "diocotron-instability-ring-interfaces", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/faraday-rotation.py b/docs/src/examples/faraday-rotation.py
index c147146..7e4f96f 100644
--- a/docs/src/examples/faraday-rotation.py
+++ b/docs/src/examples/faraday-rotation.py
@@ -19,6 +19,7 @@
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import ColdPlasma
 from struphy.linear_algebra.solver import SolverParameters
+from struphy_plots import save_figure
 
 stem = "faraday-rotation"
 length, amplitude, omega, density = 2.0 * np.pi, 0.05, 2.5, 3.15
@@ -65,27 +66,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     import plotly.graph_objects as go
@@ -145,7 +125,7 @@ def pproc(sim: Simulation, show: bool = False):
     )
     figure.update_xaxes(title_text="z", row=1, col=2)
     figure.update_yaxes(title_text="polarization angle [degrees]", row=1, col=2)
-    save(figure, stem, height=700, show=show)
+    save_figure(figure, stem, height=700, show=show)
     probes = go.Figure()
     for target in (0.0, length / 4, length / 2):
         j = int(np.argmin(np.abs(z - target)))
@@ -153,7 +133,7 @@ def pproc(sim: Simulation, show: bool = False):
     probes.update_layout(title="Local polarization remains linear", template="plotly_white",
                           xaxis_title="E_x", yaxis_title="E_y", yaxis={"scaleanchor": "x"},
                           margin={"l": 70, "r": 30, "t": 80, "b": 70})
-    save(probes, f"{stem}-polarization", show=show)
+    save_figure(probes, f"{stem}-polarization", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/gas-expansion.py b/docs/src/examples/gas-expansion.py
index e174416..a334665 100644
--- a/docs/src/examples/gas-expansion.py
+++ b/docs/src/examples/gas-expansion.py
@@ -33,6 +33,7 @@
 )
 from struphy.models import ViscousEulerSPH
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 # The gas is isothermal, p = kappa * rho, so the sound speed is c = sqrt(kappa) = 1.
 kappa = 1.0
@@ -111,27 +112,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -270,7 +250,7 @@ def frame_traces(index, webgl=False):
         ],
     )
     still_index = int(np.argmin(np.abs(times - 0.8)))
-    save(figure, "gas-expansion", height=750, frame=still_index, show=show)
+    save_figure(figure, "gas-expansion", height=750, frame=still_index, show=show)
 
     # A spatial view of every SPH particle. Fixed display lanes separate overlapping particles;
     # only x is a physical coordinate in this one-dimensional simulation.
@@ -332,7 +312,7 @@ def particle_trace(index):
     )
     particles.add_vline(x=release_point, line_dash="dash", line_color="#64748b",
                         annotation_text="initial gas edge", annotation_position="top right")
-    save(particles, "gas-expansion-particles", show=show)
+    save_figure(particles, "gas-expansion-particles", show=show)
 
     # The same data against the similarity variable xi = (x - x0) / t, where every time falls on one curve.
     similarity = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1)
@@ -420,8 +400,8 @@ def particle_trace(index):
         margin={"l": 70, "r": 30, "t": 80, "b": 60},
     )
 
-    save(similarity, "gas-expansion-similarity", show=show)
-    save(error_figure, "gas-expansion-error", show=show)
+    save_figure(similarity, "gas-expansion-similarity", show=show)
+    save_figure(error_figure, "gas-expansion-error", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/grad-b-drift.py b/docs/src/examples/grad-b-drift.py
index b2accfa..b552624 100644
--- a/docs/src/examples/grad-b-drift.py
+++ b/docs/src/examples/grad-b-drift.py
@@ -23,6 +23,7 @@
     SavingParameters, Simulation, Time, domains, equils, grids, maxwellians,
 )
 from struphy.models import Vlasov
+from struphy_plots import save_figure
 
 stem = "grad-b-drift"
 box, ripple = 40.0, 0.3
@@ -79,27 +80,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     equil = sim.equil
     time_opts = sim.time_opts
@@ -155,7 +135,7 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_yaxes(title_text="mean drift velocity in y", row=1, col=2)
     figure.update_layout(title="Grad-B drift in a straight magnetic field", template="plotly_white",
                          legend={"orientation": "h", "y": -0.2}, margin={"l": 70, "r": 35, "t": 100, "b": 120})
-    save(figure, stem, height=650, show=show)
+    save_figure(figure, stem, height=650, show=show)
     drift = go.Figure()
     for j, color in enumerate(palette):
         drift.add_scatter(x=times, y=center_y[:, j] - center_y[0, j], name=f"v⊥ = {speed[0, j]:.1f}", line={"color": color})
@@ -164,7 +144,7 @@ def pproc(sim: Simulation, show: bool = False):
     drift.update_layout(title="Guiding-center displacement extracted from full orbits", template="plotly_white",
                          xaxis_title="t", yaxis_title="Y_gc(t) − Y_gc(0)",
                          legend={"orientation": "h", "y": -0.2}, margin={"l": 80, "r": 30, "t": 80, "b": 110})
-    save(drift, f"{stem}-drift", show=show)
+    save_figure(drift, f"{stem}-drift", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/guiding-center-orbits.py b/docs/src/examples/guiding-center-orbits.py
index 415278c..5c513bf 100644
--- a/docs/src/examples/guiding-center-orbits.py
+++ b/docs/src/examples/guiding-center-orbits.py
@@ -35,6 +35,7 @@
     maxwellians,
 )
 from struphy.models import GuidingCenter
+from struphy_plots import save_figure
 
 # Eight markers on the outboard midplane of the flux surface eta1 = 0.5, all with speed 3 and
 # different pitch. The state of a guiding center is (position, v_parallel, mu), with the magnetic
@@ -95,27 +96,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     domain = sim.domain
     equil = sim.equil
@@ -307,7 +287,7 @@ def moving_traces(index):
     # (A frame names only the moving traces, so the still is built from them and the fixed ones.)
     still = go.Figure(data=list(figure.data[:first_moving]) + moving_traces(len(times) - 1), layout=figure.layout)
     still.layout.sliders[0].active = len(frames) - 1
-    save(figure, "guiding-center-orbits", height=750, still=still, show=show)
+    save_figure(figure, "guiding-center-orbits", height=750, still=still, show=show)
 
     # One panel per particle: the orbit in the poloidal plane.
     panels = make_subplots(rows=2, cols=4, subplot_titles=labels, horizontal_spacing=0.03, vertical_spacing=0.12)
@@ -403,9 +383,9 @@ def moving_traces(index):
         margin={"l": 70, "r": 30, "t": 100, "b": 60},
     )
 
-    save(panels, "guiding-center-orbits-panels", show=show)
-    save(velocity, "guiding-center-orbits-velocity", show=show)
-    save(conservation, "guiding-center-orbits-conservation", show=show)
+    save_figure(panels, "guiding-center-orbits-panels", show=show)
+    save_figure(velocity, "guiding-center-orbits-velocity", show=show)
+    save_figure(conservation, "guiding-center-orbits-conservation", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/gvec-equilibrium.py b/docs/src/examples/gvec-equilibrium.py
index 3e1c5aa..bba82e7 100644
--- a/docs/src/examples/gvec-equilibrium.py
+++ b/docs/src/examples/gvec-equilibrium.py
@@ -40,6 +40,7 @@
     maxwellians,
 )
 from struphy.models import GuidingCenter
+from struphy_plots import save_figure
 
 # Keep both the GVEC solve and the following FEEC simulation deliberately
 # small. The complete equilibrium input is defined below, so the example has
@@ -238,27 +239,6 @@ def create_simulation(use_gvec=False) -> Simulation:
 metadata_overrides = {"domain": "GVECunit (generated by pyGVEC)"}
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -389,7 +369,7 @@ def pproc(sim: Simulation, show: bool = False):
         legend={"x": 0.01, "y": 0.99, "bgcolor": "rgba(255,255,255,0.72)"},
         margin={"l": 10, "r": 30, "t": 75, "b": 10},
     )
-    save(figure, "gvec-equilibrium", height=760, show=show)
+    save_figure(figure, "gvec-equilibrium", height=760, show=show)
 
     # A poloidal slice reveals the nesting more quantitatively; the adjacent
     # radial profiles come from the very same state file used by Struphy.
@@ -469,7 +449,7 @@ def pproc(sim: Simulation, show: bool = False):
         legend={"orientation": "h", "y": -0.2},
         margin={"l": 70, "r": 70, "t": 90, "b": 100},
     )
-    save(profiles, "gvec-equilibrium-flux-surfaces", show=show)
+    save_figure(profiles, "gvec-equilibrium-flux-surfaces", show=show)
 
     # One panel per marker, the orbit projected on the (R, Z) plane, as in the tokamak example. The
     # cross-section of a stellarator turns with the toroidal angle, so the surfaces drawn behind each
@@ -544,7 +524,7 @@ def pproc(sim: Simulation, show: bool = False):
         margin={"l": 70, "r": 30, "t": 90, "b": 60},
         height=760,
     )
-    save(panels, "gvec-equilibrium-poloidal-orbits", show=show)
+    save_figure(panels, "gvec-equilibrium-poloidal-orbits", show=show)
 
     diagnostics = make_subplots(
         rows=1,
@@ -584,7 +564,7 @@ def pproc(sim: Simulation, show: bool = False):
         legend={"orientation": "h", "y": -0.2},
         margin={"l": 75, "r": 30, "t": 80, "b": 100},
     )
-    save(diagnostics, "gvec-equilibrium-orbit-diagnostics", show=show)
+    save_figure(diagnostics, "gvec-equilibrium-orbit-diagnostics", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/gyromotion.py b/docs/src/examples/gyromotion.py
index c81f7b5..ca96cbb 100644
--- a/docs/src/examples/gyromotion.py
+++ b/docs/src/examples/gyromotion.py
@@ -20,6 +20,7 @@
 from struphy.kinetic_background import maxwellians
 from struphy.models import Vlasov
 from struphy.pic.base import BoundaryParameters, LoadingParameters, SavingParameters
+from struphy_plots import save_figure
 
 B0 = 1.0  # q = m = 1, so the gyrofrequency is B0
 gyrofrequency = B0
@@ -71,25 +72,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     from plotly.subplots import make_subplots
@@ -165,7 +147,7 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_yaxes(title_text="y", range=[centre - 2 * span, centre + 0.4 * span], scaleanchor="x", row=1, col=1)
     figure.update_xaxes(title_text="t", row=1, col=2)
     figure.update_yaxes(title_text="|x − x_exact|", type="log", row=1, col=2)
-    save(figure, "gyromotion", width=1100, height=620, show=show)
+    save_figure(figure, "gyromotion", width=1100, height=620, show=show)
 
     helix = go.Figure()
     for p in range(n_markers):
@@ -188,8 +170,8 @@ def pproc(sim: Simulation, show: bool = False):
         xaxis_title="t", yaxis_title="relative change of speed (solid) and perpendicular speed (dotted)", yaxis_type="log",
         margin={"l": 75, "r": 30, "t": 80, "b": 60},
     )
-    save(helix, "gyromotion-helices", show=show)
-    save(conservation, "gyromotion-conservation", show=show)
+    save_figure(helix, "gyromotion-helices", show=show)
+    save_figure(conservation, "gyromotion-conservation", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/hall-mhd-waves.py b/docs/src/examples/hall-mhd-waves.py
index 8c317dd..540dd93 100644
--- a/docs/src/examples/hall-mhd-waves.py
+++ b/docs/src/examples/hall-mhd-waves.py
@@ -20,6 +20,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import LinearExtendedMHDuniform
+from struphy_plots import save_figure
 
 # Background field along z with n0 = B0 = 1, so the Alfvén speed is 1. With epsilon = 1 the ion cyclotron
 # frequency and the ion inertial length are 1 as well, i.e. k is measured in units of 1/d_i.
@@ -74,27 +75,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -194,7 +174,7 @@ def log_power(spectrum):
         template="plotly_white", legend={"orientation": "h", "y": -0.2},
         margin={"l": 70, "r": 40, "t": 90, "b": 110},
     )
-    save(figure, "hall-mhd-waves", width=1300, height=650, show=show)
+    save_figure(figure, "hall-mhd-waves", width=1300, height=650, show=show)
 
     # Phase velocities: the ideal-MHD Alfvén wave is not dispersive (omega / k = v_A for every k); the
     # Hall branches are, one speeding up and one slowing down as k d_i grows.
@@ -217,7 +197,7 @@ def log_power(spectrum):
         margin={"l": 75, "r": 30, "t": 80, "b": 120},
     )
     phase.update_xaxes(range=[0, k_top])
-    save(phase, "hall-mhd-waves-phase-velocity", show=show)
+    save_figure(phase, "hall-mhd-waves-phase-velocity", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/hasegawa-wakatani.py b/docs/src/examples/hasegawa-wakatani.py
index f09c7c4..8f89f7e 100644
--- a/docs/src/examples/hasegawa-wakatani.py
+++ b/docs/src/examples/hasegawa-wakatani.py
@@ -23,6 +23,7 @@
 from struphy.initial.base import GenericPerturbation
 from struphy.models import HasegawaWakatani
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 # A square periodic slab.  The third direction is inactive, making this a
 # genuinely 2D fluid calculation rather than a thin 3D one.
@@ -108,27 +109,6 @@ def initial_vorticity(x, y, z):
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     env = sim.env
     time_opts = sim.time_opts
@@ -210,7 +190,7 @@ def field_traces(index, *, colorbars=False):
     for column in (1, 2):
         figure.update_xaxes(title_text="x", range=[0, length], constrain="domain", row=1, col=column)
         figure.update_yaxes(title_text="y", range=[0, length], scaleanchor="x" if column == 1 else "x2", scaleratio=1, row=1, col=column)
-    save(figure, "hasegawa-wakatani", width=1100, height=660, show=show)
+    save_figure(figure, "hasegawa-wakatani", width=1100, height=660, show=show)
 
     # E×B kinetic energy is |grad(phi)|²/2.  Averaging phi over y selects
     # ky = 0; its remaining y-directed velocity is the zonal flow.  Spectral
@@ -243,7 +223,7 @@ def field_traces(index, *, colorbars=False):
         margin={"l": 80, "r": 30, "t": 80, "b": 60},
         legend={"orientation": "h", "y": 1.1},
     )
-    save(energy_figure, "hasegawa-wakatani-zonal-energy", show=show)
+    save_figure(energy_figure, "hasegawa-wakatani-zonal-energy", show=show)
     zonal_fraction = zonal_energy / np.maximum(total_energy, np.finfo(float).tiny)
     print(f"Zonal-flow share of the kinetic energy: {zonal_fraction[-1]:.3f} at the end, {zonal_fraction.max():.3f} at its peak")
 
diff --git a/docs/src/examples/hybrid-alfven-ion-coupling.py b/docs/src/examples/hybrid-alfven-ion-coupling.py
index a540808..6515f24 100644
--- a/docs/src/examples/hybrid-alfven-ion-coupling.py
+++ b/docs/src/examples/hybrid-alfven-ion-coupling.py
@@ -35,6 +35,7 @@
     perturbations,
 )
 from struphy.models import LinearMHDVlasovPC
+from struphy_plots import save_figure
 
 
 def create_simulation() -> Simulation:
@@ -98,27 +99,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc()  # on every rank; under MPI the output is not processed implicitly
@@ -166,7 +146,7 @@ def pproc(sim: Simulation, show: bool = False):
         margin={"l": 70, "r": 70, "t": 80, "b": 60},
     )
 
-    save(figure, "hybrid-alfven-ion-coupling", show=show)
+    save_figure(figure, "hybrid-alfven-ion-coupling", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/hybrid-current-coupling.py b/docs/src/examples/hybrid-current-coupling.py
index bd253cf..56ad5bd 100644
--- a/docs/src/examples/hybrid-current-coupling.py
+++ b/docs/src/examples/hybrid-current-coupling.py
@@ -38,6 +38,7 @@
 )
 from struphy.linear_algebra.solver import SolverParameters
 from struphy.models import LinearMHDVlasovCC
+from struphy_plots import save_figure
 
 stem = "hybrid-current-coupling"
 length = 20.0
@@ -100,27 +101,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     import plotly.graph_objects as go
@@ -190,7 +170,7 @@ def pproc(sim: Simulation, show: bool = False):
         legend={"orientation": "h", "y": -0.14},
         margin={"l": 85, "r": 35, "t": 100, "b": 130},
     )
-    save(figure, stem, height=800, show=show)
+    save_figure(figure, stem, height=800, show=show)
 
     # Animate physical field snapshots with fixed axes so amplitude changes remain visible.
     # Keep the first and last snapshots, with at most 101 frames for a compact download.
@@ -254,7 +234,7 @@ def wave_traces(index):
             ],
         }],
     )
-    save(animation, f"{stem}-wave-animation", show=show)
+    save_figure(animation, f"{stem}-wave-animation", show=show)
 
     wave = velocity.assign_coords(eta3=velocity.eta3.values * length).struphy.plot.slice(
         x="eta3",
@@ -267,7 +247,7 @@ def wave_traces(index):
         colorbar_label="U_x",
         backend="plotly",
     )
-    save(wave, f"{stem}-space-time", show=show)
+    save_figure(wave, f"{stem}-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/incompressible-shear-relaxation.py b/docs/src/examples/incompressible-shear-relaxation.py
index 1b987a4..d033f97 100644
--- a/docs/src/examples/incompressible-shear-relaxation.py
+++ b/docs/src/examples/incompressible-shear-relaxation.py
@@ -38,6 +38,7 @@
 from struphy.initial.base import GenericPerturbation
 from struphy.models import IncompressibleNavierStokesSPH
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 viscosity = 0.1  # mu
 height = 1.0  # H, the channel is 1 long in x and H high in y
@@ -107,25 +108,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -194,7 +176,7 @@ def traces(index):
     figure.update_yaxes(title_text="u_x", range=[-limit, limit], row=1, col=2)
     figure.update_xaxes(title_text="t", row=2, col=1)
     figure.update_yaxes(title_text="amplitude", type="log", row=2, col=1)
-    save(figure, "incompressible-shear-relaxation", width=1000, height=850, show=show)
+    save_figure(figure, "incompressible-shear-relaxation", width=1000, height=850, show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/itg-drift-wave.py b/docs/src/examples/itg-drift-wave.py
index 1960c91..e72d69f 100644
--- a/docs/src/examples/itg-drift-wave.py
+++ b/docs/src/examples/itg-drift-wave.py
@@ -43,6 +43,7 @@
 from struphy.initial.base import GenericPerturbation
 from struphy.models import DriftKineticElectrostaticAdiabatic
 from struphy.propagators import implicit_diffusion
+from struphy_plots import save_figure
 
 # A magnetized annular column: radius in [a1, a2], periodic in angle and length.
 a1, a2, length = 0.1, 14.5, 1506.759067
@@ -200,27 +201,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -240,7 +220,7 @@ def pproc(sim: Simulation, show: bool = False):
 
     figure = perturbation_energy.struphy.plot.timeseries(logy=True, title="ITG drift wave: density perturbation energy", backend="plotly")
 
-    save(figure, "itg-drift-wave", show=show)
+    save_figure(figure, "itg-drift-wave", show=show)
 
     # ---- further figures: the potential and its Fourier modes --------------------------------------------
     # the grid points of the cells and the periodic endpoints, which mode_amplitudes drops
@@ -268,7 +248,7 @@ def pproc(sim: Simulation, show: bool = False):
         colorbar_label="φ [a.u.]",
         backend="plotly",
     )
-    save(movie, "itg-drift-wave-potential", frame=still_index, show=show)
+    save_figure(movie, "itg-drift-wave-potential", frame=still_index, show=show)
 
     # Amplitude of each poloidal mode (radial rms at the seeded axial mode number) and its growth rate.
     amplitude = radial_rms(spectrum.sel(n=mode_toroidal)).sel(m=slice(1, MAX_POLOIDAL_MODE))
@@ -304,7 +284,7 @@ def pproc(sim: Simulation, show: bool = False):
         margin={"l": 80, "r": 30, "t": 90, "b": 60},
         legend={"orientation": "h", "y": -0.22, "x": 0.0, "xanchor": "left"},
     )
-    save(growth_figure, "itg-drift-wave-mode-growth", show=show)
+    save_figure(growth_figure, "itg-drift-wave-mode-growth", show=show)
 
     # The (m, n) spectrum at the first and last saved time.
     both = radial_rms(spectrum.isel(t=[0, -1])).sel(m=slice(0, MAX_POLOIDAL_MODE))
@@ -338,7 +318,7 @@ def pproc(sim: Simulation, show: bool = False):
         spectrum_figure.update_xaxes(title_text="poloidal m", row=1, col=column)
     spectrum_figure.update_yaxes(title_text="axial n", dtick=1, row=1, col=1)
     spectrum_figure.update_layout(title="ITG drift wave: Fourier spectrum of φ", template="plotly_white", margin={"l": 70, "r": 30, "t": 90, "b": 60})
-    save(spectrum_figure, "itg-drift-wave-spectrum", show=show)
+    save_figure(spectrum_figure, "itg-drift-wave-spectrum", show=show)
 
     # Radial structure of the seeded mode.
     interior = slice(1, -1)  # phi = 0 at the Dirichlet boundaries r = a1, a2, which a log axis cannot show
@@ -355,7 +335,7 @@ def pproc(sim: Simulation, show: bool = False):
         template="plotly_white",
         margin={"l": 80, "r": 30, "t": 80, "b": 60},
     )
-    save(radial_figure, "itg-drift-wave-radial-structure", show=show)
+    save_figure(radial_figure, "itg-drift-wave-radial-structure", show=show)
 
     # Where and when does the potential grow?
     rms = np.sqrt((phi**2).mean(("eta2", "eta3"))).isel(eta1=interior)
@@ -373,7 +353,7 @@ def pproc(sim: Simulation, show: bool = False):
             backend="plotly",
         )
     )
-    save(rms_map, "itg-drift-wave-radial-time", show=show)
+    save_figure(rms_map, "itg-drift-wave-radial-time", show=show)
 
     # The flux-surface-averaged (m = n = 0) density change: does the profile flatten?
     density = output.evaluate("diagnostics/rho", **points)
@@ -392,7 +372,7 @@ def pproc(sim: Simulation, show: bool = False):
         colorbar_label="δ⟨ρ⟩ [a.u.]",
         backend="plotly",
     )
-    save(profile_change, "itg-drift-wave-profile-change", show=show)
+    save_figure(profile_change, "itg-drift-wave-profile-change", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/langmuir-wave-dispersion.py b/docs/src/examples/langmuir-wave-dispersion.py
index bdb00dd..52ea5e4 100644
--- a/docs/src/examples/langmuir-wave-dispersion.py
+++ b/docs/src/examples/langmuir-wave-dispersion.py
@@ -36,6 +36,7 @@
     perturbations,
 )
 from struphy.models import VlasovAmpereOneSpecies
+from struphy_plots import save_figure
 
 wavenumbers = (0.3, 0.4, 0.5, 0.6)  # the scan; the example itself is the k = 0.5 run
 amplitude = 0.001
@@ -89,27 +90,6 @@ def mode_amplitude(run):
     return e_x.struphy.analysis.project_mode(dim="eta1", number=1)
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     runs = {0.5: sim.output}
     for k in wavenumbers:
@@ -168,7 +148,7 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_xaxes(title_text="wavenumber k λ_D")
     figure.update_yaxes(title_text="ω_r / ω_p", row=1, col=1)
     figure.update_yaxes(title_text="γ / ω_p", row=1, col=2)
-    save(figure, "langmuir-wave-dispersion", width=1200, height=560, show=show)
+    save_figure(figure, "langmuir-wave-dispersion", width=1200, height=560, show=show)
 
     colors = {0.3: "#168aad", 0.4: "#2a9d8f", 0.5: "#f77f00", 0.6: "#d62828"}
     signals = go.Figure()
@@ -188,7 +168,7 @@ def pproc(sim: Simulation, show: bool = False):
         xaxis_title="t [1/ω_p]", yaxis_title="amplitude of the sin(kx) mode of E_x",
         margin={"l": 75, "r": 30, "t": 80, "b": 60},
     )
-    save(signals, "langmuir-wave-dispersion-signals", show=show)
+    save_figure(signals, "langmuir-wave-dispersion-signals", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/linear-dissipative-alfven-wave.py b/docs/src/examples/linear-dissipative-alfven-wave.py
index 5b4d0cb..71721a3 100644
--- a/docs/src/examples/linear-dissipative-alfven-wave.py
+++ b/docs/src/examples/linear-dissipative-alfven-wave.py
@@ -20,6 +20,7 @@
     domains, equils, grids, perturbations,
 )
 from struphy.models import ViscoResistiveLinearMHD
+from struphy_plots import save_figure
 
 stem = "linear-dissipative-alfven-wave"
 length, amplitude, diffusivity = 2.0 * np.pi, 0.1, 0.1
@@ -58,27 +59,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     from plotly.subplots import make_subplots
@@ -130,12 +110,12 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_yaxes(title_text="energy / initial wave energy", row=2, col=1)
     figure.update_layout(title="Viscous and resistive damping of a linear Alfvén wave", template="plotly_white",
                          legend={"orientation": "h", "y": -0.18}, margin={"l": 85, "r": 30, "t": 90, "b": 140})
-    save(figure, stem, height=800, show=show)
+    save_figure(figure, stem, height=800, show=show)
     space_time = velocity.assign_coords(eta3=z).struphy.plot.slice(
         x="eta3", y="t", symmetric=True, cmap="RdBu_r", title="A standing Alfvén wave with a fading amplitude",
         xlabel="z", ylabel="t [a.u.]", colorbar_label="u_x", backend="plotly",
     )
-    save(space_time, f"{stem}-space-time", show=show)
+    save_figure(space_time, f"{stem}-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/maxwell-cavity-resonances.py b/docs/src/examples/maxwell-cavity-resonances.py
index 772f2aa..3641f4c 100644
--- a/docs/src/examples/maxwell-cavity-resonances.py
+++ b/docs/src/examples/maxwell-cavity-resonances.py
@@ -22,6 +22,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, grids, perturbations
 from struphy.models import Maxwell
+from struphy_plots import save_figure
 
 lx, ly = 1.0, 1.5
 
@@ -71,27 +72,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc(physical=True)
@@ -131,7 +111,7 @@ def pproc(sim: Simulation, show: bool = False):
     figure = spectrum.struphy.plot.power_spectrum(
         frequencies=exact_lines, omega_max=14.0, title="Resonances of a rectangular box", backend="plotly",
     )
-    save(figure, "maxwell-cavity-resonances", show=show)
+    save_figure(figure, "maxwell-cavity-resonances", show=show)
 
     error_figure = go.Figure(go.Scatter(
         x=[omega_exact for omega_exact, _ in exact], y=100.0 * np.asarray(errors), mode="markers+lines",
@@ -143,7 +123,7 @@ def pproc(sim: Simulation, show: bool = False):
         xaxis_title="exact ω [a.u.]", yaxis_title="(measured − exact) / exact [%]",
         margin={"l": 75, "r": 30, "t": 80, "b": 60},
     )
-    save(error_figure, "maxwell-cavity-resonances-frequency-error", show=show)
+    save_figure(error_figure, "maxwell-cavity-resonances-frequency-error", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/maxwell-curved-mesh.py b/docs/src/examples/maxwell-curved-mesh.py
index d02ee0b..af3986d 100644
--- a/docs/src/examples/maxwell-curved-mesh.py
+++ b/docs/src/examples/maxwell-curved-mesh.py
@@ -20,6 +20,7 @@
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, grids
 from struphy.initial.base import GenericPerturbation
 from struphy.models import Maxwell
+from struphy_plots import save_figure
 
 lx, ly = 2.0, 3.0
 distortion = 0.1
@@ -75,27 +76,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
     from scipy.interpolate import griddata
@@ -159,7 +139,7 @@ def resample(values):
     # The image shows the pulse at a fifth of the run, when the ring has spread; the page keeps the animation.
     still = go.Figure(data=figure.data, layout=figure.layout)
     still.data[0].z = resample(numeric[len(times) // 5])
-    save(figure, "maxwell-curved-mesh", width=750, height=1000, still=still, show=show)
+    save_figure(figure, "maxwell-curved-mesh", width=750, height=1000, still=still, show=show)
 
     diagnostics = make_subplots(rows=2, cols=1, vertical_spacing=0.18,
                                 subplot_titles=("Error of E_z against the exact solution", "Total energy"))
@@ -171,7 +151,7 @@ def resample(values):
     diagnostics.update_yaxes(title_text="rms error / initial rms", type="log", row=1, col=1)
     diagnostics.update_yaxes(title_text="relative change", exponentformat="e", row=2, col=1)
     diagnostics.update_xaxes(title_text="t", row=2, col=1)
-    save(diagnostics, "maxwell-curved-mesh-diagnostics", show=show)
+    save_figure(diagnostics, "maxwell-curved-mesh-diagnostics", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/maxwell-structure-preservation.py b/docs/src/examples/maxwell-structure-preservation.py
index b02eda3..240f15c 100644
--- a/docs/src/examples/maxwell-structure-preservation.py
+++ b/docs/src/examples/maxwell-structure-preservation.py
@@ -25,6 +25,7 @@
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, grids, perturbations
 from struphy.models import Maxwell
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 length = 20.0
 mode_number = 2
@@ -82,27 +83,6 @@ def e_x(run):
     ).isel(component=0)
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     runs = {}
     for label in schemes:
@@ -140,7 +120,7 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_yaxes(type="log", range=[-16.5, 2.0], exponentformat="power")
     figure.add_annotation(x=0.62, xref="paper", y=1.8, text="Heun 2 leaves the plot: unstable", showarrow=False,
                           font={"color": "#d62828"}, bgcolor="rgba(255,255,255,0.8)")
-    save(figure, "maxwell-structure-preservation", show=show)
+    save_figure(figure, "maxwell-structure-preservation", show=show)
 
     # What the drift does to the solution: the field profile at the end of the run against the exact one.
     stable = [label for label in runs if drifts[label] < 1.0]
@@ -163,7 +143,7 @@ def pproc(sim: Simulation, show: bool = False):
         legend={"orientation": "h", "y": -0.18},
         margin={"l": 80, "r": 30, "t": 80, "b": 100},
     )
-    save(profile, "maxwell-structure-preservation-profile", show=show)
+    save_figure(profile, "maxwell-structure-preservation-profile", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/maxwell-wave.py b/docs/src/examples/maxwell-wave.py
index 687a29b..ef718d5 100644
--- a/docs/src/examples/maxwell-wave.py
+++ b/docs/src/examples/maxwell-wave.py
@@ -23,6 +23,7 @@
     perturbations,
 )
 from struphy.models import Maxwell
+from struphy_plots import save_figure
 
 
 def create_simulation() -> Simulation:
@@ -68,27 +69,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     """Measure the speed of light from the (k, omega) spectrum, and save (and show) the figures."""
     from struphy_plots.analysis import power_spectrum
@@ -128,8 +108,8 @@ def pproc(sim: Simulation, show: bool = False):
         colorbar_label="E_x",
         backend="plotly",
     )
-    save(dispersion, "maxwell-wave", show=show)
-    save(space_time, "maxwell-wave-space-time", show=show)
+    save_figure(dispersion, "maxwell-wave", show=show)
+    save_figure(space_time, "maxwell-wave-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/mhd-slab-waves.py b/docs/src/examples/mhd-slab-waves.py
index d824e3b..6b41c2c 100644
--- a/docs/src/examples/mhd-slab-waves.py
+++ b/docs/src/examples/mhd-slab-waves.py
@@ -20,6 +20,7 @@
 from struphy import DerhamOptions, EnvironmentOptions, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import LinearMHD
 from struphy_plots.theory.waves import magnetosonic_speeds
+from struphy_plots import save_figure
 
 # The background: B0 = (0, 1, 1), density 0.7 and a plasma beta of 3 (thermal over magnetic pressure).
 B0x, B0y, B0z = 0.0, 1.0, 1.0
@@ -70,27 +71,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -163,7 +143,7 @@ def log_power(spectrum):
         title="Three MHD wave branches in a magnetized slab", template="plotly_white",
         legend={"orientation": "h", "y": -0.2}, margin={"l": 70, "r": 40, "t": 90, "b": 110},
     )
-    save(figure, "mhd-slab-waves", width=1300, height=650, show=show)
+    save_figure(figure, "mhd-slab-waves", width=1300, height=650, show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/ordinary-mode-dispersion.py b/docs/src/examples/ordinary-mode-dispersion.py
index c86f290..948d8b8 100644
--- a/docs/src/examples/ordinary-mode-dispersion.py
+++ b/docs/src/examples/ordinary-mode-dispersion.py
@@ -20,6 +20,7 @@
 from struphy.models import ColdPlasma
 from struphy.linear_algebra.solver import SolverParameters
 from struphy_plots.theory.waves import plasma_light_wave
+from struphy_plots import save_figure
 
 stem = "ordinary-mode-dispersion"
 length = 8.0 * np.pi
@@ -59,27 +60,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     from plotly.subplots import make_subplots
@@ -131,12 +111,12 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_yaxes(title_text="E-mode / A + vertical offset", row=1, col=2)
     figure.update_layout(title="Ordinary electromagnetic waves in a cold plasma", template="plotly_white",
                          legend={"orientation": "h", "y": -0.2}, margin={"l": 70, "r": 30, "t": 100, "b": 140})
-    save(figure, stem, height=650, show=show)
+    save_figure(figure, stem, height=650, show=show)
     space_time = field.assign_coords(eta1=x).struphy.plot.slice(
         x="eta1", y="t", symmetric=True, cmap="RdBu_r", title="Superposed ordinary waves: E_z(x, t)",
         xlabel="x", ylabel="t [a.u.]", colorbar_label="E_z", backend="plotly",
     )
-    save(space_time, f"{stem}-space-time", show=show)
+    save_figure(space_time, f"{stem}-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/orszag-tang-vortex.py b/docs/src/examples/orszag-tang-vortex.py
index f0e1975..c6c875d 100644
--- a/docs/src/examples/orszag-tang-vortex.py
+++ b/docs/src/examples/orszag-tang-vortex.py
@@ -16,6 +16,7 @@
 
 from struphy import DerhamOptions, EnvironmentOptions, FieldsBackground, Simulation, Time, domains, equils, grids, perturbations
 from struphy.models import ViscoResistiveMHD
+from struphy_plots import save_figure
 
 
 def create_simulation() -> Simulation:
@@ -63,27 +64,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     model = sim.model
     time_opts = sim.time_opts
@@ -145,7 +125,7 @@ def traces(index):
         sliders=[{"active": 0, "x": 0.12, "len": 0.88, "y": -0.07, "currentvalue": {"prefix": "t = "},
                   "steps": [{"args": [[frame.name], {"frame": {"duration": 0, "redraw": True}, "transition": {"duration": 0}, "mode": "immediate"}], "label": frame.name, "method": "animate"} for frame in figure.frames]}],
     )
-    save(figure, "orszag-tang-vortex", frame=-1, width=900, height=850, show=show)
+    save_figure(figure, "orszag-tang-vortex", frame=-1, width=900, height=850, show=show)
 
     scalars = output.scalars
     energy = scalars.en_tot
@@ -168,8 +148,8 @@ def traces(index):
         eta2=cut_index, backend="plotly",
     )
     print(f"Maximum relative total-energy drift: {float(drift.max()):.3e}; maximum ‖div B‖: {float(divergence.max()):.3e}")
-    save(diagnostics, "orszag-tang-vortex-conservation", show=show)
-    save(cut, "orszag-tang-vortex-pressure-cut", show=show)
+    save_figure(diagnostics, "orszag-tang-vortex-conservation", show=show)
+    save_figure(cut, "orszag-tang-vortex-pressure-cut", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/poisson-convergence.py b/docs/src/examples/poisson-convergence.py
index 6a2c3be..20123f0 100644
--- a/docs/src/examples/poisson-convergence.py
+++ b/docs/src/examples/poisson-convergence.py
@@ -21,6 +21,7 @@
 from struphy.initial.base import GenericPerturbation
 from struphy.linear_algebra.solver import SolverParameters
 from struphy.models import Poisson
+from struphy_plots import save_figure
 
 lx, ly = 2.0, 3.0
 kx, ky = 2 * np.pi / lx, 2 * np.pi / ly
@@ -75,27 +76,6 @@ def computed_potential(run, celldivide):
     return run.evaluate("em_fields/phi", **plane, eta3=0.0).isel(t=-1)
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
     from scipy.interpolate import griddata
@@ -143,7 +123,7 @@ def pproc(sim: Simulation, show: bool = False):
         xaxis_title="cells in x, n = 1 / h", yaxis_title="rms error of φ", xaxis_type="log", yaxis_type="log",
         legend={"x": 1.02, "y": 0.5}, margin={"l": 80, "r": 30, "t": 80, "b": 65},
     )
-    save(figure, "poisson-convergence", width=1300, height=650, show=show)
+    save_figure(figure, "poisson-convergence", width=1300, height=650, show=show)
 
     # The solution and its error on the distorted mesh, at degree 2 and 8 x 12 cells.
     run = create_simulation(2, 8, distortion, "poisson_convergence_map").output
@@ -168,7 +148,7 @@ def resample(values):
         maps.update_xaxes(title_text="x", range=[0, lx], constrain="domain", row=1, col=column)
         maps.update_yaxes(title_text="y", range=[0, ly], scaleanchor=f"x{column}" if column > 1 else "x", row=1, col=column)
     maps.update_layout(template="plotly_white", autosize=True, margin={"l": 70, "r": 30, "t": 80, "b": 60})
-    save(maps, "poisson-convergence-maps", show=show)
+    save_figure(maps, "poisson-convergence-maps", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/poisson-source.py b/docs/src/examples/poisson-source.py
index a93c3ad..7042b61 100644
--- a/docs/src/examples/poisson-source.py
+++ b/docs/src/examples/poisson-source.py
@@ -16,6 +16,7 @@
 from struphy import (EnvironmentOptions, Simulation, Time, domains, grids,
                      perturbations)
 from struphy.models import Poisson
+from struphy_plots import save_figure
 
 # Problem parameters: a single cosine mode in space, oscillating in time.
 WAVENUMBER = 2
@@ -61,27 +62,6 @@ def create_simulation() -> Simulation:
     return simulation
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output.pproc(create_vtk=False)
     domain = sim.domain
@@ -113,7 +93,7 @@ def phi_exact(x, t):
         backend="plotly",
     )
 
-    save(figure, "poisson-source", show=show)
+    save_figure(figure, "poisson-source", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/pressureless-transport.py b/docs/src/examples/pressureless-transport.py
index 8b7dab7..1e05e58 100644
--- a/docs/src/examples/pressureless-transport.py
+++ b/docs/src/examples/pressureless-transport.py
@@ -18,6 +18,7 @@
     domains, equils, grids, perturbations,
 )
 from struphy.models import VariationalPressurelessFluid
+from struphy_plots import save_figure
 
 stem = "pressureless-transport"
 length, speed, amplitude = 2.0 * np.pi, 0.5, 0.2
@@ -56,27 +57,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     from plotly.subplots import make_subplots
@@ -122,12 +102,12 @@ def pproc(sim: Simulation, show: bool = False):
     figure.update_yaxes(title_text="absolute relative error", row=2, col=1)
     figure.update_layout(title="Pressureless transport at constant velocity", template="plotly_white",
                          legend={"orientation": "h", "y": -0.18}, margin={"l": 85, "r": 30, "t": 90, "b": 140})
-    save(figure, stem, height=800, show=show)
+    save_figure(figure, stem, height=800, show=show)
     movie = rho.assign_coords(eta1=x).struphy.plot.slice(
         x="eta1", y="t", symmetric=True, cmap="RdBu_r", title="Density transported around a periodic box",
         xlabel="x", ylabel="t [a.u.]", colorbar_label="ρ", backend="plotly",
     )
-    save(movie, f"{stem}-space-time", show=show)
+    save_figure(movie, f"{stem}-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/resistive-diffusion.py b/docs/src/examples/resistive-diffusion.py
index a46c515..9a0fb2d 100644
--- a/docs/src/examples/resistive-diffusion.py
+++ b/docs/src/examples/resistive-diffusion.py
@@ -29,6 +29,7 @@
 )
 from struphy.linear_algebra.solver import NonlinearSolverParameters
 from struphy.models import ViscoResistiveMHD
+from struphy_plots import save_figure
 
 length = 2 * np.pi
 mode_number = 2
@@ -98,27 +99,6 @@ def mode_amplitude(run):
     return b_z.struphy.analysis.project_mode(dim="eta1", number=mode_number)
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -182,7 +162,7 @@ def profile_traces(index):
     figure.update_yaxes(title_text="B_z", range=[-1.2 * amplitude, 1.2 * amplitude], row=1, col=1)
     figure.update_xaxes(title_text="t", row=2, col=1)
     figure.update_yaxes(title_text="energy change / magnetic energy lost", row=2, col=1)
-    save(figure, "resistive-diffusion", width=900, height=850, show=show)
+    save_figure(figure, "resistive-diffusion", width=900, height=850, show=show)
 
     first, *others = amplitudes.values()
     decay = first.struphy.plot.timeseries(
@@ -193,7 +173,7 @@ def profile_traces(index):
         backend="plotly",
     )
     decay.fig.update_yaxes(title_text="amplitude of the sin(kx) mode")
-    save(decay, "resistive-diffusion-decay", show=show)
+    save_figure(decay, "resistive-diffusion-decay", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/resistive-x-point.py b/docs/src/examples/resistive-x-point.py
index 5304574..eba384b 100644
--- a/docs/src/examples/resistive-x-point.py
+++ b/docs/src/examples/resistive-x-point.py
@@ -28,6 +28,7 @@
 )
 from struphy.linear_algebra.solver import NonlinearSolverParameters
 from struphy.models import ViscoResistiveMHD
+from struphy_plots import save_figure
 
 period = 2 * np.pi
 box_min, box_max = -np.pi, np.pi
@@ -115,27 +116,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     time_opts = sim.time_opts
     from plotly.subplots import make_subplots
@@ -241,7 +221,7 @@ def field_traces(index):
     # The still shows every layer at the final time, the flux contours included (the frames leave them static).
     still = go.Figure(data=field_traces(len(times) - 1), layout=figure.layout)
     still.layout.sliders[0].active = len(figure.frames) - 1
-    save(figure, "resistive-x-point", width=900, height=850, still=still, show=show)
+    save_figure(figure, "resistive-x-point", width=900, height=850, still=still, show=show)
 
     reconnection = make_subplots(
         rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.16,
@@ -286,8 +266,8 @@ def field_traces(index):
         margin={"l": 80, "r": 35, "t": 90, "b": 75}, legend={"orientation": "h", "y": -0.18},
     )
 
-    save(reconnection, "resistive-x-point-reconnection", show=show)
-    save(conservation, "resistive-x-point-conservation", show=show)
+    save_figure(reconnection, "resistive-x-point-reconnection", show=show)
+    save_figure(conservation, "resistive-x-point-conservation", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/shear-alfven-wave.py b/docs/src/examples/shear-alfven-wave.py
index a469d3b..1d6913e 100644
--- a/docs/src/examples/shear-alfven-wave.py
+++ b/docs/src/examples/shear-alfven-wave.py
@@ -27,6 +27,7 @@
     perturbations,
 )
 from struphy.models import ShearAlfven
+from struphy_plots import save_figure
 
 
 def create_simulation() -> Simulation:
@@ -77,27 +78,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     domain = sim.domain
 
@@ -128,7 +108,7 @@ def pproc(sim: Simulation, show: bool = False):
         backend="plotly",
     )
 
-    save(figure, "shear-alfven-wave", show=show)
+    save_figure(figure, "shear-alfven-wave", show=show)
 
     # The same field along z, over time: waves travelling in both directions leave diagonal
     # stripes, whose slope is the wave speed.
@@ -147,7 +127,7 @@ def pproc(sim: Simulation, show: bool = False):
         colorbar_label="u₁ (logical component)",
         backend="plotly",
     )
-    save(space_time, "shear-alfven-wave-space-time", show=show)
+    save_figure(space_time, "shear-alfven-wave-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/sph-velocity-diffusion.py b/docs/src/examples/sph-velocity-diffusion.py
index c89de6c..a155a39 100644
--- a/docs/src/examples/sph-velocity-diffusion.py
+++ b/docs/src/examples/sph-velocity-diffusion.py
@@ -30,6 +30,7 @@
 )
 from struphy.models import ViscousEulerSPH
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 length = 1.0
 viscosity = 0.05
@@ -86,27 +87,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     output.pproc()
@@ -135,7 +115,7 @@ def pproc(sim: Simulation, show: bool = False):
         title="SPH velocity diffusion: sinusoidal mode against its exact decay",
         backend="plotly",
     )
-    save(figure, "sph-velocity-diffusion", show=show)
+    save_figure(figure, "sph-velocity-diffusion", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/strong-landau-damping.py b/docs/src/examples/strong-landau-damping.py
index 39793f7..e231674 100644
--- a/docs/src/examples/strong-landau-damping.py
+++ b/docs/src/examples/strong-landau-damping.py
@@ -33,6 +33,7 @@
     perturbations,
 )
 from struphy.models import VlasovAmpereOneSpecies
+from struphy_plots import save_figure
 
 # A large-amplitude cosine mode, well past the linear (weak-damping) regime.
 perturbation_amplitude = 0.5
@@ -89,27 +90,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     output = sim.output
     field_energy_array = output.scalars["electric_energy"]
@@ -125,7 +105,7 @@ def pproc(sim: Simulation, show: bool = False):
         logy=True, title="Strong Landau damping: electric field energy", backend="plotly",
     )
 
-    save(figure, "strong-landau-damping", show=show)
+    save_figure(figure, "strong-landau-damping", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/toroidal-shear-alfven.py b/docs/src/examples/toroidal-shear-alfven.py
index a2b059e..855f242 100644
--- a/docs/src/examples/toroidal-shear-alfven.py
+++ b/docs/src/examples/toroidal-shear-alfven.py
@@ -31,6 +31,7 @@
     perturbations,
 )
 from struphy.models import LinearMHD
+from struphy_plots import save_figure
 
 STEM = "toroidal-shear-alfven"
 NUM_ELEMENTS = (8, 48, 4)
@@ -139,7 +140,7 @@ def save_fixed_theta_fft_figures(output, show=False):
             margin={"l": 85, "r": 35, "t": 95, "b": 75},
             legend={"title": {"text": "Poloidal angle"}},
         )
-        save(figure, f"{STEM}-{key}", show=show)
+        save_figure(figure, f"{STEM}-{key}", show=show)
 
 
 def create_simulation() -> Simulation:
@@ -194,27 +195,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     """Save (and show) the figures from an existing simulation."""
     domain = sim.domain
@@ -312,7 +292,7 @@ def traces(index):
     )
     # A later snapshot also shows the generated toroidal component in the thumbnail.
     still = len(frame_indices) // 2
-    save(figure, STEM, frame=still, show=show)
+    save_figure(figure, STEM, frame=still, show=show)
 
     radial_probe = int(np.abs(velocity.eta1.values - 0.5).argmin())
     probe_radius = domain.params["a1"] + (domain.params["a2"] - domain.params["a1"]) * float(velocity.eta1.values[radial_probe])
@@ -519,7 +499,7 @@ def traces(index):
             template="plotly_white", margin={"l": 85, "r": 35, "t": 95, "b": 75},
             legend={"title": {"text": "Poloidal harmonic"}},
         )
-        save(radial_mode_plot, f"{STEM}-{key}", show=show)
+        save_figure(radial_mode_plot, f"{STEM}-{key}", show=show)
 
     energy = go.Figure()
     for key, label in (("en_U", "Kinetic"), ("en_B", "Magnetic"), ("en_thermal", "Compressional")):
@@ -531,16 +511,16 @@ def traces(index):
         title="Perturbation energies", xaxis_title="t", yaxis_title="Energy [normalized units]",
         template="plotly_white", margin={"l": 80, "r": 30, "t": 80, "b": 65},
     )
-    save(history, f"{STEM}-velocity-history", show=show)
+    save_figure(history, f"{STEM}-velocity-history", show=show)
     for angle_degrees, radial_profiles in radial_profile_figures:
-        save(radial_profiles, f"{STEM}-radial-profiles-theta-{int(angle_degrees)}", show=show)
-    save(radial_history, f"{STEM}-radial-history", show=show)
-    save(mode_plot, f"{STEM}-poloidal-fft", show=show)
+        save_figure(radial_profiles, f"{STEM}-radial-profiles-theta-{int(angle_degrees)}", show=show)
+    save_figure(radial_history, f"{STEM}-radial-history", show=show)
+    save_figure(mode_plot, f"{STEM}-poloidal-fft", show=show)
     save_fixed_theta_fft_figures(output, show=show)
-    save(frequency_plot, f"{STEM}-time-fft", show=show)
-    save(radius_spectrum, f"{STEM}-radial-time-fft", show=show)
-    save(filtered_probe, f"{STEM}-filtered-velocity", show=show)
-    save(energy, f"{STEM}-energy", show=show)
+    save_figure(frequency_plot, f"{STEM}-time-fft", show=show)
+    save_figure(radius_spectrum, f"{STEM}-radial-time-fft", show=show)
+    save_figure(filtered_probe, f"{STEM}-filtered-velocity", show=show)
+    save_figure(energy, f"{STEM}-energy", show=show)
     dominant = ", ".join(f"{label} {value:.3f}" for label, value in zip(labels, band.spectrum.dominant_frequency.values))
     print(f"Dominant temporal frequencies: {dominant} (resolution {frequency_resolution:.3f})")
 
diff --git a/docs/src/examples/two-stream-instability.py b/docs/src/examples/two-stream-instability.py
index 59b1a00..6b94291 100644
--- a/docs/src/examples/two-stream-instability.py
+++ b/docs/src/examples/two-stream-instability.py
@@ -33,6 +33,7 @@
     perturbations,
 )
 from struphy.models import VlasovAmpereOneSpecies
+from struphy_plots import save_figure
 
 # Two counter-streaming Maxwellians (u1 = +/-3), each seeded with the same cosine mode.
 perturbation_amplitude = 0.001
@@ -96,27 +97,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     from struphy_plots.theory.kinetic import two_stream
 
@@ -135,7 +115,7 @@ def pproc(sim: Simulation, show: bool = False):
 
     figure = field_energy.struphy.plot.timeseries(logy=True, title="Two-stream instability: electric field energy", backend="plotly")
 
-    save(figure, "two-stream-instability", show=show)
+    save_figure(figure, "two-stream-instability", show=show)
 
     output.pproc()
     f = output.evaluate("kinetic_ions/e1_v1_density/f")  # (t, e1, v1)
@@ -170,8 +150,8 @@ def pproc(sim: Simulation, show: bool = False):
         backend="plotly",
     )
 
-    save(velocity_time, "two-stream-instability-velocity-time", show=show)
-    save(phase_space, "two-stream-instability-phasespace", frame=len(phase_space.fig.frames) // 2, show=show)
+    save_figure(velocity_time, "two-stream-instability-velocity-time", show=show)
+    save_figure(phase_space, "two-stream-instability-phasespace", frame=len(phase_space.fig.frames) // 2, show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/vlasov-tokamak.py b/docs/src/examples/vlasov-tokamak.py
index 347c95c..1503c05 100644
--- a/docs/src/examples/vlasov-tokamak.py
+++ b/docs/src/examples/vlasov-tokamak.py
@@ -29,6 +29,7 @@
 from struphy.kinetic_background import maxwellians
 from struphy.models import Vlasov
 from struphy.pic.base import BoundaryParameters, LoadingParameters, SavingParameters
+from struphy_plots import save_figure
 
 # Load a small population; only a handful of markers have their full orbit saved.
 # eta1 is the radial flux coordinate on this domain, so it must reflect at the
@@ -84,27 +85,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     domain = sim.domain
     time_opts = sim.time_opts
@@ -247,7 +227,7 @@ def comet_arrays(p: int, idx: int):
         ],
     )
 
-    save(figure, "vlasov-tokamak", height=850, show=show)
+    save_figure(figure, "vlasov-tokamak", height=850, show=show)
 
     # The same orbits, projected onto each coordinate plane -- a simpler,
     # non-animated companion to the 3D view above, useful for reading off
@@ -281,7 +261,7 @@ def comet_arrays(p: int, idx: int):
         margin={"l": 60, "r": 30, "t": 80, "b": 60},
     )
 
-    save(projection_figure, "vlasov-tokamak-projections", width=1500, height=560, show=show)
+    save_figure(projection_figure, "vlasov-tokamak-projections", width=1500, height=560, show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/vortex-merger.py b/docs/src/examples/vortex-merger.py
index 3210bd3..40b39b6 100644
--- a/docs/src/examples/vortex-merger.py
+++ b/docs/src/examples/vortex-merger.py
@@ -35,6 +35,7 @@
     maxwellians,
 )
 from struphy.models import ToyDrift
+from struphy_plots import save_figure
 
 # Two Gaussian blobs of peak density `peak` and width `width`, on the circle of radius `ring_radius`,
 # with their centres `separation` apart.
@@ -108,27 +109,6 @@ def two_blobs(etas):
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     a1, a2 = sim.domain.params["a1"], sim.domain.params["a2"]
     time_opts = sim.time_opts
@@ -165,7 +145,7 @@ def pproc(sim: Simulation, show: bool = False):
         title="Vortex merger: binned charge density", xlabel="x", ylabel="y", colorbar_label="density",
         equal_aspect=True, backend="plotly",
     )
-    save(movie, "vortex-merger", height=750, frame=len(movie.fig.frames) // 2, show=show)
+    save_figure(movie, "vortex-merger", height=750, frame=len(movie.fig.frames) // 2, show=show)
 
     # The first Poisson solve initializes the field energy after t=0. Compare
     # subsequent field energies to that first solved state, not to the zero placeholder.
@@ -177,7 +157,7 @@ def pproc(sim: Simulation, show: bool = False):
     )
     drift = float(energy.struphy.analysis.relative_error().max())
     print(f"Maximum relative drift of the electrostatic energy: {drift:.2e}")
-    save(energy_figure, "vortex-merger-energy", show=show)
+    save_figure(energy_figure, "vortex-merger-energy", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/weak-landau-damping.py b/docs/src/examples/weak-landau-damping.py
index 85d09f5..9ad0512 100644
--- a/docs/src/examples/weak-landau-damping.py
+++ b/docs/src/examples/weak-landau-damping.py
@@ -32,6 +32,7 @@
     perturbations,
 )
 from struphy.models import VlasovAmpereOneSpecies
+from struphy_plots import save_figure
 
 # A single small-amplitude cosine mode perturbs an otherwise uniform Maxwellian.
 perturbation_amplitude = 0.001
@@ -89,25 +90,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
 def pproc(sim: Simulation, show: bool = False):
     domain = sim.domain
 
@@ -135,7 +117,7 @@ def field_energy_exact(t):
         backend="plotly",
     )
 
-    save(figure, "weak-landau-damping", show=show)
+    save_figure(figure, "weak-landau-damping", show=show)
 
     # The electric field along x, over time.
     output.pproc()
@@ -154,7 +136,7 @@ def field_energy_exact(t):
         colorbar_label="E_x",
         backend="plotly",
     )
-    save(space_time, "weak-landau-damping-space-time", show=show)
+    save_figure(space_time, "weak-landau-damping-space-time", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/weibel-instability.py b/docs/src/examples/weibel-instability.py
index 4e37d60..e56bcb5 100644
--- a/docs/src/examples/weibel-instability.py
+++ b/docs/src/examples/weibel-instability.py
@@ -35,6 +35,7 @@
     perturbations,
 )
 from struphy.models import VlasovMaxwellOneSpecies
+from struphy_plots import save_figure
 
 wavenumber = 1.25
 
@@ -105,27 +106,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
-
-
 def pproc(sim: Simulation, show: bool = False):
     domain = sim.domain
 
@@ -151,7 +131,7 @@ def pproc(sim: Simulation, show: bool = False):
         backend="plotly",
     )
 
-    save(figure, "weibel-instability", show=show)
+    save_figure(figure, "weibel-instability", show=show)
 
     # The magnetic field along x, over time.
     magnetic_field = b_field.isel(component=2, eta2=0, eta3=0)  # (t, e1)
@@ -179,8 +159,8 @@ def pproc(sim: Simulation, show: bool = False):
         backend="plotly",
     )
 
-    save(space_time, "weibel-instability-space-time", show=show)
-    save(anisotropy_figure, "weibel-instability-anisotropy", show=show)
+    save_figure(space_time, "weibel-instability-space-time", show=show)
+    save_figure(anisotropy_figure, "weibel-instability-anisotropy", show=show)
 
 
 if __name__ == "__main__":
diff --git a/docs/src/examples/zeldovich-caustic.py b/docs/src/examples/zeldovich-caustic.py
index ba8e13f..e08be0d 100644
--- a/docs/src/examples/zeldovich-caustic.py
+++ b/docs/src/examples/zeldovich-caustic.py
@@ -38,6 +38,7 @@
 )
 from struphy.models import PressureLessSPH
 from struphy.ode.utils import ButcherTableau
+from struphy_plots import save_figure
 
 amplitude = 0.5  # u(x) = amplitude * sin(2 pi x) on the periodic unit interval
 caustic_time = 1.0 / (2 * np.pi * amplitude)  # the density is singular here, at x = 1/2
@@ -105,25 +106,6 @@ def create_simulation() -> Simulation:
     return sim
 
 
-def save(figure, name: str, *, show: bool = False, frame: int | None = None, still=None, width=1100, height=650):
-    """Save a Plotly figure as ``<name>.html``, ``<name>.png`` and ``<name>.plotly.json``.
-
-    ``figure`` is a plot of struphy-plots drawn with ``backend="plotly"``, or a
-    ``plotly.graph_objects.Figure``. ``show`` shows it first. For an animation, the PNG shows
-    ``frame`` (default: the first), or the figure ``still`` instead. Under MPI only rank 0 writes.
-    """
-    import struphy_plots
-    from struphy_plots.plotting import PlotResult
-
-    if not struphy_plots.is_plotting_rank():
-        return
-    result = figure if isinstance(figure, PlotResult) else PlotResult(figure, None)
-    if show:
-        result.show()
-    result.save(f"{name}.html")
-    image = PlotResult(still, None) if still is not None else result
-    image.save(f"{name}.png", frame=frame, width=width, height=height, scale=2)
-    result.save(f"{name}.plotly.json")
 def pproc(sim: Simulation, show: bool = False):
     from plotly.subplots import make_subplots
 
@@ -184,7 +166,7 @@ def traces(index):
     figure.update_yaxes(title_text="ρ", range=[0, density_ceiling], row=1, col=1)
     figure.update_yaxes(title_text="u", range=[-1.1 * amplitude, 1.1 * amplitude], row=2, col=1)
     late = int(np.argmin(np.abs(times - 1.5 * caustic_time)))
-    save(figure, "zeldovich-caustic", width=900, height=850, frame=int(np.argmin(np.abs(picks - late))), show=show)
+    save_figure(figure, "zeldovich-caustic", width=900, height=850, frame=int(np.argmin(np.abs(picks - late))), show=show)
 
 
 if __name__ == "__main__":

From 9e8ca1e417278c7865f49f3ddd9f0dc18b98497e Mon Sep 17 00:00:00 2001
From: Max <max.lindqvist@ipp.mpg.de>
Date: Mon, 28 Sep 2026 08:47:51 +0200
Subject: [PATCH 4/4] Update submodule

---
 submodules/struphy-plots | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/submodules/struphy-plots b/submodules/struphy-plots
index e6dfc8d..a5e1a65 160000
--- a/submodules/struphy-plots
+++ b/submodules/struphy-plots
@@ -1 +1 @@
-Subproject commit e6dfc8d1920c45f2c7b39903eb51d6db6cf4a81c
+Subproject commit a5e1a65620001160ba5bdf743506746c0112af7b