diff --git a/.gitignore b/.gitignore
index 9b1a5785..00214deb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,5 @@ _git2_*
 docs/build
 *.info
 .vscode/spellright.dict
+.CondaPkg
+deps/build.log
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 0cd8a104..7ef9bcbc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,11 @@ uuid = "30eb0fb0-5147-11e9-3356-d75b018717ce"
 authors = ["chrisbrahms <38351086+chrisbrahms@users.noreply.github.com>"]
 version = "0.5.0"
 
+[extensions]
+PythonPlotExt = ["PythonPlot", "PythonCall"]
+PyPlotExt = ["PyPlot"]
+GLMakieExt = ["GLMakie"]
+
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
@@ -39,8 +44,6 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
-PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
-PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -52,6 +55,12 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
+[weakdeps]
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+
 [compat]
 ArgParse = "1"
 BlackBoxOptim = "0.6"
@@ -88,8 +97,6 @@ Pkg = "1.9"
 Polynomials = "2, 3, 4"
 Printf = "1.9"
 ProgressLogging = "0.1"
-PyCall = "1.92"
-PyPlot = "2.9"
 QuadGK = "2.4"
 Random = "1.9"
 Reexport = "1.2"
diff --git a/README.md b/README.md
index c2433bc2..1e367d66 100644
--- a/README.md
+++ b/README.md
@@ -36,13 +36,18 @@ There are two ways of using Luna:
 For a short introduction on how to use the simple interface, see the [Quickstart](#quickstart) or [GNLSE](#gnlse) sections below. More information, including on the internals of Luna, can be found in the [Documentation](http://lupo-lab.com/Luna.jl).
 
 ## Installation
-Luna requires Julia v1.9 or later, which can be obtained from [here](https://julialang.org/downloads/). In a Julia terminal, to install Luna simply enter the package manager with `]` and run `add Luna`:
+Luna requires Julia v1.9 or later, which can be obtained from [here](https://julialang.org/downloads/). It is recommended to always use the most recent released Julia version. In a Julia terminal, to install Luna simply enter the package manager with `]` and run `add Luna`:
 
 ```julia
 ]
 add Luna
 ```
-This will install and precompile Luna and all its dependencies.
+This will install and precompile Luna and all its dependencies. For plotting you also need to install a plotting backend. We recommend `PythonPlot` which can be installed by calling:
+
+```julia
+]
+add PythonPlot, PythonCall
+```
 
 ## Quickstart
 To run a simple simulation of ultrafast pulse propagation in a gas-filled hollow capillary fibre, you can use `prop_capillary`. As an example, take a 3-metre length of HCF with 125 μm core radius, filled with 1 bar of helium gas, and driving pulses centred at 800 nm wavelength with 120 μJ of energy and 10 fs duration. We consider a frequency grid which spans from 120 nm to 4 μm and a time window of 1 ps.
@@ -74,7 +79,9 @@ The shape of this array is `(Nω x Nz)` where `Nω` is the number of frequency s
 
 Mode-averaged propagation is activated using `modes=:HE11` (the default) or replacing the `:HE11` with a different mode designation (for mode-averaged propagation in a different mode). To run the same simulation as above with the first four modes (HE₁₁ to HE₁₄) of the capillary, set `modes` to `4` (this example also uses smaller time and frequency windows to make the simulation run a little faster):
 ```julia
-julia> prop_capillary(125e-6, 3, :He, 1; λ0=800e-9, modes=4, energy=120e-6, τfwhm=10e-15, trange=400e-15, λlims=(150e-9, 4e-6))
+julia> output_multimode = prop_capillary(125e-6, 3, :He, 1; λ0=800e-9, modes=4, energy=120e-6, τfwhm=10e-15, trange=400e-15, λlims=(150e-9, 4e-6))
+[...]
+MemoryOutput["simulation_type", "dumps", "meta", "Eω", "prop_capillary_args", "grid", "stats", "z"]
 ```
 The propagation will take much longer, and the output field `Eω` now has shape `(Nω x Nm x Nz)` with `Nm` the number of modes:
 ```julia
@@ -84,31 +91,32 @@ julia> output_multimode["Eω"]
 ```
 **NOTE:** Setting `modes=:HE11` and `modes=1` are **not** equivalent, except if only the Kerr effect is included in the simulation. The former uses mode-averaged propagation (treating all spatial dependence of the nonlinear polarisation the same as the Kerr effect) whereas the latter projects the spatially dependent nonlinear polarisation onto a single mode. This difference is especially important when photoionisation plays a major role.
 ### Plotting results
-More usefully, you can directly plot the propagation results using `Plotting.prop_2D()` (`Plotting` is imported at the same time as `prop_capillary` by the `using Luna` statement):
+More usefully, you can directly plot the propagation results using `Plotting.prop_2D()`. While `Plotting` is imported at the same time as `prop_capillary` by the `using Luna` statement, you also need to import a plotting backend, e.g. `PythonPlot`:
 ```julia
+julia> using PythonPlot
 julia> Plotting.prop_2D(output)
-PyPlot.Figure(PyObject <Figure size 2400x800 with 4 Axes>)
+Python Figure: <Figure size 1200x400 with 4 Axes>
 ```
 This should show a plot like this:
 ![Propagation example 1](assets/readme_modeAvgProp.png)
 You can also display the power spectrum at the input and output (and anywhere in between):
 ```julia
 julia> Plotting.spec_1D(output, [0, 1.5, 3]; log10=true)
-PyPlot.Figure(PyObject <Figure size 1700x1000 with 1 Axes>)
+Python Figure: <Figure size 850x500 with 1 Axes>
 ```
 which will show this:
 ![Propagation example 2](assets/readme_modeAvgSpec.png)
 `Plotting` functions accept many additional keyword arguments to quickly display relevant information. For example, you can show the bandpass-filtered UV pulse from the simulation using the `bandpass` argument:
 ```julia
 julia> Plotting.time_1D(output, [2, 2.5, 3]; trange=(-10e-15, 30e-15), bandpass=(180e-9, 220e-9))
-PyPlot.Figure(PyObject <Figure size 1700x1000 with 1 Axes>)
+Python Figure: <Figure size 850x500 with 1 Axes>
 ```
 ![Propagation example 3](assets/readme_modeAvgTime.png)
 
 For multi-mode simulations, the plotting functions will display all modes individually by default. You can display the sum over modes instead using `modes=:sum`:
 ```julia
 julia> Plotting.spec_1D(output_multimode; log10=true, modes=:sum)
-PyPlot.Figure(PyObject <Figure size 1700x1000 with 1 Axes>)
+Python Figure: <Figure size 850x500 with 1 Axes>
 ```
 ![Propagation example 4](assets/readme_multiModeSpec.png)
 (Compare this to the mode-averaged case above and note the important differences, e.g. the appearance of additional ultraviolet dispersive waves in higher-order modes.)
@@ -130,19 +138,22 @@ julia> using Luna
 julia> γ = 0.11
 julia> flength = 15e-2
 julia> βs = [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
-julia> output = prop_gnlse(γ, flength, βs; λ0=835e-9, τfwhm=50e-15, power=10e3, pulseshape=:sech, λlims=(400e-9, 2400e-9), trange=12.5e-12)
+julia> output_gnlse = prop_gnlse(γ, flength, βs; λ0=835e-9, τfwhm=50e-15, power=10e3, pulseshape=:sech, λlims=(400e-9, 2400e-9), trange=12.5e-12)
+[...]
+MemoryOutput["simulation_type", "dumps", "meta", "Eω", "prop_capillary_args", "grid", "stats", "z"]
 ```
 After this has run, you can visualise the output, with e.g.
 ```julia
-julia> Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
-PyPlot.Figure(PyObject <Figure size 2400x800 with 4 Axes>)
+julia> using PythonPlot
+julia> Plotting.prop_2D(output_gnlse, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+Python Figure: <Figure size 1200x400 with 4 Axes>
 ```
 This should show a plot like this:
 ![GNLSE propagation example](assets/readme_gnlse_scg.png)
 
 
 ## Examples
-The [examples folder](examples/) contains complete simulation examples for a variety of scenarios, both for the [simple interface](examples/simple_interface/) and the [low-level interface](examples/low_level_interface). Some of the simple interface examples require the `PyPlot` package to be present, and many of the low-level examples require other packages as well--you can install these by simply typing `] add PyPlot` at the Julia REPL or the equivalent for other packages.
+The [examples folder](examples/) contains complete simulation examples for a variety of scenarios, both for the [simple interface](examples/simple_interface/) and the [low-level interface](examples/low_level_interface). Some of the simple interface examples require the `PythonPlot` package to be present, and many of the low-level examples require other packages as well--you can install these by simply typing `] add PythonPlot` at the Julia REPL or the equivalent for other packages.
 
 ## The low-level interface
 At its core, Luna is extremely flexible, and the simple interface using `prop_capillary` only exposes part of what Luna can do. There are lots of examples in the [low-level interface examples folder](examples/low_level_interface). These are not actively maintained and are not guaranteed to run. As a side effect of its flexibility, it is quite easy to make mistakes when using the low-level interface. For example, changing from single-mode to multi-mode propagation in a fibre requires several concurrent changes to your code. If you have trouble with this interface, [open an issue](https://github.com/LupoLab/Luna/issues/new) with as much detail as possible and we will try to help you run it.
@@ -150,6 +161,9 @@ At its core, Luna is extremely flexible, and the simple interface using `prop_ca
 ## Running parameter scans
 Luna comes with a built-in interface which allows for the running of single- and multi-dimensional parameter scans with very little additional code. An example can be found in the [examples folder](examples/simple_interface/scan.jl) and more information is available in the [documentation](http://lupo-lab.com/Luna.jl/dev/scans.html).
 
+## Plotting backends
+Luna supports three plotting backends; PythonPlot, PyPlot and GLMakie. For most use cases, you should use PythonPlot. We support PyPlot for backwards compatibility. Currently, support for GLMakie is experimental.
+
 ## New to Julia?
 There are many resources to help you learn Julia. A good place to start is [Julia Academy](https://juliaacademy.com/) which has several courses for learning Julia depending on your current experience. There are additional resources linked from the [Julia website](https://julialang.org/learning/).
 
diff --git a/docs/src/scans.md b/docs/src/scans.md
index 0b7d6ff2..243a3094 100644
--- a/docs/src/scans.md
+++ b/docs/src/scans.md
@@ -4,7 +4,7 @@ Luna comes with a flexible interface to run, save and process scans over any par
 **First**, define the fixed parameters (those which are not being scanned over):
 ```julia
 using Luna
-import PyPlot: plt
+import PythonPlot: pyplot
 
 a = 125e-6
 flength = 3
@@ -99,7 +99,7 @@ julia> size(λ)
 ```
 With this data, we can now plot the energy-pressure scan:
 ```julia
-fig, axs = plt.subplots(1, length(pressures))
+fig, axs = pyplot.subplots(1, length(pressures))
 fig.set_size_inches(8, 2)
 for (pidx, pressure) in enumerate(pressures)
     ax = axs[pidx]
@@ -110,7 +110,7 @@ for (pidx, pressure) in enumerate(pressures)
     ax.set_title("Pressure: $pressure bar")
     ax.set_xlim(100, 1200)
 end
-plt.colorbar(img, ax=axs, label="Energy density (dB)")
+pyplot.colorbar(img, ax=axs, label="Energy density (dB)")
 ```
 which will produce this figure:
 ![Scan spectra](assets/scan_spectrum.png)
diff --git a/examples/low_level_interface/MI_modeAvg_env.jl b/examples/low_level_interface/MI_modeAvg_env.jl
index 35dd84b2..5ae18643 100644
--- a/examples/low_level_interface/MI_modeAvg_env.jl
+++ b/examples/low_level_interface/MI_modeAvg_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 15e-6
 gas = :Ar
@@ -40,7 +40,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output; trange=(-1e-12, 1e-12), λrange=(220e-9, 2000e-9))
 Plotting.time_1D(output; trange=(-1e-12, 1e-12))
diff --git a/examples/low_level_interface/PPT_ionisation_rate/PPT_ionisation_rate.jl b/examples/low_level_interface/PPT_ionisation_rate/PPT_ionisation_rate.jl
index 3e46d6e3..043fc322 100644
--- a/examples/low_level_interface/PPT_ionisation_rate/PPT_ionisation_rate.jl
+++ b/examples/low_level_interface/PPT_ionisation_rate/PPT_ionisation_rate.jl
@@ -1,6 +1,5 @@
-using Luna
+using Luna, PythonPlot
 import DelimitedFiles: readdlm
-import PyPlot: plt, pygui
 import Printf: @sprintf
 import GSL: hypergeom
 import SpecialFunctions: dawson, gamma
diff --git a/examples/low_level_interface/Raman/Raman_modeAvg.jl b/examples/low_level_interface/Raman/Raman_modeAvg.jl
index c6cdd9c7..ed854ee9 100644
--- a/examples/low_level_interface/Raman/Raman_modeAvg.jl
+++ b/examples/low_level_interface/Raman/Raman_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :H2
@@ -35,7 +35,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, range(0,flength,length=5), trange=(-0.5e-12, 2e-12))
diff --git a/examples/low_level_interface/Raman/Raman_modeAvg_env.jl b/examples/low_level_interface/Raman/Raman_modeAvg_env.jl
index 88ca4fe3..257849c4 100644
--- a/examples/low_level_interface/Raman/Raman_modeAvg_env.jl
+++ b/examples/low_level_interface/Raman/Raman_modeAvg_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :H2
@@ -33,7 +33,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, range(0,flength,length=5), trange=(-0.5e-12, 2e-12))
diff --git a/examples/low_level_interface/Raman/Raman_modeAvg_noTHG.jl b/examples/low_level_interface/Raman/Raman_modeAvg_noTHG.jl
index fd9ec5ac..73367502 100644
--- a/examples/low_level_interface/Raman/Raman_modeAvg_noTHG.jl
+++ b/examples/low_level_interface/Raman/Raman_modeAvg_noTHG.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :H2
@@ -34,7 +34,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, range(0,flength,length=5), trange=(-0.5e-12, 2e-12))
diff --git a/examples/low_level_interface/SFA_gauss.jl b/examples/low_level_interface/SFA_gauss.jl
index 6efefcd2..777eb102 100644
--- a/examples/low_level_interface/SFA_gauss.jl
+++ b/examples/low_level_interface/SFA_gauss.jl
@@ -1,8 +1,7 @@
 import FFTW
 import Luna: SFA, Tools, Grid, Maths, PhysData, Ionisation
 import Luna.PhysData: wlfreq
-import PyPlot: plt, pygui
-pygui(true)
+import PythonPlot: pyplot
 
 τfwhm = 30e-15
 λ0 = 800e-9
@@ -33,23 +32,23 @@ eV0 = PhysData.ħ*wlfreq(λ0)/PhysData.electron
 harmonics = eV./eV0
 
 ##
-plt.figure()
-plt.semilogy(harmonics, abs2.(Dω), label="Harmonic spectrum")
-plt.axvline(cutoffeV/eV0, linestyle="--", color="k", label="IP + 3.17 Up")
-plt.xlim(0, 1.2cutoffeV/eV0)
-plt.xlabel("Harmonic order")
-plt.ylabel("Spectral energy density")
-plt.legend(frameon=false)
+pyplot.figure()
+pyplot.semilogy(harmonics, abs2.(Dω), label="Harmonic spectrum")
+pyplot.axvline(cutoffeV/eV0, linestyle="--", color="k", label="IP + 3.17 Up")
+pyplot.xlim(0, 1.2cutoffeV/eV0)
+pyplot.xlabel("Harmonic order")
+pyplot.ylabel("Spectral energy density")
+pyplot.legend(frameon=false)
 
 ##
-plt.figure()
-plt.pcolormesh(tg*1e15, eV[1:4:end], Maths.log10_norm(abs2.(gab[1:4:end, :])))
-plt.axhline(cutoffeV, linestyle="--", color="w")
-plt.clim(-10, -4)
-plt.ylim(20, 1.2*cutoffeV)
-plt.xlim(-τfwhm*1e15, τfwhm*1e15)
-cb = plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(tg*1e15, eV[1:4:end], Maths.log10_norm(abs2.(gab[1:4:end, :])))
+pyplot.axhline(cutoffeV, linestyle="--", color="w")
+pyplot.clim(-10, -4)
+pyplot.ylim(20, 1.2*cutoffeV)
+pyplot.xlim(-τfwhm*1e15, τfwhm*1e15)
+cb = pyplot.colorbar()
 cb.set_label("Log10 SED")
-plt.xlabel("Time (fs)")
-plt.ylabel("Photon energy (eV)")
+pyplot.xlabel("Time (fs)")
+pyplot.ylabel("Photon energy (eV)")
 
diff --git a/examples/low_level_interface/arpcf_modeAvg.jl b/examples/low_level_interface/arpcf_modeAvg.jl
index e2b8fde6..11dcff5f 100644
--- a/examples/low_level_interface/arpcf_modeAvg.jl
+++ b/examples/low_level_interface/arpcf_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -42,7 +42,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output)
diff --git a/examples/low_level_interface/basic_modal.jl b/examples/low_level_interface/basic_modal.jl
index bd892e5d..11fbec35 100644
--- a/examples/low_level_interface/basic_modal.jl
+++ b/examples/low_level_interface/basic_modal.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -40,7 +40,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.prop_2D(output; bandpass=(150e-9, 300e-9), modes=:sum)
diff --git a/examples/low_level_interface/basic_modal_env.jl b/examples/low_level_interface/basic_modal_env.jl
index 2c77c775..966b8015 100644
--- a/examples/low_level_interface/basic_modal_env.jl
+++ b/examples/low_level_interface/basic_modal_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -40,7 +40,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output)
diff --git a/examples/low_level_interface/basic_modeAvg.jl b/examples/low_level_interface/basic_modeAvg.jl
index c5689fdd..e0ddb35a 100644
--- a/examples/low_level_interface/basic_modeAvg.jl
+++ b/examples/low_level_interface/basic_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -41,7 +41,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
diff --git a/examples/low_level_interface/basic_modeAvg_env.jl b/examples/low_level_interface/basic_modeAvg_env.jl
index 405b6e30..3717ef93 100644
--- a/examples/low_level_interface/basic_modeAvg_env.jl
+++ b/examples/low_level_interface/basic_modeAvg_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -33,7 +33,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
diff --git a/examples/low_level_interface/basic_modeAvg_env_THG.jl b/examples/low_level_interface/basic_modeAvg_env_THG.jl
index 72858183..a01ec518 100644
--- a/examples/low_level_interface/basic_modeAvg_env_THG.jl
+++ b/examples/low_level_interface/basic_modeAvg_env_THG.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -41,7 +41,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 import FFTW
-import PyPlot: pygui, plt
+import PythonPlot: pygui, pyplot
 
 ω = grid.ω
 t = grid.t
@@ -65,22 +65,22 @@ for ii = 1:size(Etout, 2)
 end
 
 pygui(true)
-plt.figure()
-plt.pcolormesh(FFTW.fftshift(ω, 1)./2π.*1e-15, zout, transpose(FFTW.fftshift(Ilog, 1)))
-plt.clim(-15, 0)
-plt.xlim(0.19, 1.9)
-plt.colorbar()
-
-plt.figure()
-plt.pcolormesh(t*1e15, zout, transpose(It))
-plt.colorbar()
-plt.xlim(-30, 30)
-
-plt.figure()
-plt.plot(zout.*1e2, energy.*1e6)
-plt.xlabel("Distance [cm]")
-plt.ylabel("Energy [μJ]")
-
-plt.figure()
-plt.plot(t*1e15, abs2.(Et[:, 121]))
-plt.xlim(-20, 20)
+pyplot.figure()
+pyplot.pcolormesh(FFTW.fftshift(ω, 1)./2π.*1e-15, zout, transpose(FFTW.fftshift(Ilog, 1)))
+pyplot.clim(-15, 0)
+pyplot.xlim(0.19, 1.9)
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.pcolormesh(t*1e15, zout, transpose(It))
+pyplot.colorbar()
+pyplot.xlim(-30, 30)
+
+pyplot.figure()
+pyplot.plot(zout.*1e2, energy.*1e6)
+pyplot.xlabel("Distance [cm]")
+pyplot.ylabel("Energy [μJ]")
+
+pyplot.figure()
+pyplot.plot(t*1e15, abs2.(Et[:, 121]))
+pyplot.xlim(-20, 20)
diff --git a/examples/low_level_interface/basic_modeAvg_field_noTHG.jl b/examples/low_level_interface/basic_modeAvg_field_noTHG.jl
index 33015e24..db744e42 100644
--- a/examples/low_level_interface/basic_modeAvg_field_noTHG.jl
+++ b/examples/low_level_interface/basic_modeAvg_field_noTHG.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -39,7 +39,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output)
diff --git a/examples/low_level_interface/freespace/full3D.jl b/examples/low_level_interface/freespace/full3D.jl
index 6bb58bc0..7b147d92 100644
--- a/examples/low_level_interface/freespace/full3D.jl
+++ b/examples/low_level_interface/freespace/full3D.jl
@@ -74,45 +74,45 @@ E0ωyx = FFTW.ifft(Eω[ω0idx, :, :], (1, 2));
 Iωyx = abs2.(Eωyx)
 Iωyxlog = log10.(Maths.normbymax(Iωyx));
 
-import PyPlot:pygui, plt
+import PythonPlot: pygui, pyplot
 pygui(true)
