From 2bb4faf07c12a3ea5533dda863e859c01fa4c53b Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 8 Feb 2023 08:28:22 +0000
Subject: [PATCH 01/18] fix gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 9b1a5785..2b79f68a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,4 @@ _git2_*
 docs/build
 *.info
 .vscode/spellright.dict
+.CondaPkg
\ No newline at end of file

From b14dd0e5cab875a7d118142212403bbe8420d9b2 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 8 Feb 2023 08:58:17 +0000
Subject: [PATCH 02/18] fix deps

---
 Project.toml | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index e7e34d3c..bc94eb33 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 H5Zblosc = "c8ec2601-a99c-407f-b158-e79c03c2f5f7"
@@ -37,8 +38,6 @@ Pidfile = "fa939f87-e72e-5be4-a000-7fc836dbe307"
 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"
@@ -76,8 +75,6 @@ PhysicalConstants = "0.2"
 Pidfile = "1.2"
 Polynomials = "2, 3"
 ProgressLogging = "0.1"
-PyCall = "1.92"
-PyPlot = "2.9"
 QuadGK = "2.4"
 Reexport = "1.2"
 Roots = "2"

From 7a90185ef6af47787e0f94c50b7ecaf246a5b18e Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 8 Feb 2023 08:58:34 +0000
Subject: [PATCH 03/18] stash minimal progress

---
 src/Plotting.jl | 19 +++++--------------
 1 file changed, 5 insertions(+), 14 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 1c349f99..8fd9d3f8 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -3,26 +3,17 @@ 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 GLMakie
 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 newfig()
+    f = GLMakie.Figure(resolution = (1600, 1200))
+    display(GLMakie.Screen(), f)
+    f
 end
 
