Colorbars and titles (#307)

* Include color bar when plotting intensity heatmap.
* Add title keyword argument to easily label intensity plots.
* Hide axisopts while we consider UI alternatives.

examples/01_LSWT_CoRh2O4.jl

@@ -155,7 +155,7 @@ path = q_space_path(cryst, qs, 500)
 
 energies = range(0, 6, 300)
 res = intensities(swt, path; energies, kernel)
-plot_intensities(res; units)
+plot_intensities(res; units, title="CoRh₂O₄ LSWT")
 
 # Sometimes experimental data is only available as a powder average, i.e., as an
 # average over all possible crystal orientations. Use [`powder_average`](@ref)
@@ -169,7 +169,7 @@ radii = range(0, 3, 200) # (1/Å)
 res = powder_average(cryst, radii, 2000) do qs
     intensities(swt, qs; energies, kernel)
 end
-plot_intensities(res; units, saturation=1.0)
+plot_intensities(res; units, saturation=1.0, title="CoRh₂O₄ Powder Average")
 
 # This result can be compared to experimental neutron scattering data
 # from Fig. 5 of [Ge et al.](https://doi.org/10.1103/PhysRevB.96.064413)

examples/02_LLD_CoRh2O4.jl

@@ -119,7 +119,7 @@ grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 200), [0, 1, 0], (-10, 10))
 # are above the ordering temperature, and do not have sharp Bragg peaks.
 
 res = intensities_static(sc, grid)
-plot_intensities(res; saturation=1.0)
+plot_intensities(res; saturation=1.0, title="Static Intensities at T = 16 K")
 
 # ### Dynamical structure factor
 
@@ -164,7 +164,7 @@ path = q_space_path(cryst, qs, 500)
 # Calculate and plot the intensities along this path.
 
 res = intensities(sc, path; energies, langevin.kT)
-plot_intensities(res; units)
+plot_intensities(res; units, title="Intensities at 16 K")
 
 # ### Powder averaged intensity
 
@@ -176,4 +176,4 @@ radii = range(0, 3.5, 200) # (1/Å)
 res = powder_average(cryst, radii, 350) do qs
     intensities(sc, qs; energies, langevin.kT)
 end
-plot_intensities(res; units)
+plot_intensities(res; units, title="Powder Average at 16 K")

examples/03_LSWT_SU3_FeI2.jl

@@ -214,7 +214,7 @@ swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min))
 qs = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]]
 path = q_space_path(cryst, qs, 500)
 res = intensities_bands(swt, path)
-plot_intensities(res; units)
+plot_intensities(res; units, title="Single Crystal Bands")
 
 # To make direct comparison with inelastic neutron scattering (INS) data, we
 # must account for empirical broadening of the data. Model this using a
@@ -233,7 +233,7 @@ weights = [1, 1, 1]
 res = domain_average(cryst, path; rotations, weights) do path_rotated
     intensities(swt, path_rotated; energies, kernel)
 end
-plot_intensities(res; units, colormap=:viridis)
+plot_intensities(res; units, colormap=:viridis, title="Domain Averaged Intensities")
 
 # This result can be directly compared to experimental neutron scattering data
 # from [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1)

examples/04_GSD_FeI2.jl

@@ -83,8 +83,8 @@ minimize_energy!(sys; maxiters=10)
 suggest_timestep(sys, langevin; tol=1e-2)
 langevin.dt = 0.03;
 
-# Run a Langevin trajectory for 10,000 time-steps and plot the spins. At this
-# angle, it is difficult to discern the magnetic order.
+# Run a Langevin trajectory for 10,000 time-steps and plot the spins. The
+# magnetic order is present, but may be difficult to see.
 
 for _ in 1:10_000
     step!(sys, langevin)
@@ -174,12 +174,12 @@ qs = [[0,   0, 0],  # List of wave vectors that define a path
       [0,   0, 0]] 
 qpath = q_space_path(cryst, qs, 500)
 res = intensities(sc, qpath; energies, langevin.kT)
-plot_intensities(res; colorrange=(0.0, 1.0))
+plot_intensities(res; colorrange=(0.0, 1.0), title="Intensities at T = 2.3 K")
 
 # One can also view the intensity along a [`q_space_grid`](@ref) for a fixed
 # energy value. Alternatively, use [`intensities_static`](@ref) to integrate
 # over all available energies.
 
 grid = q_space_grid(cryst, [1, 0, 0], range(-1.5, 1.5, 300), [0, 1, 0], (-1.5, 1.5); orthogonalize=true)
-res = intensities(sc, grid; energies=[3.88], langevin.kT)
-plot_intensities(res)
+res = intensities(sc, grid; energies=[3.5], langevin.kT)
+plot_intensities(res; title="Intensity slice at ω = 3.5 meV")