-plt.figure()
-plt.pcolormesh(x, y, (abs2.(E0ωyx)))
-plt.colorbar()
-plt.xlabel("X")
-plt.ylabel("Y")
-plt.title("I(ω=ω0, x, y, z=0)")
-
-plt.figure()
-plt.pcolormesh(x, y, (abs2.(Eωyx[ω0idx, :, :, end])))
-plt.colorbar()
-plt.xlabel("X")
-plt.ylabel("Y")
-plt.title("I(ω=ω0, x, y, z=L)")
-
-plt.figure()
-plt.pcolormesh(zout, ω.*1e-15/2π, Iωyxlog[:, N÷2+1, N÷2+1, :])
-plt.xlabel("Z (m)")
-plt.ylabel("f (PHz)")
-plt.title("I(ω, x=0, y=0, z)")
-plt.clim(-6, 0)
-plt.colorbar()
-
-plt.figure()
-plt.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, 1])
-plt.xlabel("X (mm)")
-plt.ylabel("Y (mm)")
-plt.title("\$\\int I(\\omega, x, y, z=0) d\\omega\$")
-plt.colorbar()
-
-plt.figure()
-plt.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, end])
-plt.xlabel("X (mm)")
-plt.ylabel("Y (mm)")
-plt.title("\$\\int I(\\omega, x, y, z=L) d\\omega\$")
-plt.colorbar()
-
-plt.figure()
-plt.plot(zout.*1e2, energy.*1e6)
-plt.xlabel("Distance [cm]")
-plt.ylabel("Energy [μJ]")
+pyplot.figure()
+pyplot.pcolormesh(x, y, (abs2.(E0ωyx)))
+pyplot.colorbar()
+pyplot.xlabel("X")
+pyplot.ylabel("Y")
+pyplot.title("I(ω=ω0, x, y, z=0)")
+
+pyplot.figure()
+pyplot.pcolormesh(x, y, (abs2.(Eωyx[ω0idx, :, :, end])))
+pyplot.colorbar()
+pyplot.xlabel("X")
+pyplot.ylabel("Y")
+pyplot.title("I(ω=ω0, x, y, z=L)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout, ω.*1e-15/2π, Iωyxlog[:, N÷2+1, N÷2+1, :])
+pyplot.xlabel("Z (m)")
+pyplot.ylabel("f (PHz)")
+pyplot.title("I(ω, x=0, y=0, z)")
+pyplot.clim(-6, 0)
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, 1])
+pyplot.xlabel("X (mm)")
+pyplot.ylabel("Y (mm)")
+pyplot.title("\$\\int I(\\omega, x, y, z=0) d\\omega\$")
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, end])
+pyplot.xlabel("X (mm)")
+pyplot.ylabel("Y (mm)")
+pyplot.title("\$\\int I(\\omega, x, y, z=L) d\\omega\$")
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.plot(zout.*1e2, energy.*1e6)
+pyplot.xlabel("Distance [cm]")
+pyplot.ylabel("Energy [μJ]")
diff --git a/examples/low_level_interface/freespace/full3D_env.jl b/examples/low_level_interface/freespace/full3D_env.jl
index acffcc0c..3cf30c86 100644
--- a/examples/low_level_interface/freespace/full3D_env.jl
+++ b/examples/low_level_interface/freespace/full3D_env.jl
@@ -72,45 +72,45 @@ E0ωyx = FFTW.ifft(Eω[ω0idx, :, :], (1, 2));
 Iωyx = abs2.(Eωyx);
 Iωyxlog = log10.(Maths.normbymax(Iωyx));
 
-import PyPlot:pygui, plt
+import PythonPlot: pygui, pyplot
 pygui(true)
-plt.figure()
-plt.pcolormesh(x, y, (abs2.(E0ωyx)))
-plt.colorbar()
-plt.xlabel("X")
-plt.ylabel("Y")
-plt.title("I(ω=ω0, x, y, z=0)")
-
-plt.figure()
-plt.pcolormesh(x, y, (abs2.(Eωyx[ω0idx, :, :, end])))
-plt.colorbar()
-plt.xlabel("X")
-plt.ylabel("Y")
-plt.title("I(ω=ω0, x, y, z=L)")
-
-plt.figure()
-plt.pcolormesh(zout, ω.*1e-15/2π, Iωyxlog[:, N÷2+1, N÷2+1, :])
-plt.xlabel("Z (m)")
-plt.ylabel("f (PHz)")
-plt.title("I(ω, x=0, y=0, z)")
-plt.clim(-6, 0)
-plt.colorbar()
-
-plt.figure()
-plt.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, 1])
-plt.xlabel("X (mm)")
-plt.ylabel("Y (mm)")
-plt.title("\$\\int I(\\omega, x, y, z=0) d\\omega\$")
-plt.colorbar()
-
-plt.figure()
-plt.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, end])
-plt.xlabel("X (mm)")
-plt.ylabel("Y (mm)")
-plt.title("\$\\int I(\\omega, x, y, z=L) d\\omega\$")
-plt.colorbar()
-
-plt.figure()
-plt.plot(zout.*1e2, energy.*1e6)
-plt.xlabel("Distance [cm]")
-plt.ylabel("Energy [μJ]")
+pyplot.figure()
+pyplot.pcolormesh(x, y, (abs2.(E0ωyx)))
+pyplot.colorbar()
+pyplot.xlabel("X")
+pyplot.ylabel("Y")
+pyplot.title("I(ω=ω0, x, y, z=0)")
+
+pyplot.figure()
+pyplot.pcolormesh(x, y, (abs2.(Eωyx[ω0idx, :, :, end])))
+pyplot.colorbar()
+pyplot.xlabel("X")
+pyplot.ylabel("Y")
+pyplot.title("I(ω=ω0, x, y, z=L)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout, ω.*1e-15/2π, Iωyxlog[:, N÷2+1, N÷2+1, :])
+pyplot.xlabel("Z (m)")
+pyplot.ylabel("f (PHz)")
+pyplot.title("I(ω, x=0, y=0, z)")
+pyplot.clim(-6, 0)
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, 1])
+pyplot.xlabel("X (mm)")
+pyplot.ylabel("Y (mm)")
+pyplot.title("\$\\int I(\\omega, x, y, z=0) d\\omega\$")
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, end])
+pyplot.xlabel("X (mm)")
+pyplot.ylabel("Y (mm)")
+pyplot.title("\$\\int I(\\omega, x, y, z=L) d\\omega\$")
+pyplot.colorbar()
+
+pyplot.figure()
+pyplot.plot(zout.*1e2, energy.*1e6)
+pyplot.xlabel("Distance [cm]")
+pyplot.ylabel("Energy [μJ]")
diff --git a/examples/low_level_interface/freespace/radial.jl b/examples/low_level_interface/freespace/radial.jl
index 003f9d59..9ba78a65 100644
--- a/examples/low_level_interface/freespace/radial.jl
+++ b/examples/low_level_interface/freespace/radial.jl
@@ -81,62 +81,62 @@ idcs = [argmin(abs.(zout .- point)) for point in points]
 Epoints = [Hankel.symmetric(Et[:, :, idxi], q) for idxi in idcs]
 rsym = Hankel.Rsymmetric(q);
 
-import PyPlot:pygui, plt
+import PythonPlot: pygui, pyplot
 pygui(true)
 Iλ0 = Iωr[ω0idx, :, :]
 λ0 = 2π*PhysData.c/grid.ω[ω0idx]
 w1 = w0*sqrt(1+(L/2*λ0/(π*w0^2))^2)
 # w1 = w0
 Iλ0_analytic = Maths.gauss.(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx)
-plt.figure()
-plt.plot(q.r*1e3, Maths.normbymax(Iλ0[:, end]))
-plt.plot(q.r*1e3, Maths.normbymax(Iλ0_analytic), "--")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
-plt.colorbar()
-plt.ylim(0, 4)
-plt.xlabel("z (m)")
-plt.ylabel("r (m)")
-plt.title("I(r, ω=ω0)")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
-plt.colorbar()
-plt.ylim(0, 4)
-plt.xlabel("z (m)")
-plt.ylabel("r (m)")
-plt.title("I(r, t=0)")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, Ir)
-# plt.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
-plt.colorbar()
-plt.ylim(0, R*1e3)
-# plt.clim(-6, 0)
-plt.xlabel("z (m)")
-plt.ylabel("r (m)")
-plt.title("\$\\int I(r, \\omega) d\\omega\$")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, grid.ω*1e-15/2π, Iω0log)
-plt.colorbar()
-plt.clim(0, -6)
-plt.xlabel("z (m)")
-plt.ylabel("f (PHz)")
-plt.title("I(r=0, ω)")
-
-plt.figure()
-plt.plot(zout.*1e2, energy.*1e6)
-plt.xlabel("Distance [cm]")
-plt.ylabel("Energy [μJ]")
+pyplot.figure()
+pyplot.plot(q.r*1e3, Maths.normbymax(Iλ0[:, end]))
+pyplot.plot(q.r*1e3, Maths.normbymax(Iλ0_analytic), "--")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
+pyplot.colorbar()
+pyplot.ylim(0, 4)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("r (m)")
+pyplot.title("I(r, ω=ω0)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
+pyplot.colorbar()
+pyplot.ylim(0, 4)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("r (m)")
+pyplot.title("I(r, t=0)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, Ir)
+# pyplot.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
+pyplot.colorbar()
+pyplot.ylim(0, R*1e3)
+# pyplot.clim(-6, 0)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("r (m)")
+pyplot.title("\$\\int I(r, \\omega) d\\omega\$")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, grid.ω*1e-15/2π, Iω0log)
+pyplot.colorbar()
+pyplot.clim(0, -6)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("f (PHz)")
+pyplot.title("I(r=0, ω)")
+
+pyplot.figure()
+pyplot.plot(zout.*1e2, energy.*1e6)
+pyplot.xlabel("Distance [cm]")
+pyplot.ylabel("Energy [μJ]")
 
 jw = Plotting.cmap_white("jet"; n=10)
-fig = plt.figure()
+fig = pyplot.figure()
 fig.set_size_inches(12, 4)
 for ii in 1:3
-    plt.subplot(1, 3, ii)
-    plt.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
-    plt.xlim(-42, 42)
-    plt.ylim(-1.8, 1.8)
+    pyplot.subplot(1, 3, ii)
+    pyplot.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
+    pyplot.xlim(-42, 42)
+    pyplot.ylim(-1.8, 1.8)
 end
diff --git a/examples/low_level_interface/freespace/radial_env.jl b/examples/low_level_interface/freespace/radial_env.jl
index 5cbe74e8..f7fc56e7 100644
--- a/examples/low_level_interface/freespace/radial_env.jl
+++ b/examples/low_level_interface/freespace/radial_env.jl
@@ -76,61 +76,61 @@ idcs = [argmin(abs.(zout .- point)) for point in points]
 Epoints = [Hankel.symmetric(Etout[:, :, idxi], q) for idxi in idcs]
 rsym = Hankel.Rsymmetric(q);
 
-import PyPlot:pygui, plt
+import PythonPlot: pygui, pyplot
 pygui(true)
 Iλ0 = Iωr[ω0idx, :, :]
 w1 = w0*sqrt(1+(L/2*λ0/(π*w0^2))^2)
 # w1 = w0
 Iλ0_analytic = Maths.gauss.(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx)
-plt.figure()
-plt.plot(q.r*1e3, Maths.normbymax(Iλ0[:, end]))
-plt.plot(q.r*1e3, Maths.normbymax(Iλ0_analytic), "--")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
-plt.colorbar()
-plt.ylim(0, 4)
-plt.xlabel("z (m)")
-plt.ylabel("r (m)")
-plt.title("I(r, ω=ω0)")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
-plt.colorbar()
-plt.ylim(0, 4)
-plt.xlabel("z (m)")
-plt.ylabel("r (m)")
-plt.title("I(r, t=0)")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, Ir)
-# plt.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
-plt.colorbar()
-plt.ylim(0, R*1e3)
-# plt.clim(-6, 0)
-plt.xlabel("z (m)")
-plt.ylabel("r (m)")
-plt.title("\$\\int I(r, \\omega) d\\omega\$")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, ω*1e-15/2π, Iω0log)
-plt.colorbar()
-plt.clim(0, -6)
-plt.xlabel("z (m)")
-plt.ylabel("f (PHz)")
-plt.title("I(r=0, ω)")
-
-plt.figure()
-plt.plot(zout.*1e2, energy.*1e6)
-plt.xlabel("Distance [cm]")
-plt.ylabel("Energy [μJ]")
+pyplot.figure()
+pyplot.plot(q.r*1e3, Maths.normbymax(Iλ0[:, end]))
+pyplot.plot(q.r*1e3, Maths.normbymax(Iλ0_analytic), "--")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
+pyplot.colorbar()
+pyplot.ylim(0, 4)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("r (m)")
+pyplot.title("I(r, ω=ω0)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
+pyplot.colorbar()
+pyplot.ylim(0, 4)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("r (m)")
+pyplot.title("I(r, t=0)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, Ir)
+# pyplot.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
+pyplot.colorbar()
+pyplot.ylim(0, R*1e3)
+# pyplot.clim(-6, 0)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("r (m)")
+pyplot.title("\$\\int I(r, \\omega) d\\omega\$")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, ω*1e-15/2π, Iω0log)
+pyplot.colorbar()
+pyplot.clim(0, -6)
+pyplot.xlabel("z (m)")
+pyplot.ylabel("f (PHz)")
+pyplot.title("I(r=0, ω)")
+
+pyplot.figure()
+pyplot.plot(zout.*1e2, energy.*1e6)
+pyplot.xlabel("Distance [cm]")
+pyplot.ylabel("Energy [μJ]")
 
 jw = Plotting.cmap_white("jet"; N=512, n=10)
-fig = plt.figure()
+fig = pyplot.figure()
 fig.set_size_inches(12, 4)
 for ii in 1:3
-    plt.subplot(1, 3, ii)
-    plt.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
-    plt.xlim(-42, 42)
-    plt.ylim(-1.8, 1.8)
+    pyplot.subplot(1, 3, ii)
+    pyplot.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
+    pyplot.xlim(-42, 42)
+    pyplot.ylim(-1.8, 1.8)
 end
diff --git a/examples/low_level_interface/freespace/radial_hcf.jl b/examples/low_level_interface/freespace/radial_hcf.jl
index c2a0f893..d899a88d 100644
--- a/examples/low_level_interface/freespace/radial_hcf.jl
+++ b/examples/low_level_interface/freespace/radial_hcf.jl
@@ -83,54 +83,54 @@ idcs = [argmin(abs.(zout .- point)) for point in points]
 Epoints = [Hankel.symmetric(Et[:, :, idxi], q) for idxi in idcs]
 rsym = Hankel.Rsymmetric(q);
 
-import PyPlot:pygui, plt
+import PythonPlot: pygui, pyplot
 pygui(true)
 
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
-plt.colorbar()
-plt.ylim(0, 4)
-plt.xlabel("z (cm)")
-plt.ylabel("r (mm)")
-plt.title("I(r, ω=ω0)")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
-plt.colorbar()
-plt.ylim(0, 4)
-plt.xlabel("z (cm)")
-plt.ylabel("r (mm)")
-plt.title("I(r, t=0)")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, q.r*1e3, Ir)
-# plt.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
-plt.colorbar()
-plt.ylim(0, R*1e3)
-# plt.clim(-6, 0)
-plt.xlabel("z (cm)")
-plt.ylabel("r (mm)")
-plt.title("\$\\int I(r, \\omega) d\\omega\$")
-
-plt.figure()
-plt.pcolormesh(zout*1e2, grid.ω*1e-15/2π, Iω0log)
-plt.colorbar()
-plt.clim(0, -6)
-plt.xlabel("z (cm)")
-plt.ylabel("f (PHz)")
-plt.title("I(r=0, ω)")
-
-plt.figure()
-plt.plot(zout.*1e2, energy.*1e6)
-plt.xlabel("Distance [cm]")
-plt.ylabel("Energy [μJ]")
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :])
+pyplot.colorbar()
+pyplot.ylim(0, 4)
+pyplot.xlabel("z (cm)")
+pyplot.ylabel("r (mm)")
+pyplot.title("I(r, ω=ω0)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :])
+pyplot.colorbar()
+pyplot.ylim(0, 4)
+pyplot.xlabel("z (cm)")
+pyplot.ylabel("r (mm)")
+pyplot.title("I(r, t=0)")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, q.r*1e3, Ir)
+# pyplot.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir)))
+pyplot.colorbar()
+pyplot.ylim(0, R*1e3)
+# pyplot.clim(-6, 0)
+pyplot.xlabel("z (cm)")
+pyplot.ylabel("r (mm)")
+pyplot.title("\$\\int I(r, \\omega) d\\omega\$")
+
+pyplot.figure()
+pyplot.pcolormesh(zout*1e2, grid.ω*1e-15/2π, Iω0log)
+pyplot.colorbar()
+pyplot.clim(0, -6)
+pyplot.xlabel("z (cm)")
+pyplot.ylabel("f (PHz)")
+pyplot.title("I(r=0, ω)")
+
+pyplot.figure()
+pyplot.plot(zout.*1e2, energy.*1e6)
+pyplot.xlabel("Distance [cm]")
+pyplot.ylabel("Energy [μJ]")
 
 jw = Plotting.cmap_white("jet", 512, 10)
-fig = plt.figure()
+fig = pyplot.figure()
 fig.set_size_inches(12, 4)
 for ii in 1:3
-    plt.subplot(1, 3, ii)
-    plt.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
-    plt.xlim(-42, 42)
-    plt.ylim(-1.8, 1.8)
+    pyplot.subplot(1, 3, ii)
+    pyplot.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw)
+    pyplot.xlim(-42, 42)
+    pyplot.ylim(-1.8, 1.8)
 end
diff --git a/examples/low_level_interface/full_modal/basic_modal_full.jl b/examples/low_level_interface/full_modal/basic_modal_full.jl
index fb3a137c..13336f3b 100644
--- a/examples/low_level_interface/full_modal/basic_modal_full.jl
+++ b/examples/low_level_interface/full_modal/basic_modal_full.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -41,6 +41,6 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
\ No newline at end of file
diff --git a/examples/low_level_interface/full_modal/basic_modal_full_bothpolarisations.jl b/examples/low_level_interface/full_modal/basic_modal_full_bothpolarisations.jl
index a4bde7ca..7d6ee23d 100644
--- a/examples/low_level_interface/full_modal/basic_modal_full_bothpolarisations.jl
+++ b/examples/low_level_interface/full_modal/basic_modal_full_bothpolarisations.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -41,6 +41,6 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
\ No newline at end of file
diff --git a/examples/low_level_interface/gnlse/simplescg_modeAvg.jl b/examples/low_level_interface/gnlse/simplescg_modeAvg.jl
index ce960fde..c640d743 100644
--- a/examples/low_level_interface/gnlse/simplescg_modeAvg.jl
+++ b/examples/low_level_interface/gnlse/simplescg_modeAvg.jl
@@ -1,7 +1,7 @@
 # supercontinuum from simple GNLSE parameters
 # Fig.3 of Dudley et. al, RMP 78 1135 (2006)
 
-using Luna
+using Luna, PythonPlot
 
 βs =  [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
 γ = 0.11
@@ -33,7 +33,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 #Plotting.stats(output)
 #Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
 #Plotting.time_1D(output, range(0.0, 1.0, length=5).*flength, trange=(-1e-12, 5e-12))
diff --git a/examples/low_level_interface/gnlse/simplescg_modeAvg_env.jl b/examples/low_level_interface/gnlse/simplescg_modeAvg_env.jl
index 6e3c7af6..aada5771 100644
--- a/examples/low_level_interface/gnlse/simplescg_modeAvg_env.jl
+++ b/examples/low_level_interface/gnlse/simplescg_modeAvg_env.jl
@@ -1,7 +1,7 @@
 # supercontinuum from simple GNLSE parameters
 # Fig.3 of Dudley et. al, RMP 78 1135 (2006)
 
-using Luna
+using Luna, PythonPlot
 
 βs =  [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
 γ = 0.11
@@ -33,7 +33,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 #Plotting.stats(output)
 Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
 #Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
diff --git a/examples/low_level_interface/gradients/gradient_modal.jl b/examples/low_level_interface/gradients/gradient_modal.jl
index c06d4d55..adcd0617 100644
--- a/examples/low_level_interface/gradients/gradient_modal.jl
+++ b/examples/low_level_interface/gradients/gradient_modal.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -37,7 +37,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
diff --git a/examples/low_level_interface/gradients/gradient_modal_env.jl b/examples/low_level_interface/gradients/gradient_modal_env.jl
index f763b0bf..fdd80031 100644
--- a/examples/low_level_interface/gradients/gradient_modal_env.jl
+++ b/examples/low_level_interface/gradients/gradient_modal_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -35,7 +35,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
diff --git a/examples/low_level_interface/gradients/gradient_modeAvg.jl b/examples/low_level_interface/gradients/gradient_modeAvg.jl
index d9b4ce16..4a566fa1 100644
--- a/examples/low_level_interface/gradients/gradient_modeAvg.jl
+++ b/examples/low_level_interface/gradients/gradient_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -38,6 +38,6 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop!, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
diff --git a/examples/low_level_interface/gradients/gradient_modeAvg_env.jl b/examples/low_level_interface/gradients/gradient_modeAvg_env.jl
index 7ec3f7ff..1286f6cc 100644
--- a/examples/low_level_interface/gradients/gradient_modeAvg_env.jl
+++ b/examples/low_level_interface/gradients/gradient_modeAvg_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -36,6 +36,6 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop!, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
\ No newline at end of file
diff --git a/examples/low_level_interface/mixtures/mixture_modeAvg.jl b/examples/low_level_interface/mixtures/mixture_modeAvg.jl
index 38a50857..1f65e915 100644
--- a/examples/low_level_interface/mixtures/mixture_modeAvg.jl
+++ b/examples/low_level_interface/mixtures/mixture_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -43,7 +43,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) # statsfun
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
diff --git a/examples/low_level_interface/plasma_ssfbs_modeAvg.jl b/examples/low_level_interface/plasma_ssfbs_modeAvg.jl
index d84421f2..7c5d7a06 100644
--- a/examples/low_level_interface/plasma_ssfbs_modeAvg.jl
+++ b/examples/low_level_interface/plasma_ssfbs_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 9e-6
 gas = :Ar
@@ -39,7 +39,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output, :λ, λrange=(500e-9,1800e-9), trange=(-500e-15,100e-15), dBmin=-30.0)
 Plotting.time_1D(output, [5e-2, 15e-2, 25e-2])