-display(figs::AbstractArray{Figure, N}) where N = [display(fig) for fig in figs]
-
 """
     cmap_white(cmap, N=512, n=8)
 

From 44d47dcf950b3863fb54647bf5630b67b2fc7a17 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 8 Feb 2023 09:47:31 +0000
Subject: [PATCH 04/18] working prop2D

---
 src/Plotting.jl | 60 +++++++++++++++++--------------------------------
 1 file changed, 21 insertions(+), 39 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 8fd9d3f8..f9700b6c 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -282,68 +282,50 @@ end
 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
+    _prop2D_fig(num, specx, z, Iω, dBmin, speclabel, speclims, t, It, trange)
 end
 
 # multi-mode 2D propagation plots
 function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
                     speclabel, speclims, trange, dBmin, bpstr;
                     kwargs...)
-    pfigs = Figure[]
+    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 = 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...)
+        pfig = _prop2D_fig(num, specx, z, Iω[:, mi, :], dBmin, speclabel, speclims, t, It[:, mi, :], trange)
         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()
+    pfig = _prop2D_fig(num, specx, z, Iωall, dBmin, speclabel, speclims, t, Itall, trange)
     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...)
+function _prop2D_fig(name, specx, z, Iω, dBmin, speclabel, speclims, t, It, trange)
+    pfig = newfig()
+    ax, hm = GLMakie.heatmap(pfig[1,1], specx, z, 10*log10.(Iω),
+                             colorrange=(dBmin,0), interpolate=true,
+                             lowclip=:white,
+                             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,
+                             lowclip=:white,
+                             axis=(; xlabel="Time (fs)", ylabel="Distance (cm)"))
+    GLMakie.xlims!(ax, trange.*1e15)
+    cb = GLMakie.Colorbar(pfig[1, 4], hm, label="Power ($unit)")
+    pfig
 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...)

From cf6c15848e32a9763b423fe9e81e224dd4f92ba3 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 8 Feb 2023 22:13:49 +0000
Subject: [PATCH 05/18] working spec1d

---
 src/Plotting.jl | 39 +++++++++++++++++++--------------------
 1 file changed, 19 insertions(+), 20 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index f9700b6c..cfc537ff 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -323,6 +323,7 @@ function _prop2D_fig(name, specx, z, Iω, dBmin, speclabel, speclims, t, It, tra
                              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
 
@@ -452,42 +453,40 @@ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
 
     specx .*= specxfac
 
-    sfig = plt.figure()
     if multimode && nmodes > 1
-        _plot_slice_mm(plt.gca(), specx, Iω, zactual, modestrs, log10; kwargs...)
+        sfig = _plot_slice_mm(specx, Iω, zactual, modestrs, speclabel, 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)
-            (log10 ? plt.semilogy : plt.plot)(specx, Iω[:, iz]; label=label[iz], kwargs...)
+            GLMakie.lines!(specx, Iω[:, iz], label=label[iz])
         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()
+    GLMakie.axislegend(framevisible=false)
+    log10 && GLMakie.ylims!(3*maximum(Iω)*log10min, 3*maximum(Iω))
+    GLMakie.xlims!(speclims...)
+    GLMakie.DataInspector()
     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))]
+dashes = [:dash, :dot, :dashdot, :dashdotdot, [0.5, 1.0, 1.5, 2.5]]
 
-function _plot_slice_mm(ax, x, y, z, modestrs, log10=false, fwhm=false; kwargs...)
-    pfun = (log10 ? ax.semilogy : ax.plot)
+function _plot_slice_mm(x, y, z, modestrs, speclabel, log10=false, fwhm=false)
+    pfig = newfig()
+    scale = (log10 ? Base.log10 : :identity)
+    GLMakie.Axis(pfig[1, 1], yscale = scale, xlabel=speclabel, ylabel="Spectral energy density")
     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]
+        line = GLMakie.lines!(x, y[:, 1, sidx], label="$zs ($(modestrs[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...)
+            GLMakie.lines!(x, y[:, midx, sidx], label="$zs ($(modestrs[midx]))",
+                   color=line[:color], linestyle=dashes[midx])
         end
     end
+    pfig
 end
 
 spectrogram(output::AbstractOutput, args...; kwargs...) = spectrogram(

From 8e31581e6a9364ddc7465151c3078cbf6588e0c1 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 8 Feb 2023 22:38:15 +0000
Subject: [PATCH 06/18] working time_1d

---
 src/Plotting.jl | 62 +++++++++++++++++++++++--------------------------
 1 file changed, 29 insertions(+), 33 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index cfc537ff..a5621a94 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -377,26 +377,26 @@ function time_1D(output, zslice=maximum(output["z"]);
 
     yfac, unit = power_unit(abs2.(Et), y)
 
-    sfig = plt.figure()
+    xlabel = "Time (fs)"
+    ylabel = y == :Et ?  "Field ($unit)" : "Power ($unit)"
     if multimode && nmodes > 1
-        _plot_slice_mm(plt.gca(), t*1e15, yfac*yt, zactual, modestrs; kwargs...)
+       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)
-            plt.plot(t*1e15, yfac*yt[:, iz]; label=label[iz], kwargs...)
+            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
-    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
+    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.
@@ -454,7 +454,7 @@ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
     specx .*= specxfac
 
     if multimode && nmodes > 1
-        sfig = _plot_slice_mm(specx, Iω, zactual, modestrs, speclabel, log10)
+        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]
@@ -474,15 +474,25 @@ end
 
 dashes = [:dash, :dot, :dashdot, :dashdotdot, [0.5, 1.0, 1.5, 2.5]]
 
-function _plot_slice_mm(x, y, z, modestrs, speclabel, log10=false, fwhm=false)
+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=speclabel, ylabel="Spectral energy density")
+    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)
-        line = GLMakie.lines!(x, y[:, 1, sidx], label="$zs ($(modestrs[1]))")
+        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
-            GLMakie.lines!(x, y[:, midx, sidx], label="$zs ($(modestrs[midx]))",
+            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], linestyle=dashes[midx])
         end
     end
@@ -597,20 +607,6 @@ function auto_fwhm_arrows(ax, x, y; color="k", arrowlength=nothing, hpad=0, line
     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;

From 408a7f44711886d09f882e655ce614f63fb3a4e1 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 11:45:59 +0000
Subject: [PATCH 07/18] spectrogram and energy plots

---
 src/Plotting.jl | 116 +++++++++++-------------------------------------
 1 file changed, 27 insertions(+), 89 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index a5621a94..611b46a5 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -526,16 +526,17 @@ function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
 
     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()
+    Ig = Maths.normbymax(Ig)
+    log && (Ig = 10*log10.(Ig))
+    clims = (log ? (dBmin, 0) : extrema(Ig))
+
+    fig = newfig()
+    ax, hm = GLMakie.heatmap(fig[1,1], tg.*1e15, specyfac*specy, Ig,
+                             colorrange=clims, interpolate=true,
+                             axis=(; xlabel="Time (fs)", ylabel=speclabel))
+    GLMakie.ylims!(ax, speclims)
+    GLMakie.Colorbar(fig[1, 2], hm)
+    GLMakie.DataInspector()
     fig
 end
 
@@ -562,86 +563,23 @@ function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))
     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)
+    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
-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 [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...)
+    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
 
 end
\ No newline at end of file

From 0b0a953c964997be96f9289fa211c37df79238d2 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 12:46:57 +0000
Subject: [PATCH 08/18] fix spectrograms

---
 src/Plotting.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 611b46a5..03150a68 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -531,7 +531,7 @@ function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
     clims = (log ? (dBmin, 0) : extrema(Ig))
 
     fig = newfig()
-    ax, hm = GLMakie.heatmap(fig[1,1], tg.*1e15, specyfac*specy, Ig,
+    ax, hm = GLMakie.heatmap(fig[1,1], tg.*1e15, specyfac*specy, Ig',
                              colorrange=clims, interpolate=true,
                              axis=(; xlabel="Time (fs)", ylabel=speclabel))
     GLMakie.ylims!(ax, speclims)

From 31b000a9a7529d7733edf1393f46e700573a77b0 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 13:03:47 +0000
Subject: [PATCH 09/18] 3d spectrogram

---
 src/Plotting.jl | 30 ++++++++++++++++++++++++------
 1 file changed, 24 insertions(+), 6 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 03150a68..65975816 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -517,6 +517,7 @@ 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)
@@ -524,18 +525,35 @@ function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
     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)
+    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))
 
     fig = newfig()
-    ax, hm = GLMakie.heatmap(fig[1,1], tg.*1e15, specyfac*specy, Ig',
-                             colorrange=clims, interpolate=true,
-                             axis=(; xlabel="Time (fs)", ylabel=speclabel))
-    GLMakie.ylims!(ax, speclims)
-    GLMakie.Colorbar(fig[1, 2], hm)
+    if surface3d
+        ax, pl = GLMakie.surface(fig[1,1], tg.*1e15, specyfac*specy, Ig',
+                         colorrange=clims, colormap=:turbo,
+                         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,
+                                axis=(; xlabel="Time (fs)", ylabel=speclabel))
+        GLMakie.ylims!(ax, speclims)
+    end
+    GLMakie.Colorbar(fig[1, 2], pl)
     GLMakie.DataInspector()
     fig
 end

From c17b1b5d1e32cab6a57342cd6a507aaf4a0bac69 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 13:14:49 +0000
Subject: [PATCH 10/18] fix deps

---
 Project.toml | 3 +++
 src/SFA.jl   | 1 -
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index bc94eb33..59f55696 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.2.0"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CoolProp = "e084ae63-2819-5025-826e-f8e611a84251"
 Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
@@ -26,6 +27,7 @@ H5Zblosc = "c8ec2601-a99c-407f-b158-e79c03c2f5f7"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Hankel = "74863788-d124-456e-a676-9b76578dd39e"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 LibGit2 = "76f85450-5226-5b5a-8eaa-529ad045b433"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -39,6 +41,7 @@ Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
diff --git a/src/SFA.jl b/src/SFA.jl
index 71aaadbe..176311dd 100644
--- a/src/SFA.jl
+++ b/src/SFA.jl
@@ -4,7 +4,6 @@ import FFTW
 import DataStructures: CircularBuffer
 import Luna: Maths, PhysData, Ionisation
 import Luna.PhysData: wlfreq, c
-import PyPlot: plt
 
 
 """