examples/07_Dipole_Dipole.jl

@@ -63,8 +63,8 @@ res3 = intensities_bands(swt, path)
 # slight corrections are needed for the third dispersion band.
 
 fig = Figure(size=(768, 300))
-plot_intensities!(fig[1, 1], res1; units)
-ax = plot_intensities!(fig[1, 2], res2; units)
+plot_intensities!(fig[1, 1], res1; units, title="Local Exchange Only")
+ax = plot_intensities!(fig[1, 2], res2; units, title="Local Exchange and Dipole-Dipole")
 for c in eachrow(res3.disp)
     lines!(ax, eachindex(c), c; linestyle=:dash, color=:black)
 end

examples/08_Momentum_Conventions.jl

@@ -86,6 +86,6 @@ res2 = intensities_bands(swt, path)
 # elastic peak at ``𝐪 = [0,0,0]``, reflecting the ferromagnetic ground state.
 
 fig = Figure(size=(768, 300))
-plot_intensities!(fig[1, 1], res1; axisopts=(; ylabel="", title="Classical dynamics"))
-plot_intensities!(fig[1, 2], res2; axisopts=(; ylabel="", title="Spin wave theory"))
+plot_intensities!(fig[1, 1], res1; title="Classical dynamics")
+plot_intensities!(fig[1, 2], res2; title="Spin wave theory")
 fig

examples/spinw_tutorials/SW08_sqrt3_kagome_AFM.jl

@@ -50,10 +50,10 @@ path = q_space_path(cryst, qs, 400)
 fig = Figure(size=(768, 300))
 swt = SpinWaveTheory(sys_enlarged; measure=ssf_perp(sys_enlarged))
 res = intensities_bands(swt, path)
-plot_intensities!(fig[1, 1], res; units, saturation=0.5, axisopts=(; title="Supercell method"))
+plot_intensities!(fig[1, 1], res; units, saturation=0.5,title="Supercell method")
 swt = SpinWaveTheorySpiral(sys; measure=ssf_perp(sys), k, axis)
 res = intensities_bands(swt, path)
-plot_intensities!(fig[1, 2], res; units, saturation=0.5, axisopts=(; title="Spiral method"))
+plot_intensities!(fig[1, 2], res; units, saturation=0.5, title="Spiral method")
 fig
 
 # Calculate and plot the powder averaged spectrum. Continuing to use the "spiral

examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl

@@ -82,8 +82,7 @@ qs = [[0, 1, -1], [0, 1, -1+1], [0, 1, -1+2], [0, 1, -1+3]]
 path = q_space_path(cryst, qs, 400)
 energies = range(0, 6, 400)
 res = intensities(swt, path; energies, kernel=gaussian(fwhm=0.25))
-axisopts = (; title=L"$ϵ_T = %$ϵT$", titlesize=20)
-plot_intensities(res; units, axisopts, saturation=0.7, colormap=:jet)
+plot_intensities(res; units, saturation=0.7, colormap=:jet, title="Scattering intensities")
 
 # Use [`ssf_custom_bm`](@ref) to calculate the imaginary part of
 # ``\mathcal{S}^{2, 3}(𝐪, ω) - \mathcal{S}^{3, 2}(𝐪, ω)``. In polarized
@@ -96,5 +95,4 @@ measure = ssf_custom_bm(sys; u=[0, 1, 0], v=[0, 0, 1]) do q, ssf
 end
 swt = SpinWaveTheorySpiral(sys; measure, k, axis)
 res = intensities(swt, path; energies, kernel=gaussian(fwhm=0.25))
-axisopts = (; title=L"$ϵ_T = %$ϵT$", titlesize=20)
-plot_intensities(res; units, axisopts, saturation=0.8, allpositive=false)
+plot_intensities(res; units, saturation=0.8, allpositive=false, title="Im[S²³(q, ω) - S³²(q, ω)]")

examples/spinw_tutorials/SW19_Different_Ions.jl

@@ -52,17 +52,17 @@ fig = Figure(size=(768,600))
 formfactors = [1 => FormFactor("Cu2"), 2 => FormFactor("Fe2")]
 swt = SpinWaveTheory(sys; measure=ssf_perp(sys; formfactors))
 res = intensities_bands(swt, path)
-plot_intensities!(fig[1, 1], res; units, axisopts=(; title="All correlations"))
+plot_intensities!(fig[1, 1], res; units, title="All correlations")
 
 formfactors = [1 => FormFactor("Cu2"), 2 => zero(FormFactor)]
 swt = SpinWaveTheory(sys; measure=ssf_perp(sys; formfactors))
 res = intensities_bands(swt, path)