diff --git a/examples/low_level_interface/polarisation/elliptical.jl b/examples/low_level_interface/polarisation/elliptical.jl
index 47a45156..4d34731b 100644
--- a/examples/low_level_interface/polarisation/elliptical.jl
+++ b/examples/low_level_interface/polarisation/elliptical.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 import Luna.PhysData: wlfreq
 import Logging
 import FFTW
@@ -6,8 +6,6 @@ import NumericalIntegration: integrate, SimpsonEven
 import LinearAlgebra: mul!, ldiv!
 Logging.disable_logging(Logging.BelowMinLevel)
 
-import PyPlot:pygui, plt
-
 a = 225e-6
 gas = :Ar
 pres = 0.4
@@ -86,21 +84,21 @@ Eoutp = Array{ComplexF64,2}(undef, Eoutd[1], Eoutd[2])
 permutedims!(Eoutp, Ewp, [2, 1])
 
 pygui(true)
-plt.figure()
-plt.plot(ω./2π.*1e-15, log10.(abs2.(Eout[:,:,end])))
-plt.plot(ω./2π.*1e-15, log10.(abs2.(Eoutp)))
-#plt.ylim(-6, 0)
-plt.xlim(0,2.0)
+pyplot.figure()
+pyplot.plot(ω./2π.*1e-15, log10.(abs2.(Eout[:,:,end])))
+pyplot.plot(ω./2π.*1e-15, log10.(abs2.(Eoutp)))
+#pyplot.ylim(-6, 0)
+pyplot.xlim(0,2.0)
 
 ##
-plt.figure()
+pyplot.figure()
 It = abs2.(Maths.hilbert(FFTW.irfft(Eout[:,:,end], length(grid.t), 1)))
-plt.plot(t/1e-15, It[:,1]./maximum(It[:,1]))
-plt.plot(t/1e-15, It[:,2]./maximum(It[:,2]))
+pyplot.plot(t/1e-15, It[:,1]./maximum(It[:,1]))
+pyplot.plot(t/1e-15, It[:,2]./maximum(It[:,2]))
 Itp = abs2.(Maths.hilbert(FFTW.irfft(Eoutp, length(grid.t), 1)))
-plt.plot(t/1e-15, Itp[:,1]./maximum(Itp[:,1]))
-plt.plot(t/1e-15, Itp[:,2]./maximum(Itp[:,2]))
-plt.xlim(-50.0,50.0)
+pyplot.plot(t/1e-15, Itp[:,1]./maximum(Itp[:,1]))
+pyplot.plot(t/1e-15, Itp[:,2]./maximum(Itp[:,2]))
+pyplot.xlim(-50.0,50.0)
 
 println("$(sum(Itp)/sum(It))")
 ##
@@ -108,16 +106,16 @@ Etout = FFTW.irfft(Eout, length(grid.t), 1)
 It = Maths.normbymax(abs2.(Maths.hilbert(Etout)))
 Ilog = log10.(Maths.normbymax(abs2.(Eout)))
 for i = 1:nmodes
-    plt.figure()
-    plt.subplot(121)
-    plt.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog[:,i,:]))
-    plt.clim(-6, 0)
-    plt.xlim(0,2.0)
-    plt.colorbar()
-    plt.subplot(122)
-    plt.pcolormesh(t.*1e15, zout, transpose(It[:,i,:]))
-    plt.xlim(-30.0,100.0)
-    plt.colorbar()
+    pyplot.figure()
+    pyplot.subplot(121)
+    pyplot.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog[:,i,:]))
+    pyplot.clim(-6, 0)
+    pyplot.xlim(0,2.0)
+    pyplot.colorbar()
+    pyplot.subplot(122)
+    pyplot.pcolormesh(t.*1e15, zout, transpose(It[:,i,:]))
+    pyplot.xlim(-30.0,100.0)
+    pyplot.colorbar()
 end
 
 S = Array{Float64,3}(undef, 4, Eoutd[1], Eoutd[3])
@@ -130,17 +128,17 @@ for i = 1:Eoutd[3]
     end
 end
 
-plt.figure()
-plt.pcolormesh(zout, ω./2π.*1e-15, S[3,:,:]./S[1,:,:])
-plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(zout, ω./2π.*1e-15, S[3,:,:]./S[1,:,:])
+pyplot.colorbar()
 
-plt.figure()
-plt.pcolormesh(zout, ω./2π.*1e-15, S[2,:,:]./S[1,:,:])
-plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(zout, ω./2π.*1e-15, S[2,:,:]./S[1,:,:])
+pyplot.colorbar()
 
-plt.figure()
-plt.pcolormesh(zout, ω./2π.*1e-15, elip)
-plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(zout, ω./2π.*1e-15, elip)
+pyplot.colorbar()
 
 Eenvo = Maths.hilbert(Etout)
 Etoutd = size(Eenvo)
@@ -162,12 +160,12 @@ for i = 1:Etoutd[3]
     Elipt[i] = elipt[k,i]
 end
 
-plt.figure()
-plt.plot(zout, Ss[1,:], label="S1")
-plt.plot(zout, Ss[2,:], label="S2")
-plt.plot(zout, Ss[3,:], label="S3")
-plt.plot(zout, Ss[4,:], label="S4")
-plt.plot(zout, Elipt, label="ε")
-plt.legend()
+pyplot.figure()
+pyplot.plot(zout, Ss[1,:], label="S1")
+pyplot.plot(zout, Ss[2,:], label="S2")
+pyplot.plot(zout, Ss[3,:], label="S3")
+pyplot.plot(zout, Ss[4,:], label="S4")
+pyplot.plot(zout, Elipt, label="ε")
+pyplot.legend()
 
 
diff --git a/examples/low_level_interface/polarisation/elliptical_env.jl b/examples/low_level_interface/polarisation/elliptical_env.jl
index cdd52c1e..342971fc 100644
--- a/examples/low_level_interface/polarisation/elliptical_env.jl
+++ b/examples/low_level_interface/polarisation/elliptical_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 import StaticArrays: SMatrix, SVector
 import LinearAlgebra: mul!, ldiv!
 
@@ -61,7 +61,7 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 import FFTW
-import PyPlot:pygui, plt
+import PythonPlot: pygui, pyplot
 
 ω = FFTW.fftshift(grid.ω)
 t = grid.t
@@ -82,17 +82,17 @@ for i = 1:Eoutd[3]
 end
 
 pygui(true)
-plt.figure()
-plt.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(S[3,:,:]./S[1,:,:],1))
-plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(S[3,:,:]./S[1,:,:],1))
+pyplot.colorbar()
 
-plt.figure()
-plt.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(S[2,:,:]./S[1,:,:],1))
-plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(S[2,:,:]./S[1,:,:],1))
+pyplot.colorbar()
 
-plt.figure()
-plt.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(elip,1))
-plt.colorbar()
+pyplot.figure()
+pyplot.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(elip,1))
+pyplot.colorbar()
 
 Etout = FFTW.ifft(Eout, 1)
 It = abs2.(Etout)
@@ -102,14 +102,14 @@ Ilog = FFTW.fftshift(log10.(Maths.normbymax(abs2.(Eout))), 1)
 
 
 for i = 1:nmodes
-    plt.figure()
-    plt.subplot(121)
-    plt.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog[:,i,:]))
-    plt.clim(-6, 0)
-    plt.xlim(0,2.0)
-    plt.colorbar()
-    plt.subplot(122)
-    plt.pcolormesh(t.*1e15, zout, transpose(It[:,i,:]))
-    plt.xlim(-30.0,100.0)
-    plt.colorbar()
+    pyplot.figure()
+    pyplot.subplot(121)
+    pyplot.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog[:,i,:]))
+    pyplot.clim(-6, 0)
+    pyplot.xlim(0,2.0)
+    pyplot.colorbar()
+    pyplot.subplot(122)
+    pyplot.pcolormesh(t.*1e15, zout, transpose(It[:,i,:]))
+    pyplot.xlim(-30.0,100.0)
+    pyplot.colorbar()
 end
diff --git a/examples/low_level_interface/polarisation/modal_nonvector_plasma.jl b/examples/low_level_interface/polarisation/modal_nonvector_plasma.jl
index 8b9b5e7d..798e9864 100644
--- a/examples/low_level_interface/polarisation/modal_nonvector_plasma.jl
+++ b/examples/low_level_interface/polarisation/modal_nonvector_plasma.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 9e-6
 gas = :Ar
@@ -40,7 +40,7 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output, :λ, λrange=(500e-9,1800e-9), trange=(-500e-15,100e-15), dBmin=-30.0)
 
diff --git a/examples/low_level_interface/polarisation/modal_vector_plasma.jl b/examples/low_level_interface/polarisation/modal_vector_plasma.jl
index 396fcb23..cbf4e16b 100644
--- a/examples/low_level_interface/polarisation/modal_vector_plasma.jl
+++ b/examples/low_level_interface/polarisation/modal_vector_plasma.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 9e-6
 gas = :Ar
@@ -40,7 +40,7 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output, :λ, λrange=(500e-9,1800e-9), trange=(-500e-15,100e-15), dBmin=-30.0)
 
diff --git a/examples/low_level_interface/polarisation/modal_vector_plasma_45deg.jl b/examples/low_level_interface/polarisation/modal_vector_plasma_45deg.jl
index 63e9aad9..21497354 100644
--- a/examples/low_level_interface/polarisation/modal_vector_plasma_45deg.jl
+++ b/examples/low_level_interface/polarisation/modal_vector_plasma_45deg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 9e-6
 gas = :Ar
@@ -42,7 +42,7 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output, :λ, λrange=(500e-9,1800e-9), trange=(-500e-15,100e-15), dBmin=-30.0)
 Plotting.time_1D(output, [5e-2, 15e-2, 25e-2])
diff --git a/examples/low_level_interface/polarisation/modal_vector_plasma_CP.jl b/examples/low_level_interface/polarisation/modal_vector_plasma_CP.jl
index 27b2ddac..39504509 100644
--- a/examples/low_level_interface/polarisation/modal_vector_plasma_CP.jl
+++ b/examples/low_level_interface/polarisation/modal_vector_plasma_CP.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 9e-6
 gas = :Ar
@@ -42,7 +42,7 @@ linop = LinearOps.make_const_linop(grid, modes, λ0)
 
 Luna.run(Eω, grid, linop, transform, FT, output)
 
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output, :λ, λrange=(500e-9,1800e-9), trange=(-500e-15,100e-15), dBmin=-30.0)
 Plotting.time_1D(output, [5e-2, 15e-2, 25e-2])
diff --git a/examples/low_level_interface/rectangular/rectangular_modal.jl b/examples/low_level_interface/rectangular/rectangular_modal.jl
index 6dd84f85..878df67e 100644
--- a/examples/low_level_interface/rectangular/rectangular_modal.jl
+++ b/examples/low_level_interface/rectangular/rectangular_modal.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 50e-6
 b = 10e-6
@@ -36,7 +36,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
diff --git a/examples/low_level_interface/rectangular/rectangular_modeAvg.jl b/examples/low_level_interface/rectangular/rectangular_modeAvg.jl
index 8316eece..871c3fa7 100644
--- a/examples/low_level_interface/rectangular/rectangular_modeAvg.jl
+++ b/examples/low_level_interface/rectangular/rectangular_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 50e-6
 b = 10e-6
@@ -39,7 +39,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
diff --git a/examples/low_level_interface/rectangular/rectangular_modeAvg_env.jl b/examples/low_level_interface/rectangular/rectangular_modeAvg_env.jl
index 370928cd..a3672724 100644
--- a/examples/low_level_interface/rectangular/rectangular_modeAvg_env.jl
+++ b/examples/low_level_interface/rectangular/rectangular_modeAvg_env.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 50e-6
 b = 10e-6
@@ -39,7 +39,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
diff --git a/examples/low_level_interface/run_all.jl b/examples/low_level_interface/run_all.jl
index b2e43b7c..d5e9c42c 100644
--- a/examples/low_level_interface/run_all.jl
+++ b/examples/low_level_interface/run_all.jl
@@ -1,7 +1,7 @@
 using Luna
 Luna.set_fftw_mode(:estimate)
 import Logging: @warn, disable_logging, Info
-import PyPlot: plt
+import PythonPlot: pyplot
 
 disable_logging(Info) # disable progress logging to avoid filling the shell window
 
@@ -14,7 +14,7 @@ function run(item)
             bt = catch_backtrace()
             @warn sprint(showerror, e, bt)
         finally
-            plt.close("all")
+            pyplot.close("all")
         end
     elseif isdir(item)
         for iitem in readdir(item)
diff --git a/examples/low_level_interface/stepindex/step_modeAvg_env.jl b/examples/low_level_interface/stepindex/step_modeAvg_env.jl
index b40e8a50..846c0046 100644
--- a/examples/low_level_interface/stepindex/step_modeAvg_env.jl
+++ b/examples/low_level_interface/stepindex/step_modeAvg_env.jl
@@ -1,6 +1,6 @@
 # propagation in step index fibre
 
-using Luna
+using Luna, PythonPlot
 
 # single mode fibre at 1030 nm
 a = 5e-6
@@ -32,7 +32,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 #Plotting.stats(output)
 Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(980e-9, 1200e-9), trange=(-2e-12, 2e-12))
 #Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
diff --git a/examples/low_level_interface/stepindex/stepscg_modeAvg.jl b/examples/low_level_interface/stepindex/stepscg_modeAvg.jl
index d2b7301c..2c715243 100644
--- a/examples/low_level_interface/stepindex/stepscg_modeAvg.jl
+++ b/examples/low_level_interface/stepindex/stepscg_modeAvg.jl
@@ -1,6 +1,6 @@
 # supercontinuum in strand of silica in air
 
-using Luna
+using Luna, PythonPlot
 
 # single mode fibre at 1030 nm
 a = 1.25e-6
@@ -31,7 +31,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 #Plotting.stats(output)
 #Plotting.prop_2D(output)
 #Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
diff --git a/examples/low_level_interface/stepindex/stepscg_modeAvg_env.jl b/examples/low_level_interface/stepindex/stepscg_modeAvg_env.jl
index 1fafbae8..8745071d 100644
--- a/examples/low_level_interface/stepindex/stepscg_modeAvg_env.jl
+++ b/examples/low_level_interface/stepindex/stepscg_modeAvg_env.jl
@@ -1,6 +1,6 @@
 # supercontinuum in strand of silica in air
 
-using Luna
+using Luna, PythonPlot
 
 # single mode fibre at 1030 nm
 a = 1.25e-6
@@ -31,7 +31,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 #Plotting.stats(output)
 Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
 #Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
diff --git a/examples/low_level_interface/stepindex/stepsimplecmp.jl b/examples/low_level_interface/stepindex/stepsimplecmp.jl
index 12a083a0..c8d68e10 100644
--- a/examples/low_level_interface/stepindex/stepsimplecmp.jl
+++ b/examples/low_level_interface/stepindex/stepsimplecmp.jl
@@ -1,6 +1,6 @@
 using Luna
 using Luna.PhysData: Polynomials
-import PyPlot: plt
+import PythonPlot: pyplot
 
 a = 1.25e-6
 flength = 15e-2
@@ -42,12 +42,12 @@ p = Polynomials.fit(1e-15*(ωs .- ω0), βs, order)
 βcoeffs = p.coeffs .* 1e-15 .^ (collect(0:order)) .* factorial.(0:order)
 
 s = SimpleFibre.SimpleMode(ω0, βcoeffs)
-plt.figure()
-plt.plot(ωs .- ω0, Modes.β_ret.(m, ωs; λ0), label="Full")
-plt.plot(ωs .- ω0, Modes.β_ret(s, ωs; λ0), "--", label="Polynomial fit")
-plt.legend()
-plt.xlabel("ω-ω0 (rad/s)")
-plt.ylabel("β (1/m)")
+pyplot.figure()
+pyplot.plot(ωs .- ω0, Modes.β_ret.(m, ωs; λ0), label="Full")
+pyplot.plot(ωs .- ω0, Modes.β_ret(s, ωs; λ0), "--", label="Polynomial fit")
+pyplot.legend()
+pyplot.xlabel("ω-ω0 (rad/s)")
+pyplot.ylabel("β (1/m)")
 
 N0, n0, n2 = Tools.getN0n0n2(ω0, :SiO2)
 γ = Tools.getγ(ω0, m, n2)
@@ -76,11 +76,11 @@ Plotting.prop_2D(outputs)
 ##
 λ, Iλm = Processing.getIω(outputm, :λ, flength)
 _, Iλs = Processing.getIω(outputs, :λ, flength)
-plt.figure()
-plt.semilogy(1e9λ, Iλm, label="step-index mode")
-plt.semilogy(1e9λ, Iλs, "--",  label="GNLSE approximation")
-plt.legend()
-plt.xlabel("Wavelength (nm)")
-plt.ylabel("SED (a.u.)")
+pyplot.figure()
+pyplot.semilogy(1e9λ, Iλm, label="step-index mode")
+pyplot.semilogy(1e9λ, Iλs, "--",  label="GNLSE approximation")
+pyplot.legend()
+pyplot.xlabel("Wavelength (nm)")
+pyplot.ylabel("SED (a.u.)")
 
 # close but not exact. This is because we cannot fully cancel the frequency dependence of neff.
\ No newline at end of file
diff --git a/examples/low_level_interface/tapers/taper_modal.jl b/examples/low_level_interface/tapers/taper_modal.jl
index 39864bdd..ee32164e 100644
--- a/examples/low_level_interface/tapers/taper_modal.jl
+++ b/examples/low_level_interface/tapers/taper_modal.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -46,7 +46,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
diff --git a/examples/low_level_interface/tapers/taper_modeAvg.jl b/examples/low_level_interface/tapers/taper_modeAvg.jl
index 9da94c05..f7fada10 100644
--- a/examples/low_level_interface/tapers/taper_modeAvg.jl
+++ b/examples/low_level_interface/tapers/taper_modeAvg.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 a = 13e-6
 gas = :Ar
@@ -45,7 +45,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
 Luna.run(Eω, grid, linop, transform, FT, output)
 
 ##
-Plotting.pygui(true)
+
 Plotting.stats(output)
 Plotting.prop_2D(output)
 Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
diff --git a/examples/simple_interface/2stage_compressor_RDWemission.jl b/examples/simple_interface/2stage_compressor_RDWemission.jl
index 2b36abf9..eb1337ca 100644
--- a/examples/simple_interface/2stage_compressor_RDWemission.jl
+++ b/examples/simple_interface/2stage_compressor_RDWemission.jl
@@ -2,7 +2,7 @@
 We then optimise the pulse compression, plot the resulting compressed pulse,
 and pass the compressed pulse on to a second simulation for UV RDW emission.
 =#
-using Luna
+using Luna, PythonPlot
 
 λ0 = 800e-9
 τfwhm = 35e-15
diff --git a/examples/simple_interface/RDWemission_DUV.jl b/examples/simple_interface/RDWemission_DUV.jl
index 7825597e..7f5fbb28 100644
--- a/examples/simple_interface/RDWemission_DUV.jl
+++ b/examples/simple_interface/RDWemission_DUV.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 radius = 125e-6 # HCF core radius
 flength = 3 # HCF length
@@ -21,4 +21,4 @@ Plotting.time_1D(duv; modes=:sum)
 Plotting.time_1D(duv; modes=1, bandpass=(220e-9, 270e-9))
 # plot the spectrogram of the pulse at the exit with a white background
 Plotting.spectrogram(duv, flength; trange=(-20e-15, 30e-15), λrange=(150e-9, 1000e-9),
-                     N=256, fw=3e-15, cmap=Plotting.cmap_white("viridis"))
\ No newline at end of file
+                     N=256, fw=3e-15)
\ No newline at end of file
diff --git a/examples/simple_interface/RDWemission_VUV_gradient.jl b/examples/simple_interface/RDWemission_VUV_gradient.jl
index ac3a537a..2ca58260 100644
--- a/examples/simple_interface/RDWemission_VUV_gradient.jl
+++ b/examples/simple_interface/RDWemission_VUV_gradient.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 radius = 125e-6 # HCF core radius
 flength = 1.5 # HCF length
diff --git a/examples/simple_interface/RamanSSFS.jl b/examples/simple_interface/RamanSSFS.jl
index 4f0bfeac..81f3a632 100644
--- a/examples/simple_interface/RamanSSFS.jl
+++ b/examples/simple_interface/RamanSSFS.jl
@@ -1,4 +1,4 @@
-using Luna
+using Luna, PythonPlot
 
 radius = 13e-6 # HCF core radius
 flength = 2 # HCF length
diff --git a/examples/simple_interface/fourwavemixing_opposite_helicity.jl b/examples/simple_interface/fourwavemixing_opposite_helicity.jl
index b85a2a8b..d7fa71e8 100644
--- a/examples/simple_interface/fourwavemixing_opposite_helicity.jl
+++ b/examples/simple_interface/fourwavemixing_opposite_helicity.jl
@@ -1,5 +1,4 @@
-using Luna
-import PyPlot: plt
+using Luna, PythonPlot
 
 #= In this example we simulate degenerate four-wave mixing between a circularly polarised
 pump pulse at 400 nm and a seed pulse at 800 nm to generate a circularly polarised idler
@@ -48,12 +47,12 @@ Iλequal = dropdims(sum(Iλequal; dims=2); dims=(2, 3))
 Iλopposite = dropdims(sum(Iλopposite; dims=2); dims=(2, 3))
 
 # plot the result
-plt.figure()
-plt.plot(λ*1e9, Iλequal*1e-3, label="Equal helicity")
-plt.plot(λ*1e9, Iλopposite*1e-3, label="Opposite helicity")
-plt.xlim(200, 950)
-plt.xlabel("Wavelength (nm)")
-plt.ylabel("Spectral energy density (μJ/nm)")
-plt.ylim(ymin=0)
-plt.legend()
+pyplot.figure()
+pyplot.plot(λ*1e9, Iλequal*1e-3, label="Equal helicity")
+pyplot.plot(λ*1e9, Iλopposite*1e-3, label="Opposite helicity")
+pyplot.xlim(200, 950)
+pyplot.xlabel("Wavelength (nm)")
+pyplot.ylabel("Spectral energy density (μJ/nm)")
+pyplot.ylim(ymin=0)
+pyplot.legend()
 