From fce7319feec593bc7bd454c823b016d92fd37b2e Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 13:14:56 +0000
Subject: [PATCH 11/18] fix fonts

---
 src/Plotting.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 65975816..9aa080b2 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -8,6 +8,8 @@ import FFTW
 import Printf: @sprintf
 import Base: display
 
+GLMakie.set_theme!(GLMakie.Theme(fontsize = 20))
+
 function newfig()
     f = GLMakie.Figure(resolution = (1600, 1200))
     display(GLMakie.Screen(), f)

From 528668605e69393b1eb50a53f3393c471fd67b44 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 13:42:41 +0000
Subject: [PATCH 12/18] add stats

---
 src/Plotting.jl | 56 ++++++++++++++++++++++---------------------------
 1 file changed, 25 insertions(+), 31 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 9aa080b2..9e288d85 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -8,7 +8,7 @@ import FFTW
 import Printf: @sprintf
 import Base: display
 
-GLMakie.set_theme!(GLMakie.Theme(fontsize = 20))
+GLMakie.set_theme!(GLMakie.Theme(fontsize = 40))
 
 function newfig()
     f = GLMakie.Figure(resolution = (1600, 1200))
@@ -57,16 +57,7 @@ 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]
+    collect(Iterators.product(1:rows,1:cols))
 end
 
 """
@@ -157,37 +148,40 @@ function stats(output; kwargs...)
     z = stats["z"]*1e2
 
     multimode, modes = get_modes(output)
+    modes = isnothing(modes) ? [""] : modes
 
     Npl = length(pstats)
     if Npl > 0
-        pfig, axs = subplotgrid(Npl, title="Pulse stats")
+        pfig = newfig()
+        idcs = subplotgrid(Npl)
         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)
+            data, ylabel = pstats[n]
+            data = data'
+            scale = (multimode ? log10 : identity)
+            ax = GLMakie.Axis(pfig[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
-        pfig.tight_layout()
+        GLMakie.DataInspector()
     end
     
     Npl = length(fstats)
     if Npl > 0
-        ffig, axs = subplotgrid(Npl, title="Other stats")
+        ffig = newfig()
+        idcs = subplotgrid(Npl)
         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)
+            data, ylabel = fstats[n]
+            data = data'
+            scale = ((multimode && should_log10(data)) ? log10 : identity)
+            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
-        ffig.tight_layout()
+        GLMakie.DataInspector()
     end
     [pfig, ffig]
 end

From ccefe09e3e16a8385a1dc18d11ea42f0139f5eab Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Thu, 9 Feb 2023 14:06:46 +0000
Subject: [PATCH 13/18] fix superscripts in stats

---
 src/Plotting.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 9e288d85..938f105e 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -123,7 +123,7 @@ function stats(output; kwargs...)
         push!(pstats, (Pfac*stats["peakpower"], "Peak power ($unit)"))
     end
     haskey(stats, "peakintensity") && push!(
-        pstats, (1e-16*stats["peakintensity"], "Peak Intensity (TW/cm\$^2\$)"))
+        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)"))
@@ -131,14 +131,14 @@ function stats(output; kwargs...)
 
     fstats = [] # fibre/waveguide/propagation statistics
     if haskey(stats, "electrondensity")
-        push!(fstats, (1e-6*stats["electrondensity"], "Electron density (cm\$^{-3}\$)"))
+        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\$^{-3}\$)"))
+        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)"))

From 86590dd75299a52e2a8fe7634aa278fda6f782f1 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 6 Nov 2024 22:30:43 +0000
Subject: [PATCH 14/18] remove resolution

---
 src/Plotting.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 5660119f..03895250 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -11,7 +11,7 @@ import Base: display
 GLMakie.set_theme!(GLMakie.Theme(fontsize = 40))
 
 function newfig()
-    f = GLMakie.Figure(resolution = (1600, 1200))
+    f = GLMakie.Figure()
     display(GLMakie.Screen(), f)
     f
 end

From 43606e18fc9c8969876f3ccff03b245fc9023e77 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 6 Nov 2024 22:30:53 +0000
Subject: [PATCH 15/18] fix gitignore

---
 .gitignore | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 2b79f68a..0f7a834c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,4 +10,5 @@ _git2_*
 docs/build
 *.info
 .vscode/spellright.dict
-.CondaPkg
\ No newline at end of file
+.CondaPkg
+deps/build.log

From afbce89b2886c811896265868404c98d71e83337 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 6 Nov 2024 23:16:24 +0000
Subject: [PATCH 16/18] try to fix multimode stats plots

---
 src/Plotting.jl | 24 +++++++++++-------------
 1 file changed, 11 insertions(+), 13 deletions(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 03895250..761e39f5 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -10,8 +10,8 @@ import Base: display
 
 GLMakie.set_theme!(GLMakie.Theme(fontsize = 40))
 
-function newfig()
-    f = GLMakie.Figure()
+function newfig(; size=(800,600))
+    f = GLMakie.Figure(; size)
     display(GLMakie.Screen(), f)
     f
 end
@@ -54,10 +54,10 @@ Create a figure with `N` subplots laid out in a grid that is as close to square
 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)
+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))
+    collect(Iterators.product(1:rows,1:cols)), cols*colw, rows*rowh
 end
 
 """