-plot_intensities!(fig[1, 2], res; units, axisopts=(; title="Cu-Cu correlations"))
+plot_intensities!(fig[1, 2], res; units, title="Cu-Cu correlations")
 
 formfactors = [1 => zero(FormFactor), 2 => FormFactor("Fe2")]
 swt = SpinWaveTheory(sys; measure=ssf_perp(sys; formfactors))
 res = intensities_bands(swt, path)
-plot_intensities!(fig[2, 2], res; units, axisopts=(; title="Fe-Fe correlations"))
+plot_intensities!(fig[2, 2], res; units, title="Fe-Fe correlations")
 
 fig
 

ext/PlottingExt.jl

@@ -1001,7 +1001,7 @@ end
 
 function get_unit_energy(units, into)
     if isnothing(units)
-        isnothing(into) || error("units parameter required")
+        isnothing(into) || error("`into` parameter requires `units` parameter")
         return (1.0, "Energy")
     else
         sym = @something into units.energy
@@ -1032,13 +1032,15 @@ panel of a Makie figure.
 
 ```julia
 fig = Figure()
-plot_intensities!(fig[1, 1], res1)
-plot_intensities!(fig[2, 1], res2)
+plot_intensities!(fig[1, 1], res1; title="Panel 1")
+plot_intensities!(fig[2, 1], res2; title="Panel 2")
 display(fig)
 ```
 """
-function Sunny.plot_intensities!(panel, res::Sunny.BandIntensities{Float64}; colormap=nothing, saturation=0.9, sensitivity=0.0025, allpositive=true, units=nothing, into=nothing, fwhm=nothing, ylims=nothing, axisopts=Dict())
+function Sunny.plot_intensities!(panel, res::Sunny.BandIntensities{Float64}; colormap=nothing, saturation=0.9, sensitivity=0.0025, allpositive=true,
+                                 units=nothing, into=nothing, fwhm=nothing, ylims=nothing, title="", axisopts=NamedTuple())
     unit_energy, ylabel = get_unit_energy(units, into)
+    axisopts = (; title, axisopts...)
 
     if res.qpts isa Sunny.QPath 
         mindisp, maxdisp = extrema(res.disp)
@@ -1099,7 +1101,8 @@ function suggest_labels_for_grid(grid::Sunny.QGrid{N}) where N
 end
 
 
-function Sunny.plot_intensities!(panel, res::Sunny.Intensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, axisopts=Dict())
+function Sunny.plot_intensities!(panel, res::Sunny.Intensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, title="", axisopts=NamedTuple())
+    axisopts = (; title, axisopts...)
     (; crystal, qpts, data, energies) = res
 
     colorrange_suggest = colorrange_from_data(; data, saturation, sensitivity=0, allpositive)
@@ -1109,16 +1112,18 @@ function Sunny.plot_intensities!(panel, res::Sunny.Intensities{Float64}; colorma
     if qpts isa Sunny.QPath
         unit_energy, ylabel = get_unit_energy(units, into)
         xticklabelrotation = maximum(length.(qpts.xticks[2])) > 3 ? π/6 : 0.0
-        ax = Makie.Axis(panel; xlabel="Momentum (r.l.u.)", ylabel, qpts.xticks, xticklabelrotation, axisopts...)
-        Makie.heatmap!(ax, axes(data, 2), collect(energies/unit_energy), data'; colormap, colorrange)
+        ax = Makie.Axis(panel[1, 1]; xlabel="Momentum (r.l.u.)", ylabel, qpts.xticks, xticklabelrotation, axisopts...)
+        hm = Makie.heatmap!(ax, axes(data, 2), collect(energies/unit_energy), data'; colormap, colorrange)
+        Makie.Colorbar(panel[1, 2], hm)
         return ax
     elseif qpts isa Sunny.QGrid{2}
         if isone(length(energies))
             aspect = grid_aspect_ratio(crystal, qpts)
             xlabel, ylabel = suggest_labels_for_grid(qpts)
             (xs, ys) = map(range, qpts.coefs_lo, qpts.coefs_hi, size(qpts.qs))
-            ax = Makie.Axis(panel; xlabel, ylabel, aspect, axisopts...)
-            Makie.heatmap!(ax, xs, ys, dropdims(data; dims=1); colormap, colorrange)
+            ax = Makie.Axis(panel[1, 1]; xlabel, ylabel, aspect, axisopts...)
+            hm = Makie.heatmap!(ax, xs, ys, dropdims(data; dims=1); colormap, colorrange)
+            Makie.Colorbar(panel[1, 2], hm)
             return ax
         else
             error("Cannot yet plot $(typeof(res))")
@@ -1128,7 +1133,8 @@ function Sunny.plot_intensities!(panel, res::Sunny.Intensities{Float64}; colorma
     end
 end
 
-function Sunny.plot_intensities!(panel, res::Sunny.StaticIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, axisopts=Dict())
+function Sunny.plot_intensities!(panel, res::Sunny.StaticIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, title="", axisopts=NamedTuple())
+    axisopts = (; title, axisopts...)
     (; crystal, qpts, data) = res
 
     colorrange_suggest = colorrange_from_data(; data=reshape(data, 1, size(data)...), saturation, sensitivity=0, allpositive)
@@ -1139,44 +1145,53 @@ function Sunny.plot_intensities!(panel, res::Sunny.StaticIntensities{Float64}; c
         aspect = grid_aspect_ratio(crystal, qpts)
         xlabel, ylabel = suggest_labels_for_grid(qpts)
         (xs, ys) = map(range, qpts.coefs_lo, qpts.coefs_hi, size(qpts.qs))
-        ax = Makie.Axis(panel; xlabel, ylabel, aspect, axisopts...)
-        Makie.heatmap!(ax, xs, ys, data; colormap, colorrange)
+        ax = Makie.Axis(panel[1, 1]; xlabel, ylabel, aspect, axisopts...)
+        hm = Makie.heatmap!(ax, xs, ys, data; colormap, colorrange)
+        Makie.Colorbar(panel[1, 2], hm)
         return ax
     elseif qpts isa Sunny.QPath
         xticklabelrotation = maximum(length.(qpts.xticks[2])) > 3 ? π/6 : 0.0
         ax = Makie.Axis(panel; xlabel="Momentum (r.l.u.)", ylabel="Intensity", qpts.xticks, xticklabelrotation, axisopts...)
         Makie.lines!(ax, data)
-        println("colorrange ", colorrange)
         Makie.ylims!(ax, colorrange)
         return ax
     else
         error("Cannot yet plot $(typeof(res))")
     end
 end
 
-function Sunny.plot_intensities!(panel, res::Sunny.PowderIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, axisopts=Dict())
+function Sunny.plot_intensities!(panel, res::Sunny.PowderIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, title="", axisopts=NamedTuple())
+    axisopts = (; title, axisopts...)
     unit_energy, ylabel = get_unit_energy(units, into)
     xlabel = isnothing(units) ? "Momentum " : "Momentum ($(Sunny.unit_strs[units.length])⁻¹)" 
 
     colorrange_suggest = colorrange_from_data(; res.data, saturation, sensitivity=0, allpositive)
     colormap = @something colormap (allpositive ? :gnuplot2 : :bwr)
     colorrange = @something colorrange colorrange_suggest
 
-    ax = Makie.Axis(panel; xlabel, ylabel, axisopts...)
-    Makie.heatmap!(ax, res.radii, collect(res.energies/unit_energy), res.data'; colormap, colorrange)
+    ax = Makie.Axis(panel[1, 1]; xlabel, ylabel, axisopts...)
+    hm = Makie.heatmap!(ax, res.radii, collect(res.energies/unit_energy), res.data'; colormap, colorrange)
+    Makie.Colorbar(panel[1, 2], hm)
     return ax
 end
 
+#=
+  * `axisopts`: An additional collection of named arguments that will be passed
+    to the `Makie.Axis` constructor. This allows to override the axis `xlabel`,
+    `ylabel`, `xticks`, etc. See [Makie
+    documentation](https://docs.makie.org/release/reference/blocks/axis#attributes)
+    for details.
+=#
 """
     plot_intensities(res::BandIntensities; colormap=nothing, allpositive=true, saturation=0.9,
                      sensitivity=0.0025, fwhm=nothing, ylims=nothing, units=nothing, into=nothing, 
-                     axisopts=Dict())
+                     title="")
 
     plot_intensities(res::Intensities; colormap=nothing, colorrange=nothing, allpositive=true, 
-                     saturation=0.9, units=nothing, into=nothing, axisopts=Dict())
+                     saturation=0.9, units=nothing, into=nothing, title="")
 
     plot_intensities(res::PowderIntensities; colormap=nothing, colorrange=nothing, allpositive=true, 
-                     saturation=0.9, units=nothing, into=nothing, axisopts=Dict())
+                     saturation=0.9, units=nothing, into=nothing, title="")
 
 Keyword arguments:
 
@@ -1197,9 +1212,7 @@ Keyword arguments:
     conversions.
   * `into`: A symbol for conversion into a new base energy unit (e.g. `:meV`,
     `:K`, etc.)
-  * `axisopts`: An additional collection of named arguments that will be passed
-    to the `Makie.Axis` constructor. This allows to override the axis `title`,
-    `xlabel`, `ylabel`, `xticks`, etc. See Makie documentation for details.
+  * `title`: An optional title for the plot.
 """
 function Sunny.plot_intensities(res::Sunny.AbstractIntensities; opts...)
     fig = Makie.Figure()