Merge pull request #354 from chrisbrahms/peakpower

Add modal sum option to peak power in Stats and Processing

src/Processing.jl

@@ -311,15 +311,20 @@ end
 fwhm(x::Vector, I::Vector; minmax=:min) = Maths.fwhm(x, I; minmax)
 
 """
-    peakpower(grid, Eω; bandpass=nothing, oversampling=1)
-    peakpower(output; bandpass=nothing, oversampling=1)
+    peakpower(grid, Eω; bandpass=nothing, oversampling=1, sumdims=nothing)
+    peakpower(output; bandpass=nothing, oversampling=1, sumdims=nothing)
 
 Extract the peak power. If `bandpass` is given, bandpass the field according to
-[`window_maybe`](@ref).
+[`window_maybe`](@ref). If `sumdims` is not `nothing`, sum the time-dependent power
+over these dimensions (e.g. modes) before taking the maximum.
 """
-function peakpower(grid, Eω; bandpass=nothing, oversampling=1)
+function peakpower(grid, Eω; bandpass=nothing, oversampling=1, sumdims=nothing)
     to, Eto = getEt(grid, Eω; oversampling=oversampling, bandpass=bandpass)
-    dropdims(maximum(abs2.(Eto); dims=1); dims=1)
+    Pt = abs2.(Eto)
+    if !isnothing(sumdims)
+        Pt = dropdims(sum(Pt; dims=sumdims); dims=sumdims)
+    end
+    dropdims(maximum(Pt; dims=1); dims=1)
 end
 
 function peakpower(output; kwargs...)

src/Stats.jl

@@ -87,8 +87,11 @@ function peakpower(grid)
     function addstat!(d, Eω, Et, z, dz)
         if ndims(Et) > 1
             d["peakpower"] = dropdims(maximum(abs2.(Et), dims=1), dims=1)
+            d["peakpower_allmodes"] = maximum(eachindex(grid.t)) do ii
+                sum(abs2, Et[ii, :])
+            end
         else
-            d["peakpower"] = maximum(abs2.(Et))
+            d["peakpower"] = maximum(abs2, Et)
         end
     end
     return addstat!
@@ -104,14 +107,19 @@ dataset is labeled as `peakpower_[label]`.
 function peakpower(grid, Eω, window::Vector{<:Real}; label)
     Etbuf, analytic! = plan_analytic(grid, Eω) # output buffer and function for inverse FT
     Eωbuf = similar(Eω) # buffer for Eω with window applied
+    Pt = zeros((length(grid.t), size(Eωbuf, 2)))
     key = "peakpower_$label"
     function addstat!(d, Eω, Et, z, dz)
         Eωbuf .= Eω .* window
         analytic!(Etbuf, Eωbuf)
         if ndims(Etbuf) > 1
-            d[key] = dropdims(maximum(abs2.(Etbuf), dims=1), dims=1)
+            Pt .= abs2.(Etbuf)
+            d[key] = dropdims(maximum(Pt, dims=1), dims=1)
+            d[key*"_allmodes"] = maximum(eachindex(grid.t)) do ii
+                sum(Pt[ii, :]; dims=2)
+            end
         else
-            d[key] = maximum(abs2.(Etbuf))
+            d[key] = maximum(abs2, Etbuf)
         end
     end
 end
@@ -141,7 +149,7 @@ Create stats function to calculate the mode-averaged peak intensity given the ef
 """
 function peakintensity(grid, aeff)
     function addstat!(d, Eω, Et, z, dz)
-        d["peakintensity"] = maximum(abs2.(Et))/aeff(z)
+        d["peakintensity"] = maximum(abs2, Et)/aeff(z)
     end
 end
 
@@ -157,9 +165,11 @@ function peakintensity(grid, modes::Modes.ModeCollection; components=:y)
     function addstat!(d, Eω, Et, z, dz)
         Modes.to_space!(Et0, Et, (0, 0), tospace; z=z)
         if npol > 1
-            d["peakintensity"] = c*ε_0/2 * maximum(sum(abs2.(Et0), dims=2))
+            d["peakintensity"] = c*ε_0/2 * maximum(eachindex(grid.t)) do ii
+                sum(abs2, Et0[ii, :])
+            end
         else
-            d["peakintensity"] = c*ε_0/2 * maximum(abs2.(Et0))
+            d["peakintensity"] = c*ε_0/2 * maximum(abs2, Et0)
         end
     end
 end
@@ -200,7 +210,7 @@ function fwhm_r(grid, modes; components=:y)
     function addstat!(d, Eω, Et, z, dz)
         function f(r)
             Modes.to_space!(Eω0, Eω, (r, 0), tospace; z=z)
-            sum(abs2.(Eω0))
+            sum(abs2, Eω0)
         end
         d["fwhm_r"] = 2*Maths.hwhm(f)
     end

test/test_processing.jl

@@ -1,6 +1,7 @@
 import Test: @test, @testset
 import FFTW
-import Luna: Grid, Processing, Maths, Fields, settings, PhysData
+using Luna
+import Luna: settings
 import Luna.PhysData: wlfreq
 import NumericalIntegration: integrate
 
@@ -275,3 +276,21 @@ end
     end
 end
 
+@testset "peak power modal sum" begin
+    λ0 = 800e-9
+    τfwhm = 10e-15
+    pHE11 = 1e8
+    pHE12 = 5e7
+
+    pulses = [
+        Pulses.GaussPulse(;λ0, τfwhm, power=pHE11, mode=:HE11)
+        Pulses.GaussPulse(;λ0, τfwhm, power=pHE12, mode=:HE12)
+    ]
+
+    out = prop_capillary(100e-6, 0.2, :He, 1e-3; λ0, pulses, λlims=(70e-9, 4e-6), trange=1e-12, modes=4)
+
+    @test all(isapprox.(Processing.peakpower(out)[1:2, 1], [pHE11, pHE12]; rtol=1e-4))
+    @test isapprox(Processing.peakpower(out; sumdims=2)[1], pHE11+pHE12; rtol=1e-4)
+
+    @test isapprox(out["stats"]["peakpower_allmodes"][1], pHE11+pHE12; rtol=1e-4)
+end