diff --git a/examples/simple_interface/gnlse_scg.jl b/examples/simple_interface/gnlse_scg.jl
index c6749000..79032e1d 100644
--- a/examples/simple_interface/gnlse_scg.jl
+++ b/examples/simple_interface/gnlse_scg.jl
@@ -14,6 +14,6 @@ power = 10000.0
 output = prop_gnlse(γ, flength, βs; λ0, τfwhm, power, pulseshape=:sech, λlims=(400e-9, 2400e-9), trange=12.5e-12, ramanmodel=:sdo, τ1=12.2e-15, τ2=32e-15)
 
 ##
-Plotting.pygui(true)
+
 Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
 Plotting.spec_1D(output, range(0.0, 1.0, length=5).*flength, λrange=(400e-9, 1300e-9))
diff --git a/examples/simple_interface/gnlse_sol.jl b/examples/simple_interface/gnlse_sol.jl
index 79c6321f..ac9bb6d0 100644
--- a/examples/simple_interface/gnlse_sol.jl
+++ b/examples/simple_interface/gnlse_sol.jl
@@ -20,5 +20,5 @@ output = prop_gnlse(γ, flength, βs; λ0, τfwhm, power=P0, pulseshape=:sech, 
                     raman=false, shock=false, fr, shotnoise=false)
 
 ##
-Plotting.pygui(true)
+
 Plotting.prop_2D(output, :ω, dBmin=-100.0,  λrange=(720e-9,1000e-9), trange=(-300e-15, 300e-15), oversampling=1)
diff --git a/examples/simple_interface/gnlse_ssfs.jl b/examples/simple_interface/gnlse_ssfs.jl
index 46e3cc65..f118e691 100644
--- a/examples/simple_interface/gnlse_ssfs.jl
+++ b/examples/simple_interface/gnlse_ssfs.jl
@@ -18,5 +18,5 @@ output = prop_gnlse(γ, flength, βs; λ0, τfwhm=1.763*τ0, power=P0, pulseshap
                     saveN=601)
 
 ##
-Plotting.pygui(true)
+
 Plotting.prop_2D(output, :ω, dBmin=-40.0,  λrange=(700e-9,1500e-9), trange=(-50e-15, 6000e-15), oversampling=1)
diff --git a/examples/simple_interface/scan.jl b/examples/simple_interface/scan.jl
index e5a49aa6..04ceeba3 100644
--- a/examples/simple_interface/scan.jl
+++ b/examples/simple_interface/scan.jl
@@ -1,5 +1,4 @@
-using Luna
-import PyPlot: plt
+using Luna, PythonPlot
 
 # Fixed parameters:
 a = 125e-6
@@ -39,7 +38,7 @@ end
     Processing.Common(λ), Iλ[:, end], zstat, edens, max_peakpower
 end
 
-fig, axs = plt.subplots(1, length(pressures))
+fig, axs = pyplot.subplots(1, length(pressures))
 fig.set_size_inches(8, 2)
 for (pidx, pressure) in enumerate(pressures)
     ax = axs[pidx]
@@ -50,9 +49,9 @@ for (pidx, pressure) in enumerate(pressures)
     ax.set_title("Pressure: $pressure bar")
     ax.set_xlim(100, 1200)
 end
-plt.colorbar(img, ax=axs, label="Energy density (dB)")
+pyplot.colorbar(img, ax=axs, label="Energy density (dB)")
 
-fig, axs = plt.subplots(1, length(pressures))
+fig, axs = pyplot.subplots(1, length(pressures))
 fig.set_size_inches(8, 2)
 cols = Plotting.cmap_colours(length(energies))
 edmax = maximum(maximum.(edens))
@@ -69,12 +68,12 @@ for (pidx, pressure) in enumerate(pressures)
     ax.set_title("Pressure: $pressure bar")
 end
 
-fig = plt.figure()
+fig = pyplot.figure()
 for (pidx, pressure) in enumerate(pressures)
-    plt.plot(energies*1e6, max_peakpower[:, pidx], label="$pressure bar")
+    pyplot.plot(energies*1e6, max_peakpower[:, pidx], label="$pressure bar")
 end
-plt.xlim(1e6.*extrema(energies))
-plt.ylim_ymin=0
-plt.xlabel("Energy (μJ)")
-plt.ylabel("Maximum peak power (W)")
-plt.legend()
+pyplot.xlim(1e6.*extrema(energies))
+pyplot.ylim_ymin=0
+pyplot.xlabel("Energy (μJ)")
+pyplot.ylabel("Maximum peak power (W)")
+pyplot.legend()
diff --git a/examples/simple_interface/selfcompression_dispersion.jl b/examples/simple_interface/selfcompression_dispersion.jl
index a0a1f27b..835bb075 100644
--- a/examples/simple_interface/selfcompression_dispersion.jl
+++ b/examples/simple_interface/selfcompression_dispersion.jl
@@ -1,5 +1,4 @@
-using Luna
-import PyPlot: plt
+using Luna, PythonPlot
 
 #= In this example we simulate the effect of pump pulse dispersion on soliton self-compression
 and RDW emission. To do this we use the "ϕ" and "propagator" keyword arguments to
