Make Plotting optional and extendible

.gitignore

@@ -10,3 +10,5 @@ _git2_*
 docs/build
 *.info
 .vscode/spellright.dict
+.CondaPkg
+deps/build.log

Project.toml

@@ -3,6 +3,11 @@ uuid = "30eb0fb0-5147-11e9-3356-d75b018717ce"
 authors = ["chrisbrahms <[email protected]>"]
 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"

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,26 +138,32 @@ 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.
 
 ## 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/).
 

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)

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

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

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

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

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

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

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)

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)

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)

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

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