Refactor plots

Project.toml

@@ -6,6 +6,7 @@ version = "0.2.0"
 [deps]
 JSServe = "824d6782-a2ef-11e9-3a09-e5662e0c26f9"
 LaserTypes = "e07c0bfa-524c-4f35-a151-c3dd916fa2f0"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 PICAnalysisTools = "049d77de-c438-40ac-b6d1-f1d179569bd2"
 PICDataStructures = "cb18667b-f199-4b8b-b886-9916fef47a26"
@@ -25,7 +26,7 @@ WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"
 [compat]
 JSServe = "1.2"
 LaserTypes = "0.1"
-PICDataStructures = "0.3"
+PICDataStructures = "0.4"
 PhysicalConstants = "0.2"
 Plots = "1.10"
 ProgressLogging = "0.1"

src/SDFResultViewer.jl

@@ -1,12 +1,14 @@
 module SDFResultViewer
 
-export add_λ_units, unit_l, unit_t, unit_E, unit_B, unit_L, pₑ,
+export add_λ_units, unit_l, unit_t, unit_E, unit_B, unit_L, unit_S, pₑ,
     batch_plot,
     mean_Lx_plot, Lx_section_plot,
     phase_space_summary_plot,
     average_linear_density_x, npart_plot,
-    section_widget,
-    compare_E_slice_with_analytic, similar_E, similar_laser
+    section_widget, approx_plot,
+    compare_E_slice_with_analytic, compare_S_slice_with_analytic,
+    poynting_summary,
+    apply_analytic, analytic_laser
 
 using SDFResults
 using PICDataStructures
@@ -18,8 +20,9 @@ using Rotations
 using RecursiveArrayTools: recursive_bottom_eltype
 using Unitful
 using Unitful: superscript
-using PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e
+using PhysicalConstants.CODATA2018: c_0, ε_0, μ_0, m_e, e
 using Statistics, StatsBase
+using LinearAlgebra
 using ProgressLogging
 # plots
 import Plots
@@ -30,6 +33,8 @@ using JSServe
 using JSServe: DOM, Slider
 using Markdown
 
+Unitful.register(@__MODULE__)
+
 include("units.jl")
 include("batch.jl")
 include("plots/angular_momentum.jl")

src/plots/fields.jl

@@ -1,33 +1,219 @@
-function E_slice_plot(E::ScalarField, dir, slice_location, t;
-    issliced=false,
-    save=false,
-    slice_dir=:x,
+function approx_plot(f;
+    fig = Figure(),
+    figurepos = fig[1,1],
+    xlabel = "",
+    ylabel = "",
+    title = "",
+    downsample_size = nothing,
+    cbar_label="",
     kwargs...)
 
-    if !issliced
-        E_slice = slice(E, slice_dir, slice_location)
+    ax = Axis(figurepos[1,1];
+        xlabel, ylabel,
+        title,
+        aspect = DataAspect()
+    )
+
+    if isnothing(downsample_size)
+        f_approx = downsample_approx(f)
     else
-        E_slice = E
+        f_approx = downsample_approx(f, downsample_size)
     end
+    @debug "Downsampled to $(size(f_approx))"
+
+    plt = fieldplot!(ax, f_approx; kwargs...)
+
+    Colorbar(figurepos[1,2], plt; width=20, label=cbar_label)
+
+    # TODO: Implement save
+    # if save
+    #     save && Plots.savefig(plt, "E$(dir_slice)_$loc.png")
+    # end
+
+    return fig
+end
+
+get_unit_label(f) = " (" * unitname(f) * ")"
+
+function compare_slice_with_analytic(numeric_f, f, laser, profile, file;
+    fig = Figure(resolution = (800, 700)),
+    component=1,
+    slice_location=0.0unit_l,
+    slice_dir=:x,
+    fieldname="",
+    f_unit)
 
-    grid = getdomain(E_slice)
-    cl = ustrip(max(abs.(extrema(E_slice))...))
+    fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l
+    t = get_time(file) |> unit_t
+    s = analytic_laser(laser, profile, file)
+
+    grid = getdomain(numeric_f)
+
+    analytic_f = apply_analytic(f, grid, t, s;
+        component,
+        propagation_dir=:x,
+        x₀=fx) |> f_unit
+
+    analytic_slice, loc = dynamic_slice(analytic_f, slice_dir, slice_location)
+    numeric_slice, loc = dynamic_slice(numeric_f, slice_dir, slice_location)
+
+    units = get_unit_label(grid)
+    xlabel, ylabel = filter(x->x≠string(slice_dir),
+        ["x", "y", "z"] .* units)
     ts = sprint(show, t, context=:compact => true)
-    loc = sprint(show, slice_location, context=:compact => true)
-    xlabel, ylabel = filter(x->x≠string(slice_dir), ["x", "y", "z"])
-
-    plt = Plots.heatmap(grid..., E_slice;
-        title = "E$dir at t=$ts, x=$loc",
-        xlabel, ylabel, colorbar_title = "E$dir",
-        seriescolor = :jet1,
-        aspect_ratio = 1,
-        framestyle = :box,
-        clims = (-cl, cl),
-        kwargs...)
-
-    if save
-        save && Plots.savefig(plt, "E$(dir_slice)_$loc.png")
-    end
+    loc_str = @lift sprint(show,
+        $loc, context=:compact => true) * " " * unitname(slice_location)
+    title = @lift fieldname * " slice at t = " * ts * " and " *
+        string(slice_dir) * " = " * $loc_str
+
+    ax1 = Axis(fig[1,1][1,1];
+        xlabel, ylabel,
+        title = "Analytic",
+        aspect = DataAspect()
+    )
+    ax2 = Axis(fig[1,2][1,1];
+        xlabel, ylabel,
+        title = "Numeric",
+        aspect = DataAspect()
+    )
+
+    plt1 = fieldplot!(ax1, analytic_slice)
+    plt2 = fieldplot!(ax2, numeric_slice)
+
+    Colorbar(fig[1,1][2,1], plt1;
+        width = Relative(3/4),
+        height = 20,
+        vertical = false, flipaxis = false,
+        tellheight = true,
+        label = fieldname * get_unit_label(numeric_f),
+    )
+    Colorbar(fig[1,2][2,1], plt2;
+        width = Relative(3/4),
+        height = 20,
+        vertical = false, flipaxis = false,
+        tellheight = true,
+        label = fieldname * get_unit_label(numeric_f),
+    )
+
+    Label(fig[0, :], title, textsize = 30, color = :black)
+
+    return fig
+end
+
+function compare_E_slice_with_analytic(laser, profile, file, dir;
+    slice_location=0.0unit_l,
+    slice_dir=:x)
+
+    E = uconvert(unit_E, unit_l, file["e$dir"])
+
+    compare_slice_with_analytic(E, LaserTypes.E, laser, profile, file;
+        component=dir_to_idx(dir),
+        slice_location,
+        slice_dir,
+        f_unit=unit_E,
+        fieldname="E$dir"
+    )
+end
+
+function compare_S_slice_with_analytic(laser, profile, file, dir;
+    slice_location=0.0unit_l,
+    slice_dir=:x)
+
+    Ex, Ey, Ez = uconvert.(unit_E, unit_l, file[:ex, :ey, :ez])
+    Bx, By, Bz = uconvert.(unit_B, unit_l, file[:bx, :by, :bz])
+
+    E = build_vector((Ex, Ey, Ez), (:x, :y, :z))
+    B = build_vector((Bx, By, Bz), (:x, :y, :z))
+    S = ustrip(E × B) / ustrip(μ_0)
+
+    compare_slice_with_analytic(S, LaserTypes.S, laser, profile, file;
+        component=dir_to_idx(dir),
+        slice_location=ustrip(slice_location),
+        slice_dir,
+        f_unit=unit_E,
+    )
+end
+
+function poynting_summary(file, laser, profile, slice_location)
+    t = get_time(file) |> unit_t
+
+    # get data
+    fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l
+    Ex, Ey, Ez, E, Bx, By, Bz, B = read_EM_fields(file)
+
+    Ex_grid = getdomain(Ex)
+    Bx_grid = getdomain(Bx)
+
+    s = analytic_laser(laser, profile, file)
+
+    analytic_E = get_analytic_E(Ex_grid, t, s, x₀=fx)
+    analytic_B = get_analytic_B(Bx_grid, t, s, x₀=fx)
+    analytic_S = uconvert(unit_S, analytic_E × analytic_B / μ_0)
+
+    S = uconvert(unit_S, E × B / μ_0)
+
+    S_approx = downsample_approx(S, 15)
+    analytic_S_approx = downsample_approx(analytic_S, 15)
+    S_slice = downsample_approx(
+        slice(S, :x, slice_location),
+        50
+    )
+    analytic_S_slice = downsample_approx(
+        slice(analytic_S, :x, slice_location),
+        50
+    )
+
+    fig = Figure(resolution=(800,800))
+    # Layout
+    lscene1 = LScene(fig[1, 1], scenekw = (camera = cam3d!, raw = false))
+    lscene2 = LScene(fig[1, 2], scenekw = (camera = cam3d!, raw = false))
+    ax1 = Axis(fig[2,1][1,1],
+        aspect = DataAspect(),
+        title="Numeric",
+        xlabel="y ($unit_l)",
+        ylabel="z ($unit_l)"
+    )
+    ax2 = Axis(fig[2,2][1,1],
+        aspect = DataAspect(),
+        title="Analytic",
+        xlabel="y ($unit_l)",
+        ylabel="z ($unit_l)"
+    )
+    # Plots
+    # 3D
+    arrowsize_factor = 0.02
+    lengthscale_factor = 10
+    fieldplot!(lscene1, S_approx;
+        arrowsize_factor,
+        lengthscale_factor
+    )
+    fieldplot!(lscene2, analytic_S_approx;
+        arrowsize_factor,
+        lengthscale_factor
+    )
+    # 2D
+    arrowsize_factor = 0.15
+    lengthscale_factor = 6
+    plt1 = fieldplot!(ax1, S_slice;
+        arrowsize_factor,
+        lengthscale_factor
+    )
+    plt2 = fieldplot!(ax2, analytic_S_slice;
+        arrowsize_factor,
+        lengthscale_factor
+    )
+    # Colorbars
+    Colorbar(fig[2,1][2,1], plt1,
+        width = Relative(3/4), height = 15,
+        vertical = false, flipaxis = false,
+        label = "S ($unit_S)",
+    )
+    Colorbar(fig[2,2][2,1], plt2,
+        width = Relative(3/4), height = 15,
+        vertical = false, flipaxis = false,
+        label = "S ($unit_S)",
+    )
 
-    return plt
+    WGLMakie.trim!(fig.layout)
+    fig
 end

src/plots/utils.jl

@@ -0,0 +1,8 @@
+struct MinMaxTicks end
+
+function MakieLayout.get_tickvalues(::MinMaxTicks, min, max)
+    ticks = MakieLayout.get_tickvalues(WilkinsonTicks(5, k_min = 3), min, max)
+    pushfirst!(ticks, min)
+    push!(ticks, max)
+    ticks
+end

src/units.jl

@@ -2,6 +2,7 @@ function add_λ_units(λ)
     expr = quote
 
         using Unitful
+        using Unitful: MW
         using PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e
 
         local λ = $λ
@@ -12,6 +13,7 @@ function add_λ_units(λ)
         @unit unit_B "mₑ ω/e" UnitMagneticField m_e * ω / e false
         @unit unit_E "mₑ c ω/e" UnitElectricField m_e * ω / e * c_0 false
         @unit unit_L "mₑ c²/ω" UnitAngularMomenta m_e * c_0^2 / ω false
+        @unit unit_S "MW ω²/c²" UnitPoyntingFlux MW * ω^2 / c_0^2 false
         @unit pₑ "mₑ c" ElectronMomenta m_e * c_0 false
 
         Unitful.register(@__MODULE__)