@@ -152,12 +152,12 @@ function stats(output; kwargs...)
 
     Npl = length(pstats)
     if Npl > 0
-        pfig = newfig()
-        idcs = subplotgrid(Npl)
+        idcs, width, height = subplotgrid(Npl)
+        pfig = newfig(size=(width, height))
         for n in 1:Npl
             data, ylabel = pstats[n]
-            data = data'
             scale = (multimode ? log10 : identity)
+            data = multimode ? max.(data, 1e-300) : data
             ax = GLMakie.Axis(pfig[idcs[n]...]; xlabel="Distance (cm)", ylabel, yscale=scale)
             for i in 1:size(data,1)
                 GLMakie.lines!(z, data[i,:], label=modes[i])
@@ -169,12 +169,12 @@ function stats(output; kwargs...)
     
     Npl = length(fstats)
     if Npl > 0
-        ffig = newfig()
-        idcs = subplotgrid(Npl)
+        idcs, width, height = subplotgrid(Npl)
+        ffig = newfig(size=(width, height))
         for n in 1:Npl
             data, ylabel = fstats[n]
-            data = data'
             scale = ((multimode && should_log10(data)) ? log10 : identity)
+            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])
@@ -468,8 +468,6 @@ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
     sfig
 end
 
-dashes = [:dash, :dot, :dashdot, :dashdotdot, [0.5, 1.0, 1.5, 2.5]]
-
 function _plot_slice_mm(x, y, z, modestrs, xlabel, ylabel, log10=false; fwlabel=false)
     pfig = newfig()
     scale = (log10 ? Base.log10 : identity)
@@ -489,7 +487,7 @@ function _plot_slice_mm(x, y, z, modestrs, xlabel, ylabel, log10=false; fwlabel=
                 label *= @sprintf(" [%.2f %s]", fw, "fs")
             end
             GLMakie.lines!(x, y[:, midx, sidx]; label,
-                   color=line[:color], linestyle=dashes[midx])
+                   color=line[:color], cycle=[:linestyle])
         end
     end
     pfig

From f1ca301f8832f735aa8fbe2972ebb995cf429830 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 6 Nov 2024 23:33:23 +0000
Subject: [PATCH 17/18] fix stats

---
 src/Plotting.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 761e39f5..2dd614d0 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -157,8 +157,10 @@ function stats(output; kwargs...)
         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
@@ -174,6 +176,7 @@ function stats(output; kwargs...)
         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)

From fb64b5f7e6b932a576d6b9269348558afcf935e6 Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Wed, 6 Nov 2024 23:38:17 +0000
Subject: [PATCH 18/18] fix aspect ratio of prop2D

---
 src/Plotting.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Plotting.jl b/src/Plotting.jl
index 2dd614d0..5953158c 100644
--- a/src/Plotting.jl
+++ b/src/Plotting.jl
@@ -307,7 +307,7 @@ function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
 end
 
 function _prop2D_fig(name, specx, z, Iω, dBmin, speclabel, speclims, t, It, trange)
-    pfig = newfig()
+    pfig = newfig(size=(1000,400))
     ax, hm = GLMakie.heatmap(pfig[1,1], specx, z, 10*log10.(Iω),
                              colorrange=(dBmin,0), interpolate=true,
                              lowclip=:white,