@@ -31,12 +30,12 @@ outputs = [prop_capillary(a, flength, gas, pressure; trange=1e-12, λlims=(150e-
 λ, _ = Processing.getIω(outputs[1], :λ)
 Iλs = [Processing.getIω(oi, :λ)[2] for oi in outputs]
 
-plt.figure()
+pyplot.figure()
 for (d, Iλ) in zip(thicknesses, Iλs)
-    plt.semilogy(λ*1e9, 1e-3*Iλ[:, end], label="$(d*1e3) mm")
+    pyplot.semilogy(λ*1e9, 1e-3*Iλ[:, end], label="$(d*1e3) mm")
 end
-plt.xlim(130, 1100)
-plt.ylim(1e-5, 10)
-plt.xlabel("Wavelength (nm)")
-plt.ylabel("Spectral energy density (μJ/nm)")
-plt.legend()
\ No newline at end of file
+pyplot.xlim(130, 1100)
+pyplot.ylim(1e-5, 10)
+pyplot.xlabel("Wavelength (nm)")
+pyplot.ylabel("Spectral energy density (μJ/nm)")
+pyplot.legend()
\ No newline at end of file
diff --git a/ext/GLMakieExt.jl b/ext/GLMakieExt.jl
new file mode 100644
index 00000000..a2e7ce57
--- /dev/null
+++ b/ext/GLMakieExt.jl
@@ -0,0 +1,638 @@
+module GLMakieExt
+import Luna: Grid, Maths, PhysData, Processing
+import Luna.PhysData: wlfreq, c, ε_0
+import Luna.Output: AbstractOutput
+import Luna.Processing: makegrid, getIω, getEω, getEt, nearest_z
+import FFTW
+import Printf: @sprintf
+import Base: display
+import GLMakie
+
+GLMakie.set_theme!(GLMakie.Theme(fontsize = 40))
+
+function newfig(; size=(800,600))
+    f = GLMakie.Figure(; size)
+    display(GLMakie.Screen(), f)
+    f
+end
+
+"""
+    cmap_white(cmap, N=512, n=8)
+
+Replace the lowest colour stop of `cmap` (after splitting into `n` stops) with white and
+create a new colourmap with `N` stops.
+"""
+function cmap_white(cmap; N=2^12, n=8)
+    insert!(GLMakie.to_colormap(:viridis), 1, GLMakie.Colors.RGBA(Float32(1),1,1,1))
+end
+
+"""
+    cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
+
+Make an array of `num` different colours that follow the colourmap `cmap` between the values
+`cmin` and `cmax`.
+"""
+function cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
+    cm = ColorMap(cmap)
+    n = collect(range(cmin, cmax; length=num))
+    cm.(n)
+end
+
+"""
+    subplotgrid(N, portrait=true, kwargs...)
+
+Create a figure with `N` subplots laid out in a grid that is as close to square as possible.
+If `portrait` is `true`, try to lay out the grid in portrait orientation (taller than wide),
+otherwise landscape (wider than tall).
+"""
+function subplotgrid(N, portrait=true; colw=350, rowh=250, title=nothing)
+    cols = ceil(Int, sqrt(N))
+    rows = ceil(Int, N/cols)
+    collect(Iterators.product(1:rows,1:cols)), cols*colw, rows*rowh
+end
+
+"""
+    get_modes(output)
+
+Determine whether `output` contains a multimode simulation, and if so, return the names
+of the modes.
+"""
+function get_modes(output)
+    t = output["simulation_type"]["transform"]
+    !startswith(t, "TransModal") && return false, nothing
+    lines = split(t, "\n")
+    modeline = findfirst(li -> startswith(li, "  modes:"), lines)
+    endline = findnext(li -> !startswith(li, " "^4), lines, modeline+1)
+    mlines = lines[modeline+1 : endline-1]
+    labels = [match(r"{([^,]*),", li).captures[1] for li in mlines]
+    angles = zeros(length(mlines))
+    for (ii, li) in enumerate(mlines)
+        m = match(r"ϕ=(-?[0-9]+.[0-9]+)π", li)
+        isnothing(m) && continue # no angle information in mode label)
+        angles[ii] = parse(Float64, m.captures[1])
+    end
+    if !all(angles .== 0)
+        for i in eachindex(labels)
+            if startswith(labels[i], "HE")
+                if angles[i] == 0
+                    θs = "x"
+                elseif angles[i] == 0.5
+                    θs = "y"
+                else
+                    θs = "$(angles[i])π"
+                end
+                labels[i] *= " ($θs)"
+            end
+        end
+    end
+    return true, labels
+end
+
+"""
+    stats(output; kwargs...)
+
+Plot all statistics available in `output`. Additional `kwargs` are passed onto `plt.plot()`
+"""
+function stats(output; kwargs...)
+    stats = output["stats"]
+
+    pstats = [] # pulse statistics
+    haskey(stats, "energy") && push!(pstats, (1e6*stats["energy"], "Energy (μJ)"))
+    for (k, v) in pairs(stats)
+        startswith(k, "energy_") || continue
+        str = "Energy "*replace(k[8:end], "_" => " ")*" (μJ)"
+        push!(pstats, (1e6*stats[k], str))
+    end
+    for (k, v) in pairs(stats)
+        startswith(k, "peakpower_") || continue
+        Pfac, unit = power_unit(stats[k])
+        str = "Peak power "*replace(k[11:end], "_" => " ")*" ($unit)"
+        push!(pstats, (Pfac*stats[k], str))
+    end
+    if haskey(stats, "peakpower")
+        Pfac, unit = power_unit(stats["peakpower"])
+        push!(pstats, (Pfac*stats["peakpower"], "Peak power ($unit)"))
+    end
+    haskey(stats, "peakintensity") && push!(
+        pstats, (1e-16*stats["peakintensity"], "Peak Intensity (TW/cm²)"))
+    haskey(stats, "fwhm_t_min") && push!(pstats, (1e15*stats["fwhm_t_min"], "min FWHM (fs)"))
+    haskey(stats, "fwhm_t_max") && push!(pstats, (1e15*stats["fwhm_t_max"], "max FWHM (fs)"))
+    haskey(stats, "fwhm_r") && push!(pstats, (1e6*stats["fwhm_r"], "Radial FWHM (μm)"))
+    haskey(stats, "ω0") && push!(pstats, (1e9*wlfreq.(stats["ω0"]), "Central wavelength (nm)"))
+
+    fstats = [] # fibre/waveguide/propagation statistics
+    if haskey(stats, "electrondensity")
+        push!(fstats, (1e-6*stats["electrondensity"], "Electron density (cm⁻³)"))
+        if haskey(stats, "density")
+            push!(fstats,
+                 (100*stats["electrondensity"]./stats["density"], "Ionisation fraction (%)"))
+        end
+    end
+    haskey(stats, "density") && push!(
+        fstats, (1e-6*stats["density"], "Density (cm⁻³)"))
+    haskey(stats, "pressure") && push!(
+        fstats, (stats["pressure"], "Pressure (bar)"))
+    haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)"))
+    haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)"))
+    haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)"))
+
+    z = stats["z"]*1e2
+
+    multimode, modes = get_modes(output)
+    modes = isnothing(modes) ? [""] : modes
+
+    Npl = length(pstats)
+    if Npl > 0
+        idcs, width, height = subplotgrid(Npl)
+        pfig = newfig(size=(width, height))
+        for n in 1:Npl
+            data, ylabel = pstats[n]
+            scale = (multimode ? log10 : identity)
+            data = (multimode && ndims(data) > 1) ? data : data'
+            data = multimode ? max.(data, 1e-300) : data
+            ax = GLMakie.Axis(pfig[idcs[n]...]; xlabel="Distance (cm)", ylabel, yscale=scale)
+            println("$(size(data))   $modes")
+            for i in 1:size(data,1)
+                GLMakie.lines!(z, data[i,:], label=modes[i])
+            end
+            multimode && (ndims(data) > 1) && GLMakie.axislegend(framevisible=false)
+        end
+        GLMakie.DataInspector()
+    end
+    
+    Npl = length(fstats)
+    if Npl > 0
+        idcs, width, height = subplotgrid(Npl)
+        ffig = newfig(size=(width, height))
+        for n in 1:Npl
+            data, ylabel = fstats[n]
+            scale = ((multimode && should_log10(data)) ? log10 : identity)
+            data = (multimode && ndims(data) > 1) ? data : data'
+            data = (multimode && should_log10(data)) ? max.(data, 1e-300) : data
+            ax = GLMakie.Axis(ffig[idcs[n]...]; xlabel="Distance (cm)", ylabel, yscale=scale)
+            for i in 1:size(data,1)
+                GLMakie.lines!(z, data[i,:], label=modes[i])
+            end
+            multimode && (ndims(data) > 1) && GLMakie.axislegend(framevisible=false)
+        end
+        GLMakie.DataInspector()
+    end
+    [pfig, ffig]
+end
+
+"""
+    should_log10(A, tolfac=10)
+
+For multi-line plots, determine whether data for different lines contained in `A` spans
+a sufficiently large range that a logarithmic scale should be used. By default, this is the
+case when there is any point where the lines are different by more than a factor of 10.
+"""
+function should_log10(A, tolfac=10)
+    mi = minimum(A; dims=2)
+    ma = maximum(A; dims=2)
+    any(ma./mi .> 10)
+end
+
+window_str(::Nothing) = ""
+window_str(win::NTuple{4, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win[2:3]...)
+window_str(win::NTuple{2, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win...)
+window_str(window) = "custom bandpass"
+
+"""
+    prop_2D(output, specaxis=:f)
+
+Make false-colour propagation plots for `output`, using spectral x-axis `specaxis` (see
+[`getIω`](@ref)). For multimode simulations, create one figure for each mode plus one for
+the sum of all modes.
+
+# Keyword arguments
+- `λrange::Tuple(Float64, Float64)` : x-axis limits for spectral plot (wavelength in metres)
+- `trange::Tuple(Float64, Float64)` : x-axis limits for time-domain plot (time in seconds)
+- `dBmin::Float64` : lower colour-scale limit for logarithmic spectral plot
+- `resolution::Real` smooth the spectral energy density as defined by [`getIω`](@ref).
+"""
+function prop_2D(output, specaxis=:f;
+                 trange=(-50e-15, 50e-15), bandpass=nothing,
+                 λrange=(150e-9, 2000e-9), dBmin=-60,
+                 resolution=nothing, modes=nothing, oversampling=4,
+                 kwargs...)
+    z = output["z"]*1e2
+    if specaxis == :λ
+            specx, Iω = getIω(output, specaxis, specrange=λrange, resolution=resolution)
+    else
+            specx, Iω = getIω(output, specaxis, resolution=resolution)
+    end
+
+    t, Et = getEt(output; trange, bandpass, oversampling)
+    It = abs2.(Et)
+
+    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
+    specx .*= specxfac
+
+    multimode, modelabels = get_modes(output)
+
+    if multimode
+        fig = _prop2D_mm(modelabels, modeidcs(modes, modelabels), t, z, specx, It, Iω,
+                         speclabel, speclims, trange, dBmin, window_str(bandpass);
+                         kwargs...)
+    else
+        fig = _prop2D_sm(t, z, specx, It, Iω,
+                         speclabel, speclims, trange, dBmin, window_str(bandpass);
+                         kwargs...)
+    end
+    fig
+end
+
+modeidcs(m::Int, ml) = [m]
+modeidcs(m::Symbol, ml) = (m == :sum) ? [] : error("modes must be :sum, a single integer, or iterable")
+modeidcs(m::Nothing, ml) = 1:length(ml)
+modeidcs(m, ml) = m
+
+# Helper function to convert λrange to the correct numbers depending on specaxis
+function getspeclims(λrange, specaxis)
+    if specaxis == :f
+        specxfac = 1e-15
+        speclims = (specxfac*c/maximum(λrange), specxfac*c/minimum(λrange))
+        speclabel = "Frequency (PHz)"
+    elseif specaxis == :ω
+        specxfac = 1e-15
+        speclims = (specxfac*wlfreq(maximum(λrange)), specxfac*wlfreq(minimum(λrange)))
+        speclabel = "Angular frequency (rad/fs)"
+    elseif specaxis == :λ
+        specxfac = 1e9
+        speclims = λrange .* specxfac
+        speclabel = "Wavelength (nm)"
+    else
+        error("Unknown specaxis $specaxis")
+    end
+    return speclims, speclabel, specxfac
+end
+
+# single-mode 2D propagation plots
+function _prop2D_sm(t, z, specx, It, Iω, speclabel, speclims, trange, dBmin, bpstr; kwargs...)
+    id = "($(string(hash(gensym()); base=16)[1:4])) "
+    num = id * "Propagation" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+    Iω = Maths.normbymax(Iω)
+    _prop2D_fig(num, specx, z, Iω, dBmin, speclabel, speclims, t, It, trange; kwargs...)
+end
+
+# multi-mode 2D propagation plots
+function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
+                    speclabel, speclims, trange, dBmin, bpstr;
+                    kwargs...)
+    pfigs = []
+    Iω = Maths.normbymax(Iω)
+    id = "($(string(hash(gensym()); base=16)[1:4])) "
+    for mi in modes
+        num = id * "Propagation ($(modelabels[mi]))" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+        pfig = _prop2D_fig(num, specx, z, Iω[:, mi, :], dBmin, speclabel, speclims, t, It[:, mi, :], trange; kwargs...)
+        push!(pfigs, pfig)
+    end
+
+    num = id * "Propagation (all modes)" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+    Iωall = dropdims(sum(Iω, dims=2), dims=2)
+    Itall = dropdims(sum(It, dims=2), dims=2)
+    pfig = _prop2D_fig(num, specx, z, Iωall, dBmin, speclabel, speclims, t, Itall, trange; kwargs...)
+    push!(pfigs, pfig)
+    return pfigs
+end
+
+function _prop2D_fig(name, specx, z, Iω, dBmin, speclabel, speclims, t, It, trange; kwargs...)
+    if haskey(kwargs, :cmap)
+        cmap = kwargs[:cmap]
+    else
+        cmap = :viridis
+    end
+    pfig = newfig(size=(1000,400))
+    ax, hm = GLMakie.heatmap(pfig[1,1], specx, z, 10*log10.(Iω),
+                             colorrange=(dBmin,0), interpolate=true,
+                             colormap=cmap,
+                             axis=(; xlabel=speclabel, ylabel="Distance (cm)"))
+    GLMakie.xlims!(ax, speclims)
+    cb = GLMakie.Colorbar(pfig[1, 2], hm, label="SED (dB)")
+
+    Pfac, unit = power_unit(It)
+    ax, hm = GLMakie.heatmap(pfig[1,3], t*1e15, z, Pfac .* It,
+                             interpolate=true,
+                             colormap=cmap,
+                             axis=(; xlabel="Time (fs)", ylabel="Distance (cm)"))
+    GLMakie.xlims!(ax, trange.*1e15)
+    cb = GLMakie.Colorbar(pfig[1, 4], hm, label="Power ($unit)")
+    GLMakie.DataInspector()
+    pfig
+end
+
+
+"""
+    time_1D(output, zslice, y=:Pt, kwargs...)
+
+Create lineplots of time-domain slice(s) of the propagation.
+
+The keyword argument `y` determines
+what is plotted: `:Pt` (power, default), `:Esq` (squared electric field) or `:Et` (electric field).
+
+The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
+a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
+
+The keyword argument `oversampling` determines the amount of oversampling done before plotting.
+
+Other `kwargs` are passed onto `plt.plot`.
+"""
+function time_1D(output, zslice=maximum(output["z"]);
+                y=:Pt, modes=nothing,
+                oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing,
+                FTL=false, propagate=nothing,
+                kwargs...)
+    t, Et, zactual = getEt(output, zslice,
+                           trange=trange, oversampling=oversampling, bandpass=bandpass,
+                           FTL=FTL, propagate=propagate)
+    if y == :Pt
+        yt = abs2.(Et)
+    elseif y == :Et
+        yt = real(Et)
+    elseif y == :Esq
+        yt = real(Et).^2
+    else
+        error("unknown time plot variable $y")
+    end
+    multimode, modestrs = get_modes(output)
+    if multimode
+        if modes == :sum
+            y == :Pt || error("Modal sum can only be plotted for power!")
+            yt = dropdims(sum(yt, dims=2), dims=2)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            yt = yt[:, modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    end
+
+    yfac, unit = power_unit(abs2.(Et), y)
+
+    xlabel = "Time (fs)"
+    ylabel = y == :Et ?  "Field ($unit)" : "Power ($unit)"
+    if multimode && nmodes > 1
+       sfig = _plot_slice_mm(t*1e15, yfac*yt, zactual, modestrs, xlabel, ylabel, fwlabel=true)
+    else
+        sfig = newfig()
+        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
+        label = multimode ? zs.*" ($modestrs)" : zs
+        GLMakie.Axis(sfig[1, 1]; xlabel, ylabel)
+        for iz in eachindex(zactual)
+            fw = Maths.fwhm(t*1e15, yfac*yt[:, iz])
+            label=(label[iz] * @sprintf(" [%.2f %s]", fw, "fs"))
+            GLMakie.lines!(t*1e15, yfac*yt[:, iz]; label)
+        end
+    end
+    GLMakie.axislegend(framevisible=false)
+    GLMakie.xlims!((1e15.*trange)...)
+    y == :Et || GLMakie.ylims!(low=0)
+    GLMakie.DataInspector()
+    sfig  
+end
+
+# Automatically find power unit depending on scale of electric field.
+function power_unit(Pt, y=:Pt)
+    units = ["kW", "MW", "GW", "TW", "PW"]
+    Pmax = maximum(Pt)
+    oom = clamp(floor(Int, log10(Pmax)/3), 1, 5) # maximum unit is PW
+    powerfac = 1/10^(oom*3)
+    if y == :Et
+        sqrt(powerfac), "$(units[oom])\$^{1/2}\$"
+    else
+        return powerfac, units[oom]
+    end
+end    
+
+"""
+    spec_1D(output, zslice, specaxis=:λ, log10=true, log10min=1e-6)
+
+Create lineplots of spectral-domain slices of the propagation.
+
+The x-axis is determined by `specaxis` (see [`getIω`](@ref)).
+
+If `log10` is true, plot on a logarithmic scale, with a y-axis range of `log10min`. 
+
+The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
+a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
+
+Other `kwargs` are passed onto `plt.plot`.
+"""
+function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
+                 modes=nothing, λrange=(150e-9, 1200e-9),
+                 log10=true, log10min=1e-6, resolution=nothing,
+                 kwargs...)
+    if specaxis == :λ
+        specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution)
+    else
+        specx, Iω, zactual = getIω(output, specaxis, zslice, resolution=resolution)
+    end
+    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
+    multimode, modestrs = get_modes(output)
+    if multimode
+        modes = isnothing(modes) ? (1:size(Iω, 2)) : modes
+        if modes == :sum
+            Iω = dropdims(sum(Iω, dims=2), dims=2)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            Iω = Iω[:, modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    end
+
+    specx .*= specxfac
+
+    if multimode && nmodes > 1
+        sfig = _plot_slice_mm(specx, Iω, zactual, modestrs, speclabel, "Spectral energy density", log10)
+    else
+        sfig = newfig()
+        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
+        label = multimode ? zs.*" ($modestrs)" : zs
+        scale = (log10 ? Base.log10 : :identity)
+        GLMakie.Axis(sfig[1, 1], yscale = scale, xlabel=speclabel, ylabel="Spectral energy density")
+        for iz in eachindex(zactual)
+            GLMakie.lines!(specx, Iω[:, iz], label=label[iz])
+        end
+    end
+    GLMakie.axislegend(framevisible=false)
+    log10 && GLMakie.ylims!(3*maximum(Iω)*log10min, 3*maximum(Iω))
+    GLMakie.xlims!(speclims...)
+    GLMakie.DataInspector()
+    sfig
+end
+
+function _plot_slice_mm(x, y, z, modestrs, xlabel, ylabel, log10=false; fwlabel=false)
+    pfig = newfig()
+    scale = (log10 ? Base.log10 : identity)
+    GLMakie.Axis(pfig[1, 1], yscale = scale, xlabel=xlabel, ylabel=ylabel)
+    for sidx = 1:size(y, 3) # iterate over z-slices
+        zs = @sprintf("%.2f cm", z[sidx]*100)
+        label = "$zs ($(modestrs[1]))"
+        if fwlabel
+            fw = Maths.fwhm(x, y[:, 1, sidx])
+            label *= @sprintf(" [%.2f %s]", fw, "fs")
+        end
+        line = GLMakie.lines!(x, y[:, 1, sidx]; label)
+        for midx = 2:size(y, 2) # iterate over modes
+            label = "$zs ($(modestrs[midx]))"
+            if fwlabel
+                fw = Maths.fwhm(x, y[:, midx, sidx])
+                label *= @sprintf(" [%.2f %s]", fw, "fs")
+            end
+            GLMakie.lines!(x, y[:, midx, sidx]; label,
+                   color=line[:color], cycle=[:linestyle])
+        end
+    end
+    pfig
+end
+
+function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
+                     trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40,
+                     surface3d=false,
+                     kwargs...)
+    ω = Maths.rfftfreq(t)[2:end]
+    tmin, tmax = extrema(trange)
+    tg = collect(range(tmin, tmax, length=N))
+    g = Maths.gabor(t, real(Et), tg, fw)
+    g = g[2:end, :]
+
+    speclims, speclabel, specyfac = getspeclims(λrange, specaxis)
+    specy, Ig = getIω(ω, g*Maths.rfftnorm(t[2]-t[1]), specaxis,
+                      specrange=speclims./specyfac)
+    
+    Ig = Maths.normbymax(Ig)
+    log && (Ig = 10*log10.(Ig))
+    clims = (log ? (dBmin, 0) : extrema(Ig))
+
+    if haskey(kwargs, :cmap)
+        cmap = kwargs[:cmap]
+    else
+        cmap = :viridis
+    end
+
+    fig = newfig()
+    if surface3d
+        ax, pl = GLMakie.surface(fig[1,1], tg.*1e15, specyfac*specy, Ig',
+                         colorrange=clims, colormap=cmap,
+                         axis=(;type=GLMakie.Axis3, azimuth = pi/4, elevation=pi/4,
+                         protrusions=75, perspectiveness=0.0, viewmode=:stretch,
+                         xlabel="Time (fs)", ylabel=speclabel, ylabeloffset=80,
+                         xlabeloffset=80, zgridvisible=false, zlabelvisible=false,
+                         zticksvisible=false, zticklabelsvisible=false,
+                         yzpanelvisible=false, xzpanelvisible=false,
+                         ygridvisible=false, xgridvisible=false,
+                         zspinesvisible=false, zautolimitmargin=(0,0),
+                         xautolimitmargin=(0.0,0.0), yautolimitmargin=(0,0),
+                         xspinesvisible=false, yspinesvisible=false))
+    else
+        ax, pl = GLMakie.heatmap(fig[1,1], tg.*1e15, specyfac*specy, Ig',
+                                colorrange=clims, interpolate=true, colormap=cmap,
+                                axis=(; xlabel="Time (fs)", ylabel=speclabel))
+        GLMakie.ylims!(ax, speclims)
+    end
+    GLMakie.Colorbar(fig[1, 2], pl)
+    GLMakie.DataInspector()
+    fig
+end
+
+function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))
+    e = Processing.energy(output; bandpass=bandpass)
+    eall = Processing.energy(output)
+
+    multimode, modestrs = get_modes(output)
+    if multimode
+        e0 = sum(eall[:, 1])
+        modes = isnothing(modes) ? (1:size(e, 1)) : modes
+        if modes == :sum
+            e = dropdims(sum(e, dims=1), dims=1)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            e = e[modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    else
+        e0 = eall[1]
+    end
+
+    z = output["z"]*100
+    println("$(size(e'))")
+
+    fig = newfig()
+    ax = GLMakie.Axis(fig[1, 1], xlabel="Distance (cm)", ylabel="Energy (μJ)")
+    rax = GLMakie.Axis(fig[1, 1], yaxisposition = :right, ylabel="Conversion efficiency (%)")
+    GLMakie.hidespines!(rax)
+    GLMakie.hidexdecorations!(rax)
+    for i in 1:size(e')[1]
+        GLMakie.lines!(ax, z, 1e6*e'[i,:])
+    end
+
+    maxe = maximum(1e6*e)
+    GLMakie.xlims!(ax, extrema(z)...)
+    GLMakie.ylims!(ax, 0, maxe)
+    GLMakie.ylims!(rax, 0, 100*maxe/1e6/e0)
+    GLMakie.xlims!(rax, extrema(z)...)
+    GLMakie.DataInspector()
+    fig
+end
+
+function add_fwhm_legends(ax, unit)
+    leg = ax.get_legend()
+    texts = leg.get_texts()
+    handles, labels = ax.get_legend_handles_labels()
+    
+    for (ii, line) in enumerate(handles)
+        xy = line.get_xydata()
+        fw = Maths.fwhm(xy[:, 1], xy[:, 2])
+        t = texts[ii]
+        s = t.get_text()
+        s *= @sprintf(" [%.2f %s]", fw, unit)
+        t.set_text(s)
+    end
+end
+
+"""
+    cornertext(ax, text;
+               corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
+
+Place a `text` in the axes `ax` in the corner defined by `corner`. Padding can be
+defined for `x` and `y` together via `pad` or separately via `xpad` and `ypad`. Further
+keyword arguments are passed to `plt.text`. 
+
+Possible values for `corner` are `ul`, `ur`, `ll`, `lr` where the first letter
+defines upper/lower and the second defines left/right.
+"""
+function cornertext(ax, text; corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
+    xpad = isnothing(xpad) ? pad : xpad
+    ypad = isnothing(ypad) ? pad : ypad
+    if corner[1] == 'u'
+        val = "top"
+        y = 1 - ypad
+    elseif corner[1] == 'l'
+        val = "bottom"
+        y = ypad
+    else
+        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
+    end
+    if corner[2] == 'l'
+        hal = "left"
+        x = xpad
+    elseif corner[2] == 'r'
+        hal = "right"
+        x = 1 - xpad
+    else
+        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
+    end
+    ax.text(x, y, text; horizontalalignment=hal, verticalalignment=val,
+                 transform=ax.transAxes, kwargs...)
+end
+
+end # module
\ No newline at end of file
diff --git a/ext/PyPlotExt.jl b/ext/PyPlotExt.jl
new file mode 100644
index 00000000..39da864a
--- /dev/null
+++ b/ext/PyPlotExt.jl
@@ -0,0 +1,663 @@
+module PyPlotExt
+import Luna: Grid, Maths, PhysData, Processing
+import Luna.PhysData: wlfreq, c, ε_0
+import Luna.Output: AbstractOutput
+import Luna.Processing: makegrid, getIω, getEω, getEt, nearest_z
+import PyPlot: ColorMap, plt, pygui, Figure
+import FFTW
+import Printf: @sprintf
+import Base: display
+
+"""
+    displayall()
+
+`display` all currently open PyPlot figures.
+"""
+function displayall()
+    for fign in plt.get_fignums()
+        fig = plt.figure(fign)
+        display(fig)
+    end
+
+end
+
+display(figs::AbstractArray{Figure, N}) where N = [display(fig) for fig in figs]
+
+"""
+    cmap_white(cmap, N=512, n=8)
+
+Replace the lowest colour stop of `cmap` (after splitting into `n` stops) with white and
+create a new colourmap with `N` stops.
+"""
+function cmap_white(cmap; N=2^12, n=8)
+    vals = collect(range(0, 1, length=n))
+    vals_i = collect(range(0, 1, length=N))
+    cm = ColorMap(cmap)
+    clist = cm(vals)
+    clist[1, :] = [1, 1, 1, 1]
+    clist_i = Array{Float64}(undef, (N, 4))
+    for ii in 1:4
+        clist_i[:, ii] .= Maths.BSpline(vals, clist[:, ii]).(vals_i)
+    end
+    ColorMap(clist_i)
+end
+
+"""
+    cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
+
+Make an array of `num` different colours that follow the colourmap `cmap` between the values
+`cmin` and `cmax`.
+"""
+function cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
+    cm = ColorMap(cmap)
+    n = collect(range(cmin, cmax; length=num))
+    cm.(n)
+end
+
+"""
+    subplotgrid(N, portrait=true, kwargs...)
+
+Create a figure with `N` subplots laid out in a grid that is as close to square as possible.
+If `portrait` is `true`, try to lay out the grid in portrait orientation (taller than wide),
+otherwise landscape (wider than tall).
+"""
+function subplotgrid(N, portrait=true; colw=4, rowh=2.5, title=nothing)
+    cols = ceil(Int, sqrt(N))
+    rows = ceil(Int, N/cols)
+    portrait && ((rows, cols) = (cols, rows))
+    fig, axs = plt.subplots(rows, cols, num=title)
+    ndims(axs) > 1 && (axs = permutedims(axs, (2, 1)))
+    if cols*rows > N
+        for axi in axs[N+1:end]
+            axi.remove()
+        end
+    end
+    fig.set_size_inches(cols*colw, rows*rowh)
+    fig, N > 1 ? axs : [axs]
+end
+
+"""
+    get_modes(output)
+
+Determine whether `output` contains a multimode simulation, and if so, return the names
+of the modes.
+"""
+function get_modes(output)
+    t = output["simulation_type"]["transform"]
+    !startswith(t, "TransModal") && return false, nothing
+    lines = split(t, "\n")
+    modeline = findfirst(li -> startswith(li, "  modes:"), lines)
+    endline = findnext(li -> !startswith(li, " "^4), lines, modeline+1)
+    mlines = lines[modeline+1 : endline-1]
+    labels = [match(r"{([^,]*),", li).captures[1] for li in mlines]
+    angles = zeros(length(mlines))
+    for (ii, li) in enumerate(mlines)
+        m = match(r"ϕ=(-?[0-9]+.[0-9]+)π", li)
+        isnothing(m) && continue # no angle information in mode label)
+        angles[ii] = parse(Float64, m.captures[1])
+    end
+    if !all(angles .== 0)
+        for i in eachindex(labels)
+            if startswith(labels[i], "HE")
+                if angles[i] == 0
+                    θs = "x"
+                elseif angles[i] == 0.5
+                    θs = "y"
+                else
+                    θs = "$(angles[i])π"
+                end
+                labels[i] *= " ($θs)"
+            end
+        end
+    end
+    return true, labels
+end
+
+"""
+    stats(output; kwargs...)
+
+Plot all statistics available in `output`. Additional `kwargs` are passed onto `plt.plot()`
+"""
+function stats(output; kwargs...)
+    stats = output["stats"]
+
+    pstats = [] # pulse statistics
+    haskey(stats, "energy") && push!(pstats, (1e6*stats["energy"], "Energy (μJ)"))
+    for (k, v) in pairs(stats)
+        startswith(k, "energy_") || continue
+        str = "Energy "*replace(k[8:end], "_" => " ")*" (μJ)"
+        push!(pstats, (1e6*stats[k], str))
+    end
+    for (k, v) in pairs(stats)
+        startswith(k, "peakpower_") || continue
+        Pfac, unit = power_unit(stats[k])
+        str = "Peak power "*replace(k[11:end], "_" => " ")*" ($unit)"
+        push!(pstats, (Pfac*stats[k], str))
+    end
+    if haskey(stats, "peakpower")
+        Pfac, unit = power_unit(stats["peakpower"])
+        push!(pstats, (Pfac*stats["peakpower"], "Peak power ($unit)"))
+    end
+    haskey(stats, "peakintensity") && push!(
+        pstats, (1e-16*stats["peakintensity"], "Peak Intensity (TW/cm\$^2\$)"))
+    haskey(stats, "fwhm_t_min") && push!(pstats, (1e15*stats["fwhm_t_min"], "min FWHM (fs)"))
+    haskey(stats, "fwhm_t_max") && push!(pstats, (1e15*stats["fwhm_t_max"], "max FWHM (fs)"))
+    haskey(stats, "fwhm_r") && push!(pstats, (1e6*stats["fwhm_r"], "Radial FWHM (μm)"))
+    haskey(stats, "ω0") && push!(pstats, (1e9*wlfreq.(stats["ω0"]), "Central wavelength (nm)"))
+
+    fstats = [] # fibre/waveguide/propagation statistics
+    if haskey(stats, "electrondensity")
+        push!(fstats, (1e-6*stats["electrondensity"], "Electron density (cm\$^{-3}\$)"))
+        if haskey(stats, "density")
+            push!(fstats,
+                 (100*stats["electrondensity"]./stats["density"], "Ionisation fraction (%)"))
+        end
+    end
+    haskey(stats, "density") && push!(
+        fstats, (1e-6*stats["density"], "Density (cm\$^{-3}\$)"))
+    haskey(stats, "pressure") && push!(
+        fstats, (stats["pressure"], "Pressure (bar)"))
+    haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)"))
+    haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)"))
+    haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)"))
+
+    z = stats["z"]*1e2
+
+    multimode, modes = get_modes(output)
+
+    Npl = length(pstats)
+    if Npl > 0
+        pfig, axs = subplotgrid(Npl, title="Pulse stats")
+        for n in 1:Npl
+            ax = axs[n]
+            data, label = pstats[n]
+            multimode && (ndims(data) > 1) && (data = data')
+            ax.plot(z, data; kwargs...)
+            ax.set_xlabel("Distance (cm)")
+            ax.set_ylabel(label)
+            multimode && (ndims(data) > 1) && ax.semilogy()
+            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
+        end
+        pfig.tight_layout()
+    end
+    
+    Npl = length(fstats)
+    if Npl > 0
+        ffig, axs = subplotgrid(Npl, title="Other stats")
+        for n in 1:Npl
+            ax = axs[n]
+            data, label = fstats[n]
+            multimode && (ndims(data) > 1) && (data = data')
+            ax.plot(z, data; kwargs...)
+            ax.set_xlabel("Distance (cm)")
+            ax.set_ylabel(label)
+            multimode && (ndims(data) > 1) && should_log10(data) && ax.semilogy()
+            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
+        end
+        ffig.tight_layout()
+    end
+    [pfig, ffig]
+end
+
+"""
+    should_log10(A, tolfac=10)
+
+For multi-line plots, determine whether data for different lines contained in `A` spans
+a sufficiently large range that a logarithmic scale should be used. By default, this is the
+case when there is any point where the lines are different by more than a factor of 10.
+"""
+function should_log10(A, tolfac=10)
+    mi = minimum(A; dims=2)
+    ma = maximum(A; dims=2)
+    any(ma./mi .> 10)
+end
+
+window_str(::Nothing) = ""
+window_str(win::NTuple{4, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win[2:3]...)
+window_str(win::NTuple{2, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win...)
+window_str(window) = "custom bandpass"
+
+"""
+    prop_2D(output, specaxis=:f)
+
+Make false-colour propagation plots for `output`, using spectral x-axis `specaxis` (see
+[`getIω`](@ref)). For multimode simulations, create one figure for each mode plus one for
+the sum of all modes.
+
+# Keyword arguments
+- `λrange::Tuple(Float64, Float64)` : x-axis limits for spectral plot (wavelength in metres)
+- `trange::Tuple(Float64, Float64)` : x-axis limits for time-domain plot (time in seconds)
+- `dBmin::Float64` : lower colour-scale limit for logarithmic spectral plot
+- `resolution::Real` smooth the spectral energy density as defined by [`getIω`](@ref).
+"""
+function prop_2D(output, specaxis=:f;
+                 trange=(-50e-15, 50e-15), bandpass=nothing,
+                 λrange=(150e-9, 2000e-9), dBmin=-60,
+                 resolution=nothing, modes=nothing, oversampling=4,
+                 kwargs...)
+    z = output["z"]*1e2
+    if specaxis == :λ
+            specx, Iω = getIω(output, specaxis, specrange=λrange, resolution=resolution)
+    else
+            specx, Iω = getIω(output, specaxis, resolution=resolution)
+    end
+
+    t, Et = getEt(output; trange, bandpass, oversampling)
+    It = abs2.(Et)
+
+    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
+    specx .*= specxfac
+
+    multimode, modelabels = get_modes(output)
+
+    if multimode
+        fig = _prop2D_mm(modelabels, modeidcs(modes, modelabels), t, z, specx, It, Iω,
+                         speclabel, speclims, trange, dBmin, window_str(bandpass);
+                         kwargs...)
+    else
+        fig = _prop2D_sm(t, z, specx, It, Iω,
+                         speclabel, speclims, trange, dBmin, window_str(bandpass);
+                         kwargs...)
+    end
+    fig
+end
+
+modeidcs(m::Int, ml) = [m]
+modeidcs(m::Symbol, ml) = (m == :sum) ? [] : error("modes must be :sum, a single integer, or iterable")
+modeidcs(m::Nothing, ml) = 1:length(ml)
+modeidcs(m, ml) = m
+
+# Helper function to convert λrange to the correct numbers depending on specaxis
+function getspeclims(λrange, specaxis)
+    if specaxis == :f
+        specxfac = 1e-15
+        speclims = (specxfac*c/maximum(λrange), specxfac*c/minimum(λrange))
+        speclabel = "Frequency (PHz)"
+    elseif specaxis == :ω
+        specxfac = 1e-15
+        speclims = (specxfac*wlfreq(maximum(λrange)), specxfac*wlfreq(minimum(λrange)))
+        speclabel = "Angular frequency (rad/fs)"
+    elseif specaxis == :λ
+        specxfac = 1e9
+        speclims = λrange .* specxfac
+        speclabel = "Wavelength (nm)"
+    else
+        error("Unknown specaxis $specaxis")
+    end
+    return speclims, speclabel, specxfac
+end
+
+# single-mode 2D propagation plots
+function _prop2D_sm(t, z, specx, It, Iω, speclabel, speclims, trange, dBmin, bpstr; kwargs...)
+    id = "($(string(hash(gensym()); base=16)[1:4])) "
+    num = id * "Propagation" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+    pfig, axs = plt.subplots(1, 2, num=num)
+    pfig.set_size_inches(12, 4)
+    Iω = Maths.normbymax(Iω)
+    _spec2D_log(axs[1], specx, z, Iω, dBmin, speclabel, speclims; kwargs...)
+
+    _time2D(axs[2], t, z, It, trange; kwargs...)
+    pfig.tight_layout()
+    return pfig
+end
+
+# multi-mode 2D propagation plots
+function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
+                    speclabel, speclims, trange, dBmin, bpstr;
+                    kwargs...)
+    pfigs = Figure[]
+    Iω = Maths.normbymax(Iω)
+    id = "($(string(hash(gensym()); base=16)[1:4])) "
+    for mi in modes
+        num = id * "Propagation ($(modelabels[mi]))" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+        pfig, axs = plt.subplots(1, 2, num=num)
+        pfig.set_size_inches(12, 4)
+        _spec2D_log(axs[1], specx, z, Iω[:, mi, :], dBmin, speclabel, speclims; kwargs...)
+
+        _time2D(axs[2], t, z, It[:, mi, :], trange; kwargs...)
+        push!(pfigs, pfig)
+    end
+
+    num = id * "Propagation (all modes)" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+    pfig, axs = plt.subplots(1, 2, num=num)
+    pfig.set_size_inches(12, 4)
+    Iωall = dropdims(sum(Iω, dims=2), dims=2)
+    _spec2D_log(axs[1], specx, z, Iωall, dBmin, speclabel, speclims; kwargs...)
+
+    Itall = dropdims(sum(It, dims=2), dims=2)
+    _time2D(axs[2], t, z, Itall, trange; kwargs...)
+    pfig.tight_layout()
+    push!(pfigs, pfig)
+
+    return pfigs
+end
+
+# a single logarithmic colour-scale spectral domain plot
+function _spec2D_log(ax, specx, z, I, dBmin, speclabel, speclims; kwargs...)
+    im = ax.pcolormesh(specx, z, 10*log10.(transpose(I)); shading="auto", kwargs...)
+    im.set_clim(dBmin, 0)
+    cb = plt.colorbar(im, ax=ax)
+    cb.set_label("SED (dB)")
+    ax.set_ylabel("Distance (cm)")
+    ax.set_xlabel(speclabel)
+    ax.set_xlim(speclims...)
+end
+
+# a single time-domain propagation plot
+function _time2D(ax, t, z, I, trange; kwargs...)
+    Pfac, unit = power_unit(I)
+    im = ax.pcolormesh(t*1e15, z, Pfac*transpose(I); shading="auto", kwargs...)
+    cb = plt.colorbar(im, ax=ax)
+    cb.set_label("Power ($unit)")
+    ax.set_xlim(trange.*1e15)
+    ax.set_xlabel("Time (fs)")
+    ax.set_ylabel("Distance (cm)")
+end
+
+"""
+    time_1D(output, zslice, y=:Pt, kwargs...)
+
+Create lineplots of time-domain slice(s) of the propagation.
+
+The keyword argument `y` determines
+what is plotted: `:Pt` (power, default), `:Esq` (squared electric field) or `:Et` (electric field).
+
+The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
+a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
+
+The keyword argument `oversampling` determines the amount of oversampling done before plotting.
+
+Other `kwargs` are passed onto `plt.plot`.
+"""
+function time_1D(output, zslice=maximum(output["z"]);
+                y=:Pt, modes=nothing,
+                oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing,
+                FTL=false, propagate=nothing,
+                kwargs...)
+    t, Et, zactual = getEt(output, zslice,
+                           trange=trange, oversampling=oversampling, bandpass=bandpass,
+                           FTL=FTL, propagate=propagate)
+    if y == :Pt
+        yt = abs2.(Et)
+    elseif y == :Et
+        yt = real(Et)
+    elseif y == :Esq
+        yt = real(Et).^2
+    else
+        error("unknown time plot variable $y")
+    end
+    multimode, modestrs = get_modes(output)
+    if multimode
+        if modes == :sum
+            y == :Pt || error("Modal sum can only be plotted for power!")
+            yt = dropdims(sum(yt, dims=2), dims=2)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            yt = yt[:, modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    end
+
+    yfac, unit = power_unit(abs2.(Et), y)
+
+    sfig = plt.figure()
+    if multimode && nmodes > 1
+        _plot_slice_mm(plt.gca(), t*1e15, yfac*yt, zactual, modestrs; kwargs...)
+    else
+        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
+        label = multimode ? zs.*" ($modestrs)" : zs
+        for iz in eachindex(zactual)
+            plt.plot(t*1e15, yfac*yt[:, iz]; label=label[iz], kwargs...)
+        end
+    end
+    plt.legend(frameon=false)
+    add_fwhm_legends(plt.gca(), "fs")
+    plt.xlabel("Time (fs)")
+    plt.xlim(1e15.*trange)
+    ylab = y == :Et ?  "Field ($unit)" : "Power ($unit)"
+    plt.ylabel(ylab)
+    y == :Et || plt.ylim(ymin=0)
+    sfig.set_size_inches(8.5, 5)
+    sfig.tight_layout()
+    sfig
+end
+
+# Automatically find power unit depending on scale of electric field.
+function power_unit(Pt, y=:Pt)
+    units = ["kW", "MW", "GW", "TW", "PW"]
+    Pmax = maximum(Pt)
+    oom = clamp(floor(Int, log10(Pmax)/3), 1, 5) # maximum unit is PW
+    powerfac = 1/10^(oom*3)
+    if y == :Et
+        sqrt(powerfac), "$(units[oom])\$^{1/2}\$"
+    else
+        return powerfac, units[oom]
+    end
+end    
+
+"""
+    spec_1D(output, zslice, specaxis=:λ, log10=true, log10min=1e-6)
+
+Create lineplots of spectral-domain slices of the propagation.
+
+The x-axis is determined by `specaxis` (see [`getIω`](@ref)).
+
+If `log10` is true, plot on a logarithmic scale, with a y-axis range of `log10min`. 
+
+The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
+a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
+
+Other `kwargs` are passed onto `plt.plot`.
+"""
+function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
+                 modes=nothing, λrange=(150e-9, 1200e-9),
+                 log10=true, log10min=1e-6, resolution=nothing,
+                 kwargs...)
+    if specaxis == :λ
+        specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution)
+    else
+        specx, Iω, zactual = getIω(output, specaxis, zslice, resolution=resolution)
+    end
+    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
+    multimode, modestrs = get_modes(output)
+    if multimode
+        modes = isnothing(modes) ? (1:size(Iω, 2)) : modes
+        if modes == :sum
+            Iω = dropdims(sum(Iω, dims=2), dims=2)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            Iω = Iω[:, modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    end
+
+    specx .*= specxfac
+
+    sfig = plt.figure()
+    if multimode && nmodes > 1
+        _plot_slice_mm(plt.gca(), specx, Iω, zactual, modestrs, log10; kwargs...)
+    else
+        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
+        label = multimode ? zs.*" ($modestrs)" : zs
+        for iz in eachindex(zactual)
+            (log10 ? plt.semilogy : plt.plot)(specx, Iω[:, iz]; label=label[iz], kwargs...)
+        end
+    end
+    plt.legend(frameon=false)
+    plt.xlabel(speclabel)
+    plt.ylabel("Spectral energy density")
+    log10 && plt.ylim(3*maximum(Iω)*log10min, 3*maximum(Iω))
+    plt.xlim(speclims...)
+    sfig.set_size_inches(8.5, 5)
+    sfig.tight_layout()
+    sfig
+end
+
+dashes = [(0, (10, 1)),
+          (0, (5, 1)),
+          (0, (1, 0.5)),
+          (0, (1, 0.5, 1, 0.5, 3, 1)),
+          (0, (5, 1, 1, 1))]
+
+function _plot_slice_mm(ax, x, y, z, modestrs, log10=false, fwhm=false; kwargs...)
+    pfun = (log10 ? ax.semilogy : ax.plot)
+    for sidx = 1:size(y, 3) # iterate over z-slices
+        zs = @sprintf("%.2f cm", z[sidx]*100)
+        line = pfun(x, y[:, 1, sidx]; label="$zs ($(modestrs[1]))", kwargs...)[1]
+        for midx = 2:size(y, 2) # iterate over modes
+            pfun(x, y[:, midx, sidx], linestyle=dashes[midx], color=line.get_color(),
+                 label="$zs ($(modestrs[midx]))"; kwargs...)
+        end
+    end
+end
+
+function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
+                     trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40,
+                     kwargs...)
+    ω = Maths.rfftfreq(t)[2:end]
+    tmin, tmax = extrema(trange)
+    tg = collect(range(tmin, tmax, length=N))
+    g = Maths.gabor(t, real(Et), tg, fw)
+    g = g[2:end, :]
+
+    specy, Ig = getIω(ω, g*Maths.rfftnorm(t[2]-t[1]), specaxis)
+    speclims, speclabel, specyfac = getspeclims(λrange, specaxis)
+
+    log && (Ig = 10*log10.(Maths.normbymax(Ig)))
+
+    fig = plt.figure()
+    plt.pcolormesh(tg.*1e15, specyfac*specy, Ig; shading="auto", kwargs...)
+    plt.ylim(speclims...)
+    plt.ylabel(speclabel)
+    plt.xlabel("Time (fs)")
+    log && plt.clim(dBmin, 0)
+    plt.colorbar()
+    fig
+end
+
+function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))
+    e = Processing.energy(output; bandpass=bandpass)
+    eall = Processing.energy(output)
+
+    multimode, modestrs = get_modes(output)
+    if multimode
+        e0 = sum(eall[:, 1])
+        modes = isnothing(modes) ? (1:size(e, 1)) : modes
+        if modes == :sum
+            e = dropdims(sum(e, dims=1), dims=1)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            e = e[modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    else
+        e0 = eall[1]
+    end
+
+    z = output["z"]*100
+
+    fig = plt.figure()
+    ax = plt.axes()
+    ax.plot(z, 1e6*e')
+    ax.set_xlim(extrema(z)...)
+    ax.set_ylim(ymin=0)
+    ax.set_xlabel("Distance (cm)")
+    ax.set_ylabel("Energy (μJ)")
+    rax = ax.twinx()
+    rax.plot(z, 100*(e/e0)', linewidth=0)
+    lims = ax.get_ylim()
+    rax.set_ylim(100/(1e6*e0).*lims)
+    rax.set_ylabel("Conversion efficiency (%)")
+    fig.set_size_inches(figsize...)
+    fig
+end
+
+
+function auto_fwhm_arrows(ax, x, y; color="k", arrowlength=nothing, hpad=0, linewidth=1,
+                                    text=nothing, units="fs", kwargs...)
+    left, right = Maths.level_xings(x, y; kwargs...)
+    fw = abs(right - left)
+    halfmax = maximum(y)/2
+    arrowlength = isnothing(arrowlength) ? 2*fw : arrowlength
+
+    ax.annotate("", xy=(left-hpad, halfmax),
+                xytext=(left-hpad-arrowlength, halfmax),
+                arrowprops=Dict("arrowstyle" => "->",
+                                "color" => color,
+                                "linewidth" => linewidth))
+    ax.annotate("", xy=(right+hpad, halfmax),
+                xytext=(right+hpad+arrowlength, halfmax),
+                arrowprops=Dict("arrowstyle" => "->",
+                                "color" => color,
+                                "linewidth" => linewidth))
+
+    if text == :left
+        ax.text(left-arrowlength/2, 1.1*halfmax, @sprintf("%.2f %s", fw, units),
+                ha="right", color=color)
+    elseif text == :right
+        ax.text(right+arrowlength/2, 1.1*halfmax, @sprintf("%.2f %s", fw, units),
+                color=color)
+    end
+end
+
+function add_fwhm_legends(ax, unit)
+    leg = ax.get_legend()
+    texts = leg.get_texts()
+    handles, labels = ax.get_legend_handles_labels()
+    
+    for (ii, line) in enumerate(handles)
+        xy = line.get_xydata()
+        fw = Maths.fwhm(xy[:, 1], xy[:, 2])
+        t = texts[ii]
+        s = t.get_text()
+        s *= @sprintf(" [%.2f %s]", fw, unit)
+        t.set_text(s)
+    end
+end
+
+"""
+    cornertext(ax, text;
+               corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
+
+Place a `text` in the axes `ax` in the corner defined by `corner`. Padding can be
+defined for `x` and `y` together via `pad` or separately via `xpad` and `ypad`. Further
+keyword arguments are passed to `plt.text`. 
+
+Possible values for `corner` are `ul`, `ur`, `ll`, `lr` where the first letter
+defines upper/lower and the second defines left/right.
+"""
+function cornertext(ax, text; corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
+    xpad = isnothing(xpad) ? pad : xpad
+    ypad = isnothing(ypad) ? pad : ypad
+    if corner[1] == 'u'
+        val = "top"
+        y = 1 - ypad
+    elseif corner[1] == 'l'
+        val = "bottom"
+        y = ypad
+    else
+        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
+    end
+    if corner[2] == 'l'
+        hal = "left"
+        x = xpad
+    elseif corner[2] == 'r'
+        hal = "right"
+        x = 1 - xpad
+    else
+        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
+    end
+    ax.text(x, y, text; horizontalalignment=hal, verticalalignment=val,
+                 transform=ax.transAxes, kwargs...)
+end
+
+end # module
\ No newline at end of file
diff --git a/ext/PythonPlotExt.jl b/ext/PythonPlotExt.jl
new file mode 100644
index 00000000..97168ad9
--- /dev/null
+++ b/ext/PythonPlotExt.jl
@@ -0,0 +1,669 @@
+module PythonPlotExt
+import Luna: Grid, Maths, PhysData, Processing
+import Luna.PhysData: wlfreq, c, ε_0
+import Luna.Output: AbstractOutput
+import Luna.Processing: makegrid, getIω, getEω, getEt, nearest_z
+import FFTW
+import Printf: @sprintf
+import Base: display
+import PythonPlot: ColorMap, pygui, Figure, pyplot
+import PythonCall: pyconvert
+
+"""
+    displayall()
+
+`display` all currently open PythonPlot figures.
+"""
+function displayall()
+    for fign in pyplot.get_fignums()
+        fig = pyplot.figure(fign)
+        display(fig)
+    end
+
+end
+
+display(figs::AbstractArray{Figure, N}) where N = [display(fig) for fig in figs]
+
+"""
+    cmap_white(cmap, N=512, n=8)
+
+Replace the lowest colour stop of `cmap` (after splitting into `n` stops) with white and
+create a new colourmap with `N` stops.
+"""
+function cmap_white(cmap; N=2^12, n=8)
+    vals = collect(range(0, 1, length=n))
+    vals_i = collect(range(0, 1, length=N))
+    cm = ColorMap(cmap)
+    clist = pyconvert(Array, cm(vals))
+    clist[1, :] = [1, 1, 1, 1]
+    clist_i = Array{Float64}(undef, (N, 3))
+    for ii in 1:3
+        clist_i[:, ii] .= Maths.BSpline(vals, clist[:, ii]).(vals_i)
+    end
+    ColorMap(clist_i)
+end
+
+"""
+    cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
+
+Make an array of `num` different colours that follow the colourmap `cmap` between the values
+`cmin` and `cmax`.
+"""
+function cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
+    cm = ColorMap(cmap)
+    n = collect(range(cmin, cmax; length=num))
+    cm.(n)
+end
+
+"""
+    subplotgrid(N, portrait=true, kwargs...)
+
+Create a figure with `N` subplots laid out in a grid that is as close to square as possible.
+If `portrait` is `true`, try to lay out the grid in portrait orientation (taller than wide),
+otherwise landscape (wider than tall).
+"""
+function subplotgrid(N, portrait=true; colw=4, rowh=2.5, title=nothing)
+    cols = ceil(Int, sqrt(N))
+    rows = ceil(Int, N/cols)
+    portrait && ((rows, cols) = (cols, rows))
+    fig, axs = pyplot.subplots(rows, cols, num=title)
+    axs = pyconvert(Any, axs)
+    ndims(axs) > 1 && (axs = permutedims(axs, (2, 1)))
+    if cols*rows > N
+        for axi in axs[N+1:end]
+            axi.remove()
+        end
+    end
+    fig.set_size_inches(cols*colw, rows*rowh)
+    fig, N > 1 ? axs : [axs]
+end
+
+"""
+    get_modes(output)
+
+Determine whether `output` contains a multimode simulation, and if so, return the names
+of the modes.
+"""
+function get_modes(output)
+    t = output["simulation_type"]["transform"]
+    !startswith(t, "TransModal") && return false, nothing
+    lines = split(t, "\n")
+    modeline = findfirst(li -> startswith(li, "  modes:"), lines)
+    endline = findnext(li -> !startswith(li, " "^4), lines, modeline+1)
+    mlines = lines[modeline+1 : endline-1]
+    labels = [match(r"{([^,]*),", li).captures[1] for li in mlines]
+    angles = zeros(length(mlines))
+    for (ii, li) in enumerate(mlines)
+        m = match(r"ϕ=(-?[0-9]+.[0-9]+)π", li)
+        isnothing(m) && continue # no angle information in mode label)
+        angles[ii] = parse(Float64, m.captures[1])
+    end
+    if !all(angles .== 0)
+        for i in eachindex(labels)
+            if startswith(labels[i], "HE")
+                if angles[i] == 0
+                    θs = "x"
+                elseif angles[i] == 0.5
+                    θs = "y"
+                else
+                    θs = "$(angles[i])π"
+                end
+                labels[i] *= " ($θs)"
+            end
+        end
+    end
+    return true, labels
+end
+
+"""
+    stats(output; kwargs...)
+
+Plot all statistics available in `output`. Additional `kwargs` are passed onto `pyplot.plot()`
+"""
+function stats(output; kwargs...)
+    stats = output["stats"]
+
+    pstats = [] # pulse statistics
+    haskey(stats, "energy") && push!(pstats, (1e6*stats["energy"], "Energy (μJ)"))
+    for (k, v) in pairs(stats)
+        startswith(k, "energy_") || continue
+        str = "Energy "*replace(k[8:end], "_" => " ")*" (μJ)"
+        push!(pstats, (1e6*stats[k], str))
+    end
+    for (k, v) in pairs(stats)
+        startswith(k, "peakpower_") || continue
+        Pfac, unit = power_unit(stats[k])
+        str = "Peak power "*replace(k[11:end], "_" => " ")*" ($unit)"
+        push!(pstats, (Pfac*stats[k], str))
+    end
+    if haskey(stats, "peakpower")
+        Pfac, unit = power_unit(stats["peakpower"])
+        push!(pstats, (Pfac*stats["peakpower"], "Peak power ($unit)"))
+    end
+    haskey(stats, "peakintensity") && push!(
+        pstats, (1e-16*stats["peakintensity"], "Peak Intensity (TW/cm\$^2\$)"))
+    haskey(stats, "fwhm_t_min") && push!(pstats, (1e15*stats["fwhm_t_min"], "min FWHM (fs)"))
+    haskey(stats, "fwhm_t_max") && push!(pstats, (1e15*stats["fwhm_t_max"], "max FWHM (fs)"))
+    haskey(stats, "fwhm_r") && push!(pstats, (1e6*stats["fwhm_r"], "Radial FWHM (μm)"))
+    haskey(stats, "ω0") && push!(pstats, (1e9*wlfreq.(stats["ω0"]), "Central wavelength (nm)"))
+
+    fstats = [] # fibre/waveguide/propagation statistics
+    if haskey(stats, "electrondensity")
+        push!(fstats, (1e-6*stats["electrondensity"], "Electron density (cm\$^{-3}\$)"))
+        if haskey(stats, "density")
+            push!(fstats,
+                 (100*stats["electrondensity"]./stats["density"], "Ionisation fraction (%)"))
+        end
+    end
+    haskey(stats, "density") && push!(
+        fstats, (1e-6*stats["density"], "Density (cm\$^{-3}\$)"))
+    haskey(stats, "pressure") && push!(
+        fstats, (stats["pressure"], "Pressure (bar)"))
+    haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)"))
+    haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)"))
+    haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)"))
+
+    z = stats["z"]*1e2
+
+    multimode, modes = get_modes(output)
+
+    Npl = length(pstats)
+    if Npl > 0
+        pfig, axs = subplotgrid(Npl, title="Pulse stats")
+        for n in 1:Npl
+            ax = axs[n]
+            data, label = pstats[n]
+            multimode && (ndims(data) > 1) && (data = data')
+            ax.plot(z, data; kwargs...)
+            ax.set_xlabel("Distance (cm)")
+            ax.set_ylabel(label)
+            multimode && (ndims(data) > 1) && ax.semilogy()
+            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
+        end
+        pfig.tight_layout()
+    end
+    
+    Npl = length(fstats)
+    if Npl > 0
+        ffig, axs = subplotgrid(Npl, title="Other stats")
+        for n in 1:Npl
+            ax = axs[n]
+            data, label = fstats[n]
+            multimode && (ndims(data) > 1) && (data = data')
+            ax.plot(z, data; kwargs...)
+            ax.set_xlabel("Distance (cm)")
+            ax.set_ylabel(label)
+            multimode && (ndims(data) > 1) && should_log10(data) && ax.semilogy()
+            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
+        end
+        ffig.tight_layout()
+    end
+    [pfig, ffig]
+end
+
+"""
+    should_log10(A, tolfac=10)
+
+For multi-line plots, determine whether data for different lines contained in `A` spans
+a sufficiently large range that a logarithmic scale should be used. By default, this is the
+case when there is any point where the lines are different by more than a factor of 10.
+"""
+function should_log10(A, tolfac=10)
+    mi = minimum(A; dims=2)
+    ma = maximum(A; dims=2)
+    any(ma./mi .> 10)
+end
+
+window_str(::Nothing) = ""
+window_str(win::NTuple{4, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win[2:3]...)
+window_str(win::NTuple{2, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win...)
+window_str(window) = "custom bandpass"
+
+"""
+    prop_2D(output, specaxis=:f)
+
+Make false-colour propagation plots for `output`, using spectral x-axis `specaxis` (see
+[`getIω`](@ref)). For multimode simulations, create one figure for each mode plus one for
+the sum of all modes.
+
+# Keyword arguments
+- `λrange::Tuple(Float64, Float64)` : x-axis limits for spectral plot (wavelength in metres)
+- `trange::Tuple(Float64, Float64)` : x-axis limits for time-domain plot (time in seconds)
+- `dBmin::Float64` : lower colour-scale limit for logarithmic spectral plot
+- `resolution::Real` smooth the spectral energy density as defined by [`getIω`](@ref).
+"""
+function prop_2D(output, specaxis=:f;
+                 trange=(-50e-15, 50e-15), bandpass=nothing,
+                 λrange=(150e-9, 2000e-9), dBmin=-60,
+                 resolution=nothing, modes=nothing, oversampling=4,
+                 kwargs...)
+    z = output["z"]*1e2
+    if specaxis == :λ
+            specx, Iω = getIω(output, specaxis, specrange=λrange, resolution=resolution)
+    else
+            specx, Iω = getIω(output, specaxis, resolution=resolution)
+    end
+
+    t, Et = getEt(output; trange, bandpass, oversampling)
+    It = abs2.(Et)
+
+    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
+    specx .*= specxfac
+
+    multimode, modelabels = get_modes(output)
+
+    if multimode
+        fig = _prop2D_mm(modelabels, modeidcs(modes, modelabels), t, z, specx, It, Iω,
+                         speclabel, speclims, trange, dBmin, window_str(bandpass);
+                         kwargs...)
+    else
+        fig = _prop2D_sm(t, z, specx, It, Iω,
+                         speclabel, speclims, trange, dBmin, window_str(bandpass);
+                         kwargs...)
+    end
+    fig
+end
+
+modeidcs(m::Int, ml) = [m]
+modeidcs(m::Symbol, ml) = (m == :sum) ? [] : error("modes must be :sum, a single integer, or iterable")
+modeidcs(m::Nothing, ml) = 1:length(ml)
+modeidcs(m, ml) = m
+
+# Helper function to convert λrange to the correct numbers depending on specaxis
+function getspeclims(λrange, specaxis)
+    if specaxis == :f
+        specxfac = 1e-15
+        speclims = (specxfac*c/maximum(λrange), specxfac*c/minimum(λrange))
+        speclabel = "Frequency (PHz)"
+    elseif specaxis == :ω
+        specxfac = 1e-15
+        speclims = (specxfac*wlfreq(maximum(λrange)), specxfac*wlfreq(minimum(λrange)))
+        speclabel = "Angular frequency (rad/fs)"
+    elseif specaxis == :λ
+        specxfac = 1e9
+        speclims = λrange .* specxfac
+        speclabel = "Wavelength (nm)"
+    else
+        error("Unknown specaxis $specaxis")
+    end
+    return speclims, speclabel, specxfac
+end
+
+# single-mode 2D propagation plots
+function _prop2D_sm(t, z, specx, It, Iω, speclabel, speclims, trange, dBmin, bpstr; kwargs...)
+    id = "($(string(hash(gensym()); base=16)[1:4])) "
+    num = id * "Propagation" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+    pfig, axs = pyplot.subplots(1, 2, num=num)
+    axs = pyconvert(Any, axs)
+    pfig.set_size_inches(12, 4)
+    Iω = Maths.normbymax(Iω)
+    _spec2D_log(axs[1], specx, z, Iω, dBmin, speclabel, speclims; kwargs...)
+
+    _time2D(axs[2], t, z, It, trange; kwargs...)
+    pfig.tight_layout()
+    return pfig
+end
+
+# multi-mode 2D propagation plots
+function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
+                    speclabel, speclims, trange, dBmin, bpstr;
+                    kwargs...)
+    pfigs = []
+    Iω = Maths.normbymax(Iω)
+    id = "($(string(hash(gensym()); base=16)[1:4])) "
+    for mi in modes
+        num = id * "Propagation ($(modelabels[mi]))" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+        pfig, axs = pyplot.subplots(1, 2, num=num)
+        axs = pyconvert(Any, axs)
+        pfig.set_size_inches(12, 4)
+        _spec2D_log(axs[1], specx, z, Iω[:, mi, :], dBmin, speclabel, speclims; kwargs...)
+
+        _time2D(axs[2], t, z, It[:, mi, :], trange; kwargs...)
+        push!(pfigs, pfig)
+    end
+
+    num = id * "Propagation (all modes)" * ((length(bpstr) > 0) ? ", $bpstr" : "")
+    pfig, axs = pyplot.subplots(1, 2, num=num)
+    axs = pyconvert(Any, axs)
+    pfig.set_size_inches(12, 4)
+    Iωall = dropdims(sum(Iω, dims=2), dims=2)
+    _spec2D_log(axs[1], specx, z, Iωall, dBmin, speclabel, speclims; kwargs...)
+
+    Itall = dropdims(sum(It, dims=2), dims=2)
+    _time2D(axs[2], t, z, Itall, trange; kwargs...)
+    pfig.tight_layout()
+    push!(pfigs, pfig)
+
+    return pfigs
+end
+
+# a single logarithmic colour-scale spectral domain plot
+function _spec2D_log(ax, specx, z, I, dBmin, speclabel, speclims; kwargs...)
+    im = ax.pcolormesh(specx, z, 10*log10.(transpose(I)); shading="auto", kwargs...)
+    im.set_clim(dBmin, 0)
+    cb = pyplot.colorbar(im, ax=ax)
+    cb.set_label("SED (dB)")
+    ax.set_ylabel("Distance (cm)")
+    ax.set_xlabel(speclabel)
+    ax.set_xlim(speclims...)
+end
+
+# a single time-domain propagation plot
+function _time2D(ax, t, z, I, trange; kwargs...)
+    Pfac, unit = power_unit(I)
+    im = ax.pcolormesh(t*1e15, z, Pfac*transpose(I); shading="auto", kwargs...)
+    cb = pyplot.colorbar(im, ax=ax)
+    cb.set_label("Power ($unit)")
+    ax.set_xlim(trange.*1e15)
+    ax.set_xlabel("Time (fs)")
+    ax.set_ylabel("Distance (cm)")
+end
+
+"""
+    time_1D(output, zslice, y=:Pt, kwargs...)
+
+Create lineplots of time-domain slice(s) of the propagation.
+
+The keyword argument `y` determines
+what is plotted: `:Pt` (power, default), `:Esq` (squared electric field) or `:Et` (electric field).
+
+The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
+a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
+
+The keyword argument `oversampling` determines the amount of oversampling done before plotting.
+
+Other `kwargs` are passed onto `pyplot.plot`.
+"""
+function time_1D(output, zslice=maximum(output["z"]);
+                y=:Pt, modes=nothing,
+                oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing,
+                FTL=false, propagate=nothing,
+                kwargs...)
+    t, Et, zactual = getEt(output, zslice,
+                           trange=trange, oversampling=oversampling, bandpass=bandpass,
+                           FTL=FTL, propagate=propagate)
+    if y == :Pt
+        yt = abs2.(Et)
+    elseif y == :Et
+        yt = real(Et)
+    elseif y == :Esq
+        yt = real(Et).^2
+    else
+        error("unknown time plot variable $y")
+    end
+    multimode, modestrs = get_modes(output)
+    if multimode
+        if modes == :sum
+            y == :Pt || error("Modal sum can only be plotted for power!")
+            yt = dropdims(sum(yt, dims=2), dims=2)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            yt = yt[:, modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    end
+
+    yfac, unit = power_unit(abs2.(Et), y)
+
+    sfig = pyplot.figure()
+    if multimode && nmodes > 1
+        _plot_slice_mm(pyplot.gca(), t*1e15, yfac*yt, zactual, modestrs; kwargs...)
+    else
+        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
+        label = multimode ? zs.*" ($modestrs)" : zs
+        for iz in eachindex(zactual)
+            pyplot.plot(t*1e15, yfac*yt[:, iz]; label=label[iz], kwargs...)
+        end
+    end
+    pyplot.legend(frameon=false)
+    add_fwhm_legends(pyplot.gca(), "fs")
+    pyplot.xlabel("Time (fs)")
+    pyplot.xlim(1e15.*trange)
+    ylab = y == :Et ?  "Field ($unit)" : "Power ($unit)"
+    pyplot.ylabel(ylab)
+    y == :Et || pyplot.ylim(ymin=0)
+    sfig.set_size_inches(8.5, 5)
+    sfig.tight_layout()
+    sfig
+end
+
+# Automatically find power unit depending on scale of electric field.
+function power_unit(Pt, y=:Pt)
+    units = ["kW", "MW", "GW", "TW", "PW"]
+    Pmax = maximum(Pt)
+    oom = clamp(floor(Int, log10(Pmax)/3), 1, 5) # maximum unit is PW
+    powerfac = 1/10^(oom*3)
+    if y == :Et
+        sqrt(powerfac), "$(units[oom])\$^{1/2}\$"
+    else
+        return powerfac, units[oom]
+    end
+end    
+
+"""
+    spec_1D(output, zslice, specaxis=:λ, log10=true, log10min=1e-6)
+
+Create lineplots of spectral-domain slices of the propagation.
+
+The x-axis is determined by `specaxis` (see [`getIω`](@ref)).
+
+If `log10` is true, plot on a logarithmic scale, with a y-axis range of `log10min`. 
+
+The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
+a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
+
+Other `kwargs` are passed onto `pyplot.plot`.
+"""
+function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
+                 modes=nothing, λrange=(150e-9, 1200e-9),
+                 log10=true, log10min=1e-6, resolution=nothing,
+                 kwargs...)
+    if specaxis == :λ
+        specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution)
+    else
+        specx, Iω, zactual = getIω(output, specaxis, zslice, resolution=resolution)
+    end
+    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
+    multimode, modestrs = get_modes(output)
+    if multimode
+        modes = isnothing(modes) ? (1:size(Iω, 2)) : modes
+        if modes == :sum
+            Iω = dropdims(sum(Iω, dims=2), dims=2)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            Iω = Iω[:, modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    end
+
+    specx .*= specxfac
+
+    sfig = pyplot.figure()
+    if multimode && nmodes > 1
+        _plot_slice_mm(pyplot.gca(), specx, Iω, zactual, modestrs, log10; kwargs...)
+    else
+        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
+        label = multimode ? zs.*" ($modestrs)" : zs
+        for iz in eachindex(zactual)
+            (log10 ? pyplot.semilogy : pyplot.plot)(specx, Iω[:, iz]; label=label[iz], kwargs...)
+        end
+    end
+    pyplot.legend(frameon=false)
+    pyplot.xlabel(speclabel)
+    pyplot.ylabel("Spectral energy density")
+    log10 && pyplot.ylim(3*maximum(Iω)*log10min, 3*maximum(Iω))
+    pyplot.xlim(speclims...)
+    sfig.set_size_inches(8.5, 5)
+    sfig.tight_layout()
+    sfig
+end
+
+dashes = [(0, (10, 1)),
+          (0, (5, 1)),
+          (0, (1, 0.5)),
+          (0, (1, 0.5, 1, 0.5, 3, 1)),
+          (0, (5, 1, 1, 1))]
+
+function _plot_slice_mm(ax, x, y, z, modestrs, log10=false, fwhm=false; kwargs...)
+    pfun = (log10 ? ax.semilogy : ax.plot)
+    for sidx = 1:size(y, 3) # iterate over z-slices
+        zs = @sprintf("%.2f cm", z[sidx]*100)
+        line = pfun(x, y[:, 1, sidx]; label="$zs ($(modestrs[1]))", kwargs...)[0]
+        for midx = 2:size(y, 2) # iterate over modes
+            pfun(x, y[:, midx, sidx], linestyle=dashes[midx], color=line.get_color(),
+                 label="$zs ($(modestrs[midx]))"; kwargs...)
+        end
+    end
+end
+
+function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
+                     trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40,
+                     kwargs...)
+    ω = Maths.rfftfreq(t)[2:end]
+    tmin, tmax = extrema(trange)
+    tg = collect(range(tmin, tmax, length=N))
+    g = Maths.gabor(t, real(Et), tg, fw)
+    g = g[2:end, :]
+
+    specy, Ig = getIω(ω, g*Maths.rfftnorm(t[2]-t[1]), specaxis)
+    speclims, speclabel, specyfac = getspeclims(λrange, specaxis)
+
+    log && (Ig = 10*log10.(Maths.normbymax(Ig)))
+
+    fig = pyplot.figure()
+    pyplot.pcolormesh(tg.*1e15, specyfac*specy, Ig; shading="auto", kwargs...)
+    pyplot.ylim(speclims...)
+    pyplot.ylabel(speclabel)
+    pyplot.xlabel("Time (fs)")
+    log && pyplot.clim(dBmin, 0)
+    pyplot.colorbar()
+    fig
+end
+
+function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))
+    e = Processing.energy(output; bandpass=bandpass)
+    eall = Processing.energy(output)
+
+    multimode, modestrs = get_modes(output)
+    if multimode
+        e0 = sum(eall[:, 1])
+        modes = isnothing(modes) ? (1:size(e, 1)) : modes
+        if modes == :sum
+            e = dropdims(sum(e, dims=1), dims=1)
+            modestrs = join(modestrs, "+")
+            nmodes = 1
+        else
+            isnothing(modes) && (modes = 1:length(modestrs))
+            e = e[modes, :]
+            modestrs = modestrs[modes]
+            nmodes = length(modes)
+        end
+    else
+        e0 = eall[1]
+    end
+
+    z = output["z"]*100
+
+    fig = pyplot.figure()
+    ax = pyplot.axes()
+    ax.plot(z, 1e6*e')
+    ax.set_xlim(extrema(z)...)
+    ax.set_ylim(ymin=0)
+    ax.set_xlabel("Distance (cm)")
+    ax.set_ylabel("Energy (μJ)")
+    rax = ax.twinx()
+    rax.plot(z, 100*(e/e0)', linewidth=0)
+    lims = ax.get_ylim()
+    rax.set_ylim(100/(1e6*e0).*lims)
+    rax.set_ylabel("Conversion efficiency (%)")
+    fig.set_size_inches(figsize...)
+    fig
+end
+
+
+function auto_fwhm_arrows(ax, x, y; color="k", arrowlength=nothing, hpad=0, linewidth=1,
+                                    text=nothing, units="fs", kwargs...)
+    left, right = Maths.level_xings(x, y; kwargs...)
+    fw = abs(right - left)
+    halfmax = maximum(y)/2
+    arrowlength = isnothing(arrowlength) ? 2*fw : arrowlength
+
+    ax.annotate("", xy=(left-hpad, halfmax),
+                xytext=(left-hpad-arrowlength, halfmax),
+                arrowprops=Dict("arrowstyle" => "->",
+                                "color" => color,
+                                "linewidth" => linewidth))
+    ax.annotate("", xy=(right+hpad, halfmax),
+                xytext=(right+hpad+arrowlength, halfmax),
+                arrowprops=Dict("arrowstyle" => "->",
+                                "color" => color,
+                                "linewidth" => linewidth))
+
+    if text == :left
+        ax.text(left-arrowlength/2, 1.1*halfmax, @sprintf("%.2f %s", fw, units),
+                ha="right", color=color)
+    elseif text == :right
+        ax.text(right+arrowlength/2, 1.1*halfmax, @sprintf("%.2f %s", fw, units),
+                color=color)
+    end
+end
+
+function add_fwhm_legends(ax, unit)
+    leg = ax.get_legend()
+    texts = leg.get_texts()
+    handles, labels = ax.get_legend_handles_labels()
+    handles = pyconvert(Any, handles)
+    for (ii, line) in enumerate(handles)
+        xy = line.get_xydata()
+        xy = pyconvert(Any, xy)
+        fw = Maths.fwhm(xy[:, 1], xy[:, 2])
+        t = texts[ii-1]
+        s = pyconvert(Any, t.get_text())
+        s *= @sprintf(" [%.2f %s]", fw, unit)
+        t.set_text(s)
+    end
+end
+
+"""
+    cornertext(ax, text;
+               corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
+
+Place a `text` in the axes `ax` in the corner defined by `corner`. Padding can be
+defined for `x` and `y` together via `pad` or separately via `xpad` and `ypad`. Further
+keyword arguments are passed to `pyplot.text`. 
+
+Possible values for `corner` are `ul`, `ur`, `ll`, `lr` where the first letter
+defines upper/lower and the second defines left/right.
+"""
+function cornertext(ax, text; corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
+    xpad = isnothing(xpad) ? pad : xpad
+    ypad = isnothing(ypad) ? pad : ypad
+    if corner[1] == 'u'
+        val = "top"
+        y = 1 - ypad
+    elseif corner[1] == 'l'
+        val = "bottom"
+        y = ypad
+    else
+        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
+    end
+    if corner[2] == 'l'
+        hal = "left"
+        x = xpad
+    elseif corner[2] == 'r'
+        hal = "right"
+        x = 1 - xpad
+    else
+        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
+    end
+    ax.text(x, y, text; horizontalalignment=hal, verticalalignment=val,
+                 transform=ax.transAxes, kwargs...)
+end
+
+end # module
\ No newline at end of file
diff --git a/src/Plotting.jl b/src/Plotting.jl
index c73bcf33..c9ce7cec 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -3,26 +3,20 @@ import Luna: Grid, Maths, PhysData, Processing
 import Luna.PhysData: wlfreq, c, ε_0
 import Luna.Output: AbstractOutput
 import Luna.Processing: makegrid, getIω, getEω, getEt, nearest_z
-import PyPlot: ColorMap, plt, pygui, Figure
 import FFTW
 import Printf: @sprintf
 import Base: display
 
-"""
-    displayall()
-
-`display` all currently open PyPlot figures.
-"""
-function displayall()
-    for fign in plt.get_fignums()
-        fig = plt.figure(fign)
-        display(fig)
-    end
-
+function getext()
+    ext = Base.get_extension(@__MODULE__, :PythonPlotExt)
+    !isnothing(ext) && return ext
+    ext = Base.get_extension(@__MODULE__, :PyPlotExt)
+    !isnothing(ext) && return ext
+    ext = Base.get_extension(@__MODULE__, :GLMakieExt)
+    !isnothing(ext) && return ext
+    error("Please load a plotting backend: PythonPlot, PyPlot, or GLMakie.")
 end
 
-display(figs::AbstractArray{Figure, N}) where N = [display(fig) for fig in figs]
-
 """
     cmap_white(cmap, N=512, n=8)
 
@@ -30,193 +24,18 @@ Replace the lowest colour stop of `cmap` (after splitting into `n` stops) with w
 create a new colourmap with `N` stops.
 """
 function cmap_white(cmap; N=2^12, n=8)
-    vals = collect(range(0, 1, length=n))
-    vals_i = collect(range(0, 1, length=N))
-    cm = ColorMap(cmap)
-    clist = cm(vals)
-    clist[1, :] = [1, 1, 1, 1]
-    clist_i = Array{Float64}(undef, (N, 4))
-    for ii in 1:4
-        clist_i[:, ii] .= Maths.BSpline(vals, clist[:, ii]).(vals_i)
-    end
-    ColorMap(clist_i)
-end
-
-"""
-    cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
-
-Make an array of `num` different colours that follow the colourmap `cmap` between the values
-`cmin` and `cmax`.
-"""
-function cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
-    cm = ColorMap(cmap)
-    n = collect(range(cmin, cmax; length=num))
-    cm.(n)
-end
-
-"""
-    subplotgrid(N, portrait=true, kwargs...)
-
-Create a figure with `N` subplots laid out in a grid that is as close to square as possible.
-If `portrait` is `true`, try to lay out the grid in portrait orientation (taller than wide),
-otherwise landscape (wider than tall).
-"""
-function subplotgrid(N, portrait=true; colw=4, rowh=2.5, title=nothing)
-    cols = ceil(Int, sqrt(N))
-    rows = ceil(Int, N/cols)
-    portrait && ((rows, cols) = (cols, rows))
-    fig, axs = plt.subplots(rows, cols, num=title)
-    ndims(axs) > 1 && (axs = permutedims(axs, (2, 1)))
-    if cols*rows > N
-        for axi in axs[N+1:end]
-            axi.remove()
-        end
-    end
-    fig.set_size_inches(cols*colw, rows*rowh)
-    fig, N > 1 ? axs : [axs]
-end
-
-"""
-    get_modes(output)
-
-Determine whether `output` contains a multimode simulation, and if so, return the names
-of the modes.
-"""
-function get_modes(output)
-    t = output["simulation_type"]["transform"]
-    !startswith(t, "TransModal") && return false, nothing
-    lines = split(t, "\n")
-    modeline = findfirst(li -> startswith(li, "  modes:"), lines)
-    endline = findnext(li -> !startswith(li, " "^4), lines, modeline+1)
-    mlines = lines[modeline+1 : endline-1]
-    labels = [match(r"{([^,]*),", li).captures[1] for li in mlines]
-    angles = zeros(length(mlines))
-    for (ii, li) in enumerate(mlines)
-        m = match(r"ϕ=(-?[0-9]+.[0-9]+)π", li)
-        isnothing(m) && continue # no angle information in mode label)
-        angles[ii] = parse(Float64, m.captures[1])
-    end
-    if !all(angles .== 0)
-        for i in eachindex(labels)
-            if startswith(labels[i], "HE")
-                if angles[i] == 0
-                    θs = "x"
-                elseif angles[i] == 0.5
-                    θs = "y"
-                else
-                    θs = "$(angles[i])π"
-                end
-                labels[i] *= " ($θs)"
-            end
-        end
-    end
-    return true, labels
+    getext().cmap_white(cmap; N, n)
 end
 
 """
     stats(output; kwargs...)
 
-Plot all statistics available in `output`. Additional `kwargs` are passed onto `plt.plot()`
+Plot all statistics available in `output`. Additional `kwargs` are passed onto `pyplot.plot()`
 """
 function stats(output; kwargs...)
-    stats = output["stats"]
-
-    pstats = [] # pulse statistics
-    haskey(stats, "energy") && push!(pstats, (1e6*stats["energy"], "Energy (μJ)"))
-    for (k, v) in pairs(stats)
-        startswith(k, "energy_") || continue
-        str = "Energy "*replace(k[8:end], "_" => " ")*" (μJ)"
-        push!(pstats, (1e6*stats[k], str))
-    end
-    for (k, v) in pairs(stats)
-        startswith(k, "peakpower_") || continue
-        Pfac, unit = power_unit(stats[k])
-        str = "Peak power "*replace(k[11:end], "_" => " ")*" ($unit)"
-        push!(pstats, (Pfac*stats[k], str))
-    end
-    if haskey(stats, "peakpower")
-        Pfac, unit = power_unit(stats["peakpower"])
-        push!(pstats, (Pfac*stats["peakpower"], "Peak power ($unit)"))
-    end
-    haskey(stats, "peakintensity") && push!(
-        pstats, (1e-16*stats["peakintensity"], "Peak Intensity (TW/cm\$^2\$)"))
-    haskey(stats, "fwhm_t_min") && push!(pstats, (1e15*stats["fwhm_t_min"], "min FWHM (fs)"))
-    haskey(stats, "fwhm_t_max") && push!(pstats, (1e15*stats["fwhm_t_max"], "max FWHM (fs)"))
-    haskey(stats, "fwhm_r") && push!(pstats, (1e6*stats["fwhm_r"], "Radial FWHM (μm)"))
-    haskey(stats, "ω0") && push!(pstats, (1e9*wlfreq.(stats["ω0"]), "Central wavelength (nm)"))
-
-    fstats = [] # fibre/waveguide/propagation statistics
-    if haskey(stats, "electrondensity")
-        push!(fstats, (1e-6*stats["electrondensity"], "Electron density (cm\$^{-3}\$)"))
-        if haskey(stats, "density")
-            push!(fstats,
-                 (100*stats["electrondensity"]./stats["density"], "Ionisation fraction (%)"))
-        end
-    end
-    haskey(stats, "density") && push!(
-        fstats, (1e-6*stats["density"], "Density (cm\$^{-3}\$)"))
-    haskey(stats, "pressure") && push!(
-        fstats, (stats["pressure"], "Pressure (bar)"))
-    haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)"))
-    haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)"))
-    haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)"))
-
-    z = stats["z"]*1e2
-
-    multimode, modes = get_modes(output)
-
-    Npl = length(pstats)
-    if Npl > 0
-        pfig, axs = subplotgrid(Npl, title="Pulse stats")
-        for n in 1:Npl
-            ax = axs[n]
-            data, label = pstats[n]
-            multimode && (ndims(data) > 1) && (data = data')
-            ax.plot(z, data; kwargs...)
-            ax.set_xlabel("Distance (cm)")
-            ax.set_ylabel(label)
-            multimode && (ndims(data) > 1) && ax.semilogy()
-            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
-        end
-        pfig.tight_layout()
-    end
-    
-    Npl = length(fstats)
-    if Npl > 0
-        ffig, axs = subplotgrid(Npl, title="Other stats")
-        for n in 1:Npl
-            ax = axs[n]
-            data, label = fstats[n]
-            multimode && (ndims(data) > 1) && (data = data')
-            ax.plot(z, data; kwargs...)
-            ax.set_xlabel("Distance (cm)")
-            ax.set_ylabel(label)
-            multimode && (ndims(data) > 1) && should_log10(data) && ax.semilogy()
-            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
-        end
-        ffig.tight_layout()
-    end
-    [pfig, ffig]
+    getext().stats(output; kwargs...)
 end
 
-"""
-    should_log10(A, tolfac=10)
-
-For multi-line plots, determine whether data for different lines contained in `A` spans
-a sufficiently large range that a logarithmic scale should be used. By default, this is the
-case when there is any point where the lines are different by more than a factor of 10.
-"""
-function should_log10(A, tolfac=10)
-    mi = minimum(A; dims=2)
-    ma = maximum(A; dims=2)
-    any(ma./mi .> 10)
-end
-
-window_str(::Nothing) = ""
-window_str(win::NTuple{4, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win[2:3]...)
-window_str(win::NTuple{2, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win...)
-window_str(window) = "custom bandpass"
-
 """
     prop_2D(output, specaxis=:f)
 
@@ -235,123 +54,8 @@ function prop_2D(output, specaxis=:f;
                  λrange=(150e-9, 2000e-9), dBmin=-60,
                  resolution=nothing, modes=nothing, oversampling=4,
                  kwargs...)
-    z = output["z"]*1e2
-    if specaxis == :λ
-            specx, Iω = getIω(output, specaxis, specrange=λrange, resolution=resolution)
-    else
-            specx, Iω = getIω(output, specaxis, resolution=resolution)
-    end
-
-    t, Et = getEt(output; trange, bandpass, oversampling)
-    It = abs2.(Et)
-
-    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
-    specx .*= specxfac
-
-    multimode, modelabels = get_modes(output)
-
-    if multimode
-        fig = _prop2D_mm(modelabels, modeidcs(modes, modelabels), t, z, specx, It, Iω,
-                         speclabel, speclims, trange, dBmin, window_str(bandpass);
-                         kwargs...)
-    else
-        fig = _prop2D_sm(t, z, specx, It, Iω,
-                         speclabel, speclims, trange, dBmin, window_str(bandpass);
-                         kwargs...)
-    end
-    fig
-end
-
-modeidcs(m::Int, ml) = [m]
-modeidcs(m::Symbol, ml) = (m == :sum) ? [] : error("modes must be :sum, a single integer, or iterable")
-modeidcs(m::Nothing, ml) = 1:length(ml)
-modeidcs(m, ml) = m
-
-# Helper function to convert λrange to the correct numbers depending on specaxis
-function getspeclims(λrange, specaxis)
-    if specaxis == :f
-        specxfac = 1e-15
-        speclims = (specxfac*c/maximum(λrange), specxfac*c/minimum(λrange))
-        speclabel = "Frequency (PHz)"
-    elseif specaxis == :ω
-        specxfac = 1e-15
-        speclims = (specxfac*wlfreq(maximum(λrange)), specxfac*wlfreq(minimum(λrange)))
-        speclabel = "Angular frequency (rad/fs)"
-    elseif specaxis == :λ
-        specxfac = 1e9
-        speclims = λrange .* specxfac
-        speclabel = "Wavelength (nm)"
-    else
-        error("Unknown specaxis $specaxis")
-    end
-    return speclims, speclabel, specxfac
-end
-
-# single-mode 2D propagation plots
-function _prop2D_sm(t, z, specx, It, Iω, speclabel, speclims, trange, dBmin, bpstr; kwargs...)
-    id = "($(string(hash(gensym()); base=16)[1:4])) "
-    num = id * "Propagation" * ((length(bpstr) > 0) ? ", $bpstr" : "")
-    pfig, axs = plt.subplots(1, 2, num=num)
-    pfig.set_size_inches(12, 4)
-    Iω = Maths.normbymax(Iω)
-    _spec2D_log(axs[1], specx, z, Iω, dBmin, speclabel, speclims; kwargs...)
-
-    _time2D(axs[2], t, z, It, trange; kwargs...)
-    pfig.tight_layout()
-    return pfig
-end
-
-# multi-mode 2D propagation plots
-function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
-                    speclabel, speclims, trange, dBmin, bpstr;
-                    kwargs...)
-    pfigs = Figure[]
-    Iω = Maths.normbymax(Iω)
-    id = "($(string(hash(gensym()); base=16)[1:4])) "
-    for mi in modes
-        num = id * "Propagation ($(modelabels[mi]))" * ((length(bpstr) > 0) ? ", $bpstr" : "")
-        pfig, axs = plt.subplots(1, 2, num=num)
-        pfig.set_size_inches(12, 4)
-        _spec2D_log(axs[1], specx, z, Iω[:, mi, :], dBmin, speclabel, speclims; kwargs...)
-
-        _time2D(axs[2], t, z, It[:, mi, :], trange; kwargs...)
-        push!(pfigs, pfig)
-    end
-
-    num = id * "Propagation (all modes)" * ((length(bpstr) > 0) ? ", $bpstr" : "")
-    pfig, axs = plt.subplots(1, 2, num=num)
-    pfig.set_size_inches(12, 4)
-    Iωall = dropdims(sum(Iω, dims=2), dims=2)
-    _spec2D_log(axs[1], specx, z, Iωall, dBmin, speclabel, speclims; kwargs...)
-
-    Itall = dropdims(sum(It, dims=2), dims=2)
-    _time2D(axs[2], t, z, Itall, trange; kwargs...)
-    pfig.tight_layout()
-    push!(pfigs, pfig)
-
-    return pfigs
-end
-
-# a single logarithmic colour-scale spectral domain plot
-function _spec2D_log(ax, specx, z, I, dBmin, speclabel, speclims; kwargs...)
-    im = ax.pcolormesh(specx, z, 10*log10.(transpose(I)); shading="auto", kwargs...)
-    im.set_clim(dBmin, 0)
-    cb = plt.colorbar(im, ax=ax)
-    cb.set_label("SED (dB)")
-    ax.set_ylabel("Distance (cm)")
-    ax.set_xlabel(speclabel)
-    ax.set_xlim(speclims...)
-end
-
-# a single time-domain propagation plot
-function _time2D(ax, t, z, I, trange; kwargs...)
-    Pfac, unit = power_unit(I)
-    im = ax.pcolormesh(t*1e15, z, Pfac*transpose(I); shading="auto", kwargs...)
-    cb = plt.colorbar(im, ax=ax)
-    cb.set_label("Power ($unit)")
-    ax.set_xlim(trange.*1e15)
-    ax.set_xlabel("Time (fs)")
-    ax.set_ylabel("Distance (cm)")
+    getext().prop_2D(output, specaxis; trange, bandpass,
+                     λrange, dBmin, resolution, modes, oversampling, kwargs...)
 end
 
 """
@@ -367,77 +71,17 @@ a single index, a `range` or `:sum`. In the latter case, the sum of modes is plo
 
 The keyword argument `oversampling` determines the amount of oversampling done before plotting.
 
-Other `kwargs` are passed onto `plt.plot`.
+Other `kwargs` are passed onto `pyplot.plot`.
 """
 function time_1D(output, zslice=maximum(output["z"]);
                 y=:Pt, modes=nothing,
                 oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing,
                 FTL=false, propagate=nothing,
                 kwargs...)
-    t, Et, zactual = getEt(output, zslice,
-                           trange=trange, oversampling=oversampling, bandpass=bandpass,
-                           FTL=FTL, propagate=propagate)
-    if y == :Pt
-        yt = abs2.(Et)
-    elseif y == :Et
-        yt = real(Et)
-    elseif y == :Esq
-        yt = real(Et).^2
-    else
-        error("unknown time plot variable $y")
-    end
-    multimode, modestrs = get_modes(output)
-    if multimode
-        if modes == :sum
-            y == :Pt || error("Modal sum can only be plotted for power!")
-            yt = dropdims(sum(yt, dims=2), dims=2)
-            modestrs = join(modestrs, "+")
-            nmodes = 1
-        else
-            isnothing(modes) && (modes = 1:length(modestrs))
-            yt = yt[:, modes, :]
-            modestrs = modestrs[modes]
-            nmodes = length(modes)
-        end
-    end
-
-    yfac, unit = power_unit(abs2.(Et), y)
-
-    sfig = plt.figure()
-    if multimode && nmodes > 1
-        _plot_slice_mm(plt.gca(), t*1e15, yfac*yt, zactual, modestrs; kwargs...)
-    else
-        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
-        label = multimode ? zs.*" ($modestrs)" : zs
-        for iz in eachindex(zactual)
-            plt.plot(t*1e15, yfac*yt[:, iz]; label=label[iz], kwargs...)
-        end
-    end
-    plt.legend(frameon=false)
-    add_fwhm_legends(plt.gca(), "fs")
-    plt.xlabel("Time (fs)")
-    plt.xlim(1e15.*trange)
-    ylab = y == :Et ?  "Field ($unit)" : "Power ($unit)"
-    plt.ylabel(ylab)
-    y == :Et || plt.ylim(ymin=0)
-    sfig.set_size_inches(8.5, 5)
-    sfig.tight_layout()
-    sfig
+    getext().time_1D(output, zslice; y, modes, oversampling, trange, bandpass,
+                     FTL, propagate, kwargs...)
 end
 
-# Automatically find power unit depending on scale of electric field.
-function power_unit(Pt, y=:Pt)
-    units = ["kW", "MW", "GW", "TW", "PW"]
-    Pmax = maximum(Pt)
-    oom = clamp(floor(Int, log10(Pmax)/3), 1, 5) # maximum unit is PW
-    powerfac = 1/10^(oom*3)
-    if y == :Et
-        sqrt(powerfac), "$(units[oom])\$^{1/2}\$"
-    else
-        return powerfac, units[oom]
-    end
-end    
-
 """
     spec_1D(output, zslice, specaxis=:λ, log10=true, log10min=1e-6)
 
@@ -450,72 +94,16 @@ If `log10` is true, plot on a logarithmic scale, with a y-axis range of `log10mi
 The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
 a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.
 
-Other `kwargs` are passed onto `plt.plot`.
+Other `kwargs` are passed onto `pyplot.plot`.
 """
 function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
                  modes=nothing, λrange=(150e-9, 1200e-9),
                  log10=true, log10min=1e-6, resolution=nothing,
                  kwargs...)
-    if specaxis == :λ
-        specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution)
-    else
-        specx, Iω, zactual = getIω(output, specaxis, zslice, resolution=resolution)
-    end
-    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
-    multimode, modestrs = get_modes(output)
-    if multimode
-        modes = isnothing(modes) ? (1:size(Iω, 2)) : modes
-        if modes == :sum
-            Iω = dropdims(sum(Iω, dims=2), dims=2)
-            modestrs = join(modestrs, "+")
-            nmodes = 1
-        else
-            isnothing(modes) && (modes = 1:length(modestrs))
-            Iω = Iω[:, modes, :]
-            modestrs = modestrs[modes]
-            nmodes = length(modes)
-        end
-    end
-
-    specx .*= specxfac
-
-    sfig = plt.figure()
-    if multimode && nmodes > 1
-        _plot_slice_mm(plt.gca(), specx, Iω, zactual, modestrs, log10; kwargs...)
-    else
-        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
-        label = multimode ? zs.*" ($modestrs)" : zs
-        for iz in eachindex(zactual)
-            (log10 ? plt.semilogy : plt.plot)(specx, Iω[:, iz]; label=label[iz], kwargs...)
-        end
-    end
-    plt.legend(frameon=false)
-    plt.xlabel(speclabel)
-    plt.ylabel("Spectral energy density")
-    log10 && plt.ylim(3*maximum(Iω)*log10min, 3*maximum(Iω))
-    plt.xlim(speclims...)
-    sfig.set_size_inches(8.5, 5)
-    sfig.tight_layout()
-    sfig
+    getext().spec_1D(output, zslice, specaxis;
+                     modes, λrange, log10, log10min, resolution, kwargs...)
 end
 
-dashes = [(0, (10, 1)),
-          (0, (5, 1)),
-          (0, (1, 0.5)),
-          (0, (1, 0.5, 1, 0.5, 3, 1)),
-          (0, (5, 1, 1, 1))]
-
-function _plot_slice_mm(ax, x, y, z, modestrs, log10=false, fwhm=false; kwargs...)
-    pfun = (log10 ? ax.semilogy : ax.plot)
-    for sidx = 1:size(y, 3) # iterate over z-slices
-        zs = @sprintf("%.2f cm", z[sidx]*100)
-        line = pfun(x, y[:, 1, sidx]; label="$zs ($(modestrs[1]))", kwargs...)[1]
-        for midx = 2:size(y, 2) # iterate over modes
-            pfun(x, y[:, midx, sidx], linestyle=dashes[midx], color=line.get_color(),
-                 label="$zs ($(modestrs[midx]))"; kwargs...)
-        end
-    end
-end
 
 spectrogram(output::AbstractOutput, args...; kwargs...) = spectrogram(
     makegrid(output), output, args...; kwargs...)
@@ -534,146 +122,14 @@ function spectrogram(grid::Grid.AbstractGrid, output, zslice, specaxis=:λ;
 end
 
 function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
-                     trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40,
-                     kwargs...)
-    ω = Maths.rfftfreq(t)[2:end]
-    tmin, tmax = extrema(trange)
-    tg = collect(range(tmin, tmax, length=N))
-    g = Maths.gabor(t, real(Et), tg, fw)
-    g = g[2:end, :]
-
-    specy, Ig = getIω(ω, g*Maths.rfftnorm(t[2]-t[1]), specaxis)
-    speclims, speclabel, specyfac = getspeclims(λrange, specaxis)
-
-    log && (Ig = 10*log10.(Maths.normbymax(Ig)))
-
-    fig = plt.figure()
-    plt.pcolormesh(tg.*1e15, specyfac*specy, Ig; shading="auto", kwargs...)
-    plt.ylim(speclims...)
-    plt.ylabel(speclabel)
-    plt.xlabel("Time (fs)")
-    log && plt.clim(dBmin, 0)
-    plt.colorbar()
-    fig
+    trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40,
+    kwargs...)
+    getext().spectrogram(t, Et, specaxis;
+                         trange, N, fw, λrange, log, dBmin, kwargs...)
 end
 
 function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))
-    e = Processing.energy(output; bandpass=bandpass)
-    eall = Processing.energy(output)
-
-    multimode, modestrs = get_modes(output)
-    if multimode
-        e0 = sum(eall[:, 1])
-        modes = isnothing(modes) ? (1:size(e, 1)) : modes
-        if modes == :sum
-            e = dropdims(sum(e, dims=1), dims=1)
-            modestrs = join(modestrs, "+")
-            nmodes = 1
-        else
-            isnothing(modes) && (modes = 1:length(modestrs))
-            e = e[modes, :]
-            modestrs = modestrs[modes]
-            nmodes = length(modes)
-        end
-    else
-        e0 = eall[1]
-    end
-
-    z = output["z"]*100
-
-    fig = plt.figure()
-    ax = plt.axes()
-    ax.plot(z, 1e6*e')
-    ax.set_xlim(extrema(z)...)
-    ax.set_ylim(ymin=0)
-    ax.set_xlabel("Distance (cm)")
-    ax.set_ylabel("Energy (μJ)")
-    rax = ax.twinx()
-    rax.plot(z, 100*(e/e0)', linewidth=0)
-    lims = ax.get_ylim()
-    rax.set_ylim(100/(1e6*e0).*lims)
-    rax.set_ylabel("Conversion efficiency (%)")
-    fig.set_size_inches(figsize...)
-    fig
-end
-
-
-function auto_fwhm_arrows(ax, x, y; color="k", arrowlength=nothing, hpad=0, linewidth=1,
-                                    text=nothing, units="fs", kwargs...)
-    left, right = Maths.level_xings(x, y; kwargs...)
-    fw = abs(right - left)
-    halfmax = maximum(y)/2
-    arrowlength = isnothing(arrowlength) ? 2*fw : arrowlength
-
-    ax.annotate("", xy=(left-hpad, halfmax),
-                xytext=(left-hpad-arrowlength, halfmax),
-                arrowprops=Dict("arrowstyle" => "->",
-                                "color" => color,
-                                "linewidth" => linewidth))
-    ax.annotate("", xy=(right+hpad, halfmax),
-                xytext=(right+hpad+arrowlength, halfmax),
-                arrowprops=Dict("arrowstyle" => "->",
-                                "color" => color,
-                                "linewidth" => linewidth))
-
-    if text == :left
-        ax.text(left-arrowlength/2, 1.1*halfmax, @sprintf("%.2f %s", fw, units),
-                ha="right", color=color)
-    elseif text == :right
-        ax.text(right+arrowlength/2, 1.1*halfmax, @sprintf("%.2f %s", fw, units),
-                color=color)
-    end
-end
-
-function add_fwhm_legends(ax, unit)
-    leg = ax.get_legend()
-    texts = leg.get_texts()
-    handles, labels = ax.get_legend_handles_labels()
-    
-    for (ii, line) in enumerate(handles)
-        xy = line.get_xydata()
-        fw = Maths.fwhm(xy[:, 1], xy[:, 2])
-        t = texts[ii]
-        s = t.get_text()
-        s *= @sprintf(" [%.2f %s]", fw, unit)
-        t.set_text(s)
-    end
-end
-
-"""
-    cornertext(ax, text;
-               corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
-
-Place a `text` in the axes `ax` in the corner defined by `corner`. Padding can be
-defined for `x` and `y` together via `pad` or separately via `xpad` and `ypad`. Further
-keyword arguments are passed to `plt.text`. 
-
-Possible values for `corner` are `ul`, `ur`, `ll`, `lr` where the first letter
-defines upper/lower and the second defines left/right.
-"""
-function cornertext(ax, text; corner="ul", pad=0.02, xpad=nothing, ypad=nothing, kwargs...)
-    xpad = isnothing(xpad) ? pad : xpad
-    ypad = isnothing(ypad) ? pad : ypad
-    if corner[1] == 'u'
-        val = "top"
-        y = 1 - ypad
-    elseif corner[1] == 'l'
-        val = "bottom"
-        y = ypad
-    else
-        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
-    end
-    if corner[2] == 'l'
-        hal = "left"
-        x = xpad
-    elseif corner[2] == 'r'
-        hal = "right"
-        x = 1 - xpad
-    else
-        error("Invalid corner $corner. Must be one of ul, ur, ll, lr")
-    end
-    ax.text(x, y, text; horizontalalignment=hal, verticalalignment=val,
-                 transform=ax.transAxes, kwargs...)
+    getext().energy(output; modes, bandpass, figsize)
 end
 
 end
\ No newline at end of file
diff --git a/src/SFA.jl b/src/SFA.jl
index 71aaadbe..935ffde5 100644
--- a/src/SFA.jl
+++ b/src/SFA.jl
@@ -4,8 +4,6 @@ import FFTW
 import DataStructures: CircularBuffer
 import Luna: Maths, PhysData, Ionisation
 import Luna.PhysData: wlfreq, c
-import PyPlot: plt
-
 
 """
     approx_dipole(gas)
diff --git a/test/test_ionisation.jl b/test/test_ionisation.jl
index 74dbc0ad..5ac40a94 100644
--- a/test/test_ionisation.jl
+++ b/test/test_ionisation.jl
@@ -75,3 +75,4 @@ ratefun!(outneg, -E)
     adk_avg_kw = Ionisation.ionrate_ADK.(gas, E0; cycle_average=true)
     @test all(isapprox.(adk_avg, adk_avg_kw; rtol=0.05))
 end
+