peak intensity

LupoLab · Apr 14, 2020 · eb2385c · jtravs · Apr 14, 2020 · eb2385c
1 parent 499f70a
commit eb2385c
Show file tree

Hide file tree

Showing 3 changed files with 46 additions and 1 deletion.
diff --git a/src/Stats.jl b/src/Stats.jl
@@ -87,6 +87,37 @@ function peakpower(grid)
     return addstat!
 end
 
+"""
+    peakintensity(grid, mode)
+
+Create stats function to calculate the peak intensity for a single (averaged) mode.
+"""
+function peakintensity(grid, mode::Modes.AbstractMode)
+    fac(z) = sum(abs2.(Modes.Exy(mode, (0, 0); z=z)))
+    function addstat!(d, Eω, Et, z, dz)
+        d["peakintensity"] = c*ε_0/2*fac(z)*maximum(abs2.(Et))
+    end
+end
+
+"""
+    peakintensity(grid, mode)
+
+Create stats function to calculate the peak intensity for a single mode.
+"""
+function peakintensity(grid, modes::NTuple{N, Modes.AbstractMode}; components=:y) where N
+    tospace = Modes.ToSpace(modes, components=components)
+    npol = tospace.npol
+    Et0 = zeros(ComplexF64, (length(grid.t), npol))
+    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))
+        else
+            d["peakintensity"] = c*ε_0/2 * maximum(abs2.(Et0))
+        end
+    end
+end
+
 """
     fwhm_t(grid)
 

diff --git a/test/test_main.jl b/test/test_main.jl
@@ -58,6 +58,7 @@ statsfun = Stats.collect_stats(grid, Eω,
                                Stats.energy(grid, energyfunω),
                                Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"),
                                Stats.peakpower(grid),
+                               Stats.peakintensity(grid, m),
                                Stats.fwhm_t(grid),
                                Stats.electrondensity(grid, ionrate, densityfun, aeff),
                                Stats.density(densityfun))
@@ -118,6 +119,12 @@ plt.plot(output["stats"]["z"].*1e2, output["stats"]["peakpower"].*1e-9)
 plt.xlabel("Distance (cm)")
 plt.ylabel("Peak power (GW)")
 
+##
+plt.figure()
+plt.plot(output["stats"]["z"].*1e2, output["stats"]["peakintensity"].*1e-4.*1e-12)
+plt.xlabel("Distance (cm)")
+plt.ylabel("Peak intensity (TW/cm\$^2\$)")
+
 ##
 plt.figure()
 plt.plot(output["stats"]["z"].*1e2, output["stats"]["fwhm_t_min"].*1e15)

diff --git a/test/test_modal.jl b/test/test_modal.jl
@@ -54,14 +54,15 @@ statsfun = Stats.collect_stats(grid, Eω,
                                Stats.energy(grid, energyfunω),
                                Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"),
                                Stats.peakpower(grid),
+                               Stats.peakintensity(grid, modes),
                                Stats.fwhm_t(grid),
                                Stats.fwhm_r(grid, modes),
                                Stats.electrondensity(grid, ionrate, densityfun, modes)
                                )
 output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun)
 linop = LinearOps.make_const_linop(grid, modes, λ0)
 
-Luna.run(Eω, grid, linop, transform, FT, output)
+Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)
 
 ω = grid.ω
 t = grid.t
@@ -119,6 +120,12 @@ plt.xlabel("Distance (cm)")
 plt.ylabel("Peak power (GW)")
 plt.legend()
 
+##
+plt.figure()
+plt.plot(output["stats"]["z"].*1e2, output["stats"]["peakintensity"].*1e-4.*1e-12)
+plt.xlabel("Distance (cm)")
+plt.ylabel("Peak intensity (TW/cm\$^2\$)")
+
 ##
 plt.figure()
 for i = 1:length(modes)