diff --git a/.gitignore b/.gitignore index b1234cda..56651bc0 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ deps/deps.jl Manifest.toml _git2_* docs/build +*.info diff --git a/docs/src/scans.md b/docs/src/scans.md index da95c764..7d0aa9d8 100644 --- a/docs/src/scans.md +++ b/docs/src/scans.md @@ -5,11 +5,7 @@ The first step is to write a standard Luna script to run a *single* simulation i ```julia # Lots of setup work... gausspulse = let τ=30e-15, λ=λ0 - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end + end in1 = (func=gausspulse, energy=1e-6) inputs = (in1, ) @@ -38,15 +34,7 @@ Then, make the following changes: @scan begin # Lots of setup work... -gausspulse = let τ=$τ, λ=λ0 - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end -end -in1 = (func=gausspulse, energy=$energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=$τ, energy=$energy) # More setup work... Luna.run(...) # run the simulation end diff --git a/examples/MI_modeAvg_env.jl b/examples/MI_modeAvg_env.jl index 174cf691..690403ca 100644 --- a/examples/MI_modeAvg_env.jl +++ b/examples/MI_modeAvg_env.jl @@ -1,31 +1,20 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 15e-6 gas = :Ar pres = 25 -τ = 600e-15 +τfwhm = 600e-15 λ0 = 800e-9 +flength = 80e-2 +energy = 10e-6 -grid = Grid.EnvGrid(80e-2, 800e-9, (220e-9, 3000e-9), 4e-12) +grid = Grid.EnvGrid(flength, λ0, (220e-9, 3000e-9), 4e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -40,16 +29,22 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) # Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=10e-6) -inputs = (in1, ) +inputs = (Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy), Fields.ShotNoise()) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) -Luna.shotnoise!(Eω, grid) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, + Stats.ω0(grid), + Stats.energy(grid, energyfunω), + Stats.peakpower(grid), + Stats.fwhm_t(grid), + Stats.density(densityfun)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = grid.ω t = grid.t f = FFTW.fftshift(ω, 1)./2π.*1e-15 @@ -68,7 +63,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end pygui(true) diff --git a/examples/Raman/Raman_modeAvg.jl b/examples/Raman/Raman_modeAvg.jl index 287280a7..1efcd4b7 100644 --- a/examples/Raman/Raman_modeAvg.jl +++ b/examples/Raman/Raman_modeAvg.jl @@ -1,28 +1,20 @@ using Luna -import Logging -import FFTW -Logging.disable_logging(Logging.BelowMinLevel) a = 13e-6 gas = :H2 pres = 5 +flength = 200e-2 -τ = 20e-15 +τfwhm = 20e-15 λ0 = 800e-9 energy = 1e-6 -grid = Grid.RealGrid(200e-2, 800e-9, (180e-9, 3000e-9), 4e-12) +grid = Grid.RealGrid(flength, λ0, (180e-9, 3000e-9), 4e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -34,10 +26,9 @@ linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -46,10 +37,12 @@ statsfun = Stats.collect_stats(grid, Eω, Stats.peakintensity(grid, aeff), Stats.fwhm_t(grid), Stats.density(densityfun)) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t @@ -69,7 +62,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end import PyPlot:pygui, plt diff --git a/examples/Raman/Raman_modeAvg_env.jl b/examples/Raman/Raman_modeAvg_env.jl index 66553560..6492ad5b 100644 --- a/examples/Raman/Raman_modeAvg_env.jl +++ b/examples/Raman/Raman_modeAvg_env.jl @@ -1,27 +1,20 @@ using Luna -import Logging -import FFTW -Logging.disable_logging(Logging.BelowMinLevel) a = 13e-6 gas = :H2 pres = 5 +flength = 200e-2 -τ = 20e-15 +τfwhm = 20e-15 λ0 = 800e-9 energy = 1e-6 -grid = Grid.EnvGrid(200e-2, 800e-9, (180e-9, 3000e-9), 4e-12) +grid = Grid.EnvGrid(flength, λ0, (180e-9, 3000e-9), 4e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - Et = @. sqrt(It) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -33,22 +26,22 @@ linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω), - Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"), Stats.peakpower(grid), Stats.fwhm_t(grid), Stats.density(densityfun)) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t f = FFTW.fftshift(ω, 1)./2π.*1e-15 @@ -68,7 +61,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end import PyPlot:pygui, plt diff --git a/examples/Raman/Raman_modeAvg_noTHG.jl b/examples/Raman/Raman_modeAvg_noTHG.jl index 7497639c..f8112d65 100644 --- a/examples/Raman/Raman_modeAvg_noTHG.jl +++ b/examples/Raman/Raman_modeAvg_noTHG.jl @@ -1,29 +1,20 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes, Raman -import Logging -import FFTW -Logging.disable_logging(Logging.BelowMinLevel) +using Luna a = 13e-6 gas = :H2 pres = 5 +flength = 200e-2 -τ = 20e-15 +τfwhm = 20e-15 λ0 = 800e-9 energy = 1e-6 -grid = Grid.RealGrid(200e-2, 800e-9, (180e-9, 3000e-9), 4e-12) +grid = Grid.RealGrid(flength, λ0, (180e-9, 3000e-9), 4e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -35,16 +26,17 @@ linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t @@ -64,7 +56,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end import PyPlot:pygui, plt diff --git a/examples/arpcf_modeAvg.jl b/examples/arpcf_modeAvg.jl index 8e388790..ef52546d 100644 --- a/examples/arpcf_modeAvg.jl +++ b/examples/arpcf_modeAvg.jl @@ -1,34 +1,22 @@ using Luna -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -# import DSP.Unwrap: unwrap - -import PyPlot:pygui, plt a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) m = Antiresonant.ZeisbergerMode(a, gas, pres, wallthickness=200e-9, loss=false) aeff = let m=m z -> Modes.Aeff(m, z=z) end -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -44,11 +32,10 @@ linop, βfun, β1, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -59,9 +46,13 @@ statsfun = Stats.collect_stats(grid, Eω, Stats.fwhm_t(grid), Stats.electrondensity(grid, ionrate, densityfun, aeff), Stats.density(densityfun)) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) + +import FFTW +import PyPlot:pygui, plt + ## ω = grid.ω t = grid.t @@ -80,8 +71,8 @@ Itlog = log10.(Maths.normbymax(It)) zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) -energy = [energyfun(t, Etout[:, ii]) for ii=1:size(Etout, 2)] -energyω = [energyfunω(ω, Eout[:, ii]) for ii=1:size(Eout, 2)] +energy = [energyfun(Etout[:, ii]) for ii=1:size(Etout, 2)] +energyω = [energyfunω(Eout[:, ii]) for ii=1:size(Eout, 2)] pygui(true) ## diff --git a/examples/basic_modal.jl b/examples/basic_modal.jl index 90daf847..bdb79834 100644 --- a/examples/basic_modal.jl +++ b/examples/basic_modal.jl @@ -1,38 +1,24 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 modes = ( Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), Capillary.MarcatilliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false) ) -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -42,10 +28,9 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) # Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=false) statsfun = Stats.collect_stats(grid, Eω, @@ -58,11 +43,14 @@ statsfun = Stats.collect_stats(grid, Eω, Stats.fwhm_r(grid, modes), Stats.electrondensity(grid, ionrate, densityfun, modes) ) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output, status_period=5) +import FFTW +import PyPlot:pygui, plt + ω = grid.ω t = grid.t diff --git a/examples/basic_modal_env.jl b/examples/basic_modal_env.jl index fdf66fbd..6813d111 100644 --- a/examples/basic_modal_env.jl +++ b/examples/basic_modal_env.jl @@ -1,19 +1,13 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 modes = ( Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), @@ -21,17 +15,11 @@ modes = ( ) nmodes = length(modes) -grid = Grid.EnvGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.EnvGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -41,18 +29,26 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=false) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, + Stats.ω0(grid), + Stats.energy(grid, energyfunω), + Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"), + Stats.peakpower(grid), + Stats.fwhm_t(grid), + Stats.density(densityfun)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = FFTW.fftshift(grid.ω) t = grid.t diff --git a/examples/basic_modeAvg.jl b/examples/basic_modeAvg.jl index ed2f2ae4..6ecd7cad 100644 --- a/examples/basic_modeAvg.jl +++ b/examples/basic_modeAvg.jl @@ -1,29 +1,22 @@ using Luna -import Logging -import FFTW -Logging.disable_logging(Logging.BelowMinLevel) a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff = let m=m z -> Modes.Aeff(m, z=z) end -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -39,11 +32,10 @@ linop, βfun, β1, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -59,6 +51,7 @@ output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) import PyPlot:pygui, plt +import FFTW ω = grid.ω t = grid.t @@ -77,8 +70,8 @@ Itlog = log10.(Maths.normbymax(It)) zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) -energy = [energyfun(t, Etout[:, ii]) for ii=1:size(Etout, 2)] -energyω = [energyfunω(ω, Eout[:, ii]) for ii=1:size(Eout, 2)] +energy = [energyfun(Etout[:, ii]) for ii=1:size(Etout, 2)] +energyω = [energyfunω(Eout[:, ii]) for ii=1:size(Eout, 2)] pygui(true) ## diff --git a/examples/basic_modeAvg_env.jl b/examples/basic_modeAvg_env.jl index 15aa7517..571a8bca 100644 --- a/examples/basic_modeAvg_env.jl +++ b/examples/basic_modeAvg_env.jl @@ -1,29 +1,22 @@ using Luna -import Logging -import FFTW -Logging.disable_logging(Logging.BelowMinLevel) a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 -grid = Grid.EnvGrid(15e-2, 800e-9, (160e-9, 3000e-9), 5e-12) +grid = Grid.EnvGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false, model=:full) aeff = let m=m z -> Modes.Aeff(m, z=z) end -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -38,11 +31,10 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) # Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -51,11 +43,12 @@ statsfun = Stats.collect_stats(grid, Eω, Stats.peakpower(grid), Stats.fwhm_t(grid), Stats.density(densityfun)) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) import PyPlot:pygui, plt +import FFTW ω = grid.ω t = grid.t @@ -75,9 +68,9 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end -energyω = [energyfunω(ω, Eout[:, ii]) for ii=1:size(Eout, 2)] +energyω = [energyfunω(Eout[:, ii]) for ii=1:size(Eout, 2)] pygui(true) plt.figure() diff --git a/examples/basic_modeAvg_env_THG.jl b/examples/basic_modeAvg_env_THG.jl index cb3908f7..0b9afde5 100644 --- a/examples/basic_modeAvg_env_THG.jl +++ b/examples/basic_modeAvg_env_THG.jl @@ -1,34 +1,20 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - -import PyPlot: pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +flength = 0.5e-2 +energy = 1e-6 -grid = Grid.EnvGrid(0.5e-2, 800e-9, (160e-9, 3000e-9), 1e-12, thg=true) -#grid = Grid.EnvGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12, thg=true) +grid = Grid.EnvGrid(flength, λ0, (160e-9, 3000e-9), 1e-12, thg=true) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -39,16 +25,24 @@ normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) responses = (Nonlinear.Kerr_env_thg(PhysData.γ3_gas(gas), 2π*PhysData.c/λ0, grid.to),) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, + Stats.ω0(grid), + Stats.energy(grid, energyfunω), + Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"), + Stats.peakpower(grid), + Stats.fwhm_t(grid), + Stats.density(densityfun)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot: pygui, plt + ω = grid.ω t = grid.t @@ -67,7 +61,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Etout energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end pygui(true) diff --git a/examples/basic_modeAvg_field_noTHG.jl b/examples/basic_modeAvg_field_noTHG.jl index bed4599f..a9ddbadd 100644 --- a/examples/basic_modeAvg_field_noTHG.jl +++ b/examples/basic_modeAvg_field_noTHG.jl @@ -1,28 +1,20 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes -import FFTW -import Logging -Logging.disable_logging(Logging.BelowMinLevel) +using Luna a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -37,13 +29,20 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field_nothg(PhysData.γ3_gas(gas),length(grid.to)), Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, + Stats.ω0(grid), + Stats.energy(grid, energyfunω), + Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"), + Stats.peakpower(grid), + Stats.peakintensity(grid, aeff), + Stats.fwhm_t(grid), + Stats.electrondensity(grid, ionrate, densityfun, aeff), + Stats.density(densityfun)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) @@ -53,6 +52,8 @@ t = grid.t zout = output.data["z"] Eout = output.data["Eω"] +import FFTW + Etout = FFTW.irfft(Eout, length(grid.t), 1) Ilog = log10.(Maths.normbymax(abs2.(Eout))) @@ -66,7 +67,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end import PyPlot:pygui, plt diff --git a/examples/basic_scan.jl b/examples/basic_scan.jl index 37741cab..c0ee5947 100644 --- a/examples/basic_scan.jl +++ b/examples/basic_scan.jl @@ -1,5 +1,4 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps +using Luna import Luna.Scans: @scanvar, @scan, @scaninit import Logging: @info @@ -8,7 +7,6 @@ import Logging: @info @scanvar energy = range(0.1e-6, 1.5e-6, length=16) @scanvar τ = range(25e-15, 35e-15, length=11) - @scan begin a = 13e-6 gas = :Ar @@ -20,15 +18,7 @@ grid = Grid.RealGrid(1e-2, 800e-9, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) -energyfun = NonlinearRHS.energy_mode_avg(m) - -gausspulse = let τ=$τ, λ=λ0 - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end -end +energyfun, energyfunω = Fields.energyfuncs(grid) println("τ: $($τ * 1e15)") println("E: $($energy * 1e6)") @@ -46,12 +36,11 @@ linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun) -in1 = (func=gausspulse, energy=$energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=$τ, energy=$energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) output = Output.@ScanHDF5Output(0, grid.zmax, 101, (length(grid.ω),), statsfun) println(output["meta"]["scanvars"]) diff --git a/examples/freespace/full3D.jl b/examples/freespace/full3D.jl index 5ace2cd5..7d014e6b 100644 --- a/examples/freespace/full3D.jl +++ b/examples/freespace/full3D.jl @@ -1,13 +1,9 @@ -import Luna -import Luna: Grid, Maths, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Hankel -import Logging +using Luna +Luna.set_fftw_mode(:estimate) import FFTW import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) import Luna.PhysData: wlfreq -import PyPlot:pygui, plt - gas = :Ar pres = 4 @@ -26,13 +22,7 @@ xygrid = Grid.FreeGrid(R, N) x = xygrid.x y = xygrid.y -energyfun, energyfunω = NonlinearRHS.energy_free(grid, xygrid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) .* Maths.gauss.(xygrid.r, w0/2) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid, xygrid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -45,13 +35,12 @@ responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) linop = LinearOps.make_const_linop(grid, xygrid, PhysData.ref_index_fun(gas, pres)) normfun = NonlinearRHS.const_norm_free(grid, xygrid, PhysData.ref_index_fun(gas, pres)) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0) -Eω, transform, FT = Luna.setup(grid, xygrid, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, xygrid, densityfun, normfun, responses, inputs) -# statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 21, (length(grid.ω), N, N)) +# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 21) Luna.run(Eω, grid, linop, transform, FT, output, max_dz=Inf, init_dz=1e-1) @@ -73,7 +62,7 @@ Iωyx = abs2.(Eωyx); Iyx = zeros(Float64, (length(y), length(x), length(zout))); energy = zeros(length(zout)); for ii = 1:size(Etyx, 4) - energy[ii] = energyfun(t, Etyx[:, :, :, ii]); + energy[ii] = energyfun(Etyx[:, :, :, ii]); Iyx[:, :, ii] = (grid.ω[2]-grid.ω[1]) .* sum(Iωyx[:, :, :, ii], dims=1); end @@ -84,6 +73,7 @@ E0ωyx = FFTW.ifft(Eω[ω0idx, :, :], (1, 2)); Iωyx = abs2.(Eωyx) Iωyxlog = log10.(Maths.normbymax(Iωyx)); +import PyPlot:pygui, plt pygui(true) plt.figure() plt.pcolormesh(x, y, (abs2.(E0ωyx))) diff --git a/examples/freespace/full3D_env.jl b/examples/freespace/full3D_env.jl index b89b8094..8ed40ab8 100644 --- a/examples/freespace/full3D_env.jl +++ b/examples/freespace/full3D_env.jl @@ -1,11 +1,7 @@ -import Luna -import Luna: Grid, Maths, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps -import Logging +using Luna +Luna.set_fftw_mode(:estimate) import FFTW import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt gas = :Ar pres = 4 @@ -21,44 +17,28 @@ R = 6e-3 N = 128 grid = Grid.EnvGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12) +xygrid = Grid.FreeGrid(R, N) -Dx = 2R/N -n = collect(range(0, length=N)) -x = @. (n-N/2) * Dx -y = copy(x) - -r = sqrt.(reshape(x, (1, 1, N)).^2 .+ reshape(y, (1, N)).^2) - -xr = Array{ComplexF64}(undef, length(grid.t), length(y), length(x)) -FT = FFTW.plan_fft(xr, (1, 2, 3), flags=FFTW.MEASURE) - -energyfun = NonlinearRHS.energy_free_env(x, y) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) .* Maths.gauss.(r, w0/2) - Et = @. sqrt(It) -end +x = xygrid.x +y = xygrid.y +energyfun, energyfunω = Fields.energyfuncs(grid, xygrid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 -responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) +responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) -linop = LinearOps.make_const_linop(grid, x, y, PhysData.ref_index_fun(gas, pres)) -normfun = NonlinearRHS.norm_free(grid.ω, x, y, PhysData.ref_index_fun(gas, pres)) +linop = LinearOps.make_const_linop(grid, xygrid, PhysData.ref_index_fun(gas, pres)) +normfun = NonlinearRHS.const_norm_free(grid, xygrid, PhysData.ref_index_fun(gas, pres)) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0) -Eω, transform, FTo = Luna.setup(grid, FT, x, y, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, xygrid, densityfun, normfun, responses, inputs) -# statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 21, (length(grid.ω), N, N)) +# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 21) -xwin = Maths.planck_taper(x, -R, -0.8R, 0.8R, R) -xywin = reshape(xwin, (1, length(xwin))) .* reshape(xwin, (1, 1, length(xwin))) - -Luna.run(Eω, grid, linop, transform, FT, output, max_dz=Inf, init_dz=0.1) +Luna.run(Eω, grid, linop, transform, FT, output, max_dz=Inf, init_dz=1e-1) ω = FFTW.fftshift(grid.ω) t = grid.t @@ -80,7 +60,7 @@ Iωyx = abs2.(Eωyx); Iyx = zeros(Float64, (length(y), length(x), length(zout))); energy = zeros(length(zout)); for ii = 1:size(Etyx, 4) - energy[ii] = energyfun(t, Etyx[:, :, :, ii]); + energy[ii] = energyfun(Etyx[:, :, :, ii]); Iyx[:, :, ii] = (grid.ω[2]-grid.ω[1]) .* sum(Iωyx[:, :, :, ii], dims=1); end @@ -91,6 +71,7 @@ E0ωyx = FFTW.ifft(Eω[ω0idx, :, :], (1, 2)); Iωyx = abs2.(Eωyx); Iωyxlog = log10.(Maths.normbymax(Iωyx)); +import PyPlot:pygui, plt pygui(true) plt.figure() plt.pcolormesh(x, y, (abs2.(E0ωyx))) diff --git a/examples/freespace/radial.jl b/examples/freespace/radial.jl index f64ddc28..f751c06c 100644 --- a/examples/freespace/radial.jl +++ b/examples/freespace/radial.jl @@ -1,13 +1,8 @@ -import Luna -import Luna: Grid, Maths, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Plotting +using Luna import Luna.PhysData: wlfreq -import Logging import FFTW import Hankel import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt gas = :Ar pres = 1.2 @@ -25,26 +20,7 @@ N = 1024 grid = Grid.RealGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12) q = Hankel.QDHT(R, N, dim=2) -energyfun, energyfun_ω = NonlinearRHS.energy_radial(grid, q) - -function prop(E, z) - Eω = FFTW.rfft(E, 1) - Eωk = q * Eω - kzsq = @. (grid.ω/PhysData.c)^2 - (q.k^2)' - kzsq[kzsq .< 0] .= 0 - kz = sqrt.(kzsq) - @. Eωk *= exp(-1im * z * (kz - grid.ω/PhysData.c)) - Eω = q \ Eωk - E = FFTW.irfft(Eω, length(grid.t), 1) - return E -end - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) .* Maths.gauss(q.r, w0/2)' - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - return prop(Et, -0.3) -end +energyfun, energyfun_ω = Fields.energyfuncs(grid, q) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -52,20 +28,19 @@ densityfun(z) = dens0 ionpot = PhysData.ionisation_potential(gas) ionrate = Ionisation.ionrate_fun!_PPTcached(gas, λ0) -responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)), - Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) +responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) + #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) linop = LinearOps.make_const_linop(grid, q, PhysData.ref_index_fun(gas, pres)) normfun = NonlinearRHS.const_norm_radial(grid, q, PhysData.ref_index_fun(gas, pres)) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0, propz=-0.3) -Eω, transform, FT = Luna.setup(grid, q, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, q, densityfun, normfun, responses, inputs) -# statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω), length(q.r))) +# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201) Luna.run(Eω, grid, linop, transform, FT, output) ω = grid.ω @@ -92,7 +67,7 @@ Ir = zeros(Float64, (length(q.r), length(zout))) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 3) - energy[ii] = energyfun(t, Etout[:, :, ii]); + energy[ii] = energyfun(Etout[:, :, ii]); Ir[:, ii] = integrate(grid.ω, Iωr[:, :, ii], SimpsonEven()); end @@ -105,12 +80,13 @@ idcs = [argmin(abs.(zout .- point)) for point in points] Epoints = [Hankel.symmetric(Et[:, :, idxi], q) for idxi in idcs] rsym = Hankel.Rsymmetric(q); +import PyPlot:pygui, plt pygui(true) Iλ0 = Iωr[ω0idx, :, :] λ0 = 2π*PhysData.c/grid.ω[ω0idx] w1 = w0*sqrt(1+(L/2*λ0/(π*w0^2))^2) # w1 = w0 -Iλ0_analytic = Maths.gauss(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx) +Iλ0_analytic = Maths.gauss.(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx) plt.figure() plt.plot(q.r*1e3, Maths.normbymax(Iλ0[:, end])) plt.plot(q.r*1e3, Maths.normbymax(Iλ0_analytic), "--") diff --git a/examples/freespace/radial_env.jl b/examples/freespace/radial_env.jl index b322f8db..25949f8e 100644 --- a/examples/freespace/radial_env.jl +++ b/examples/freespace/radial_env.jl @@ -1,12 +1,7 @@ -import Luna -import Luna: Grid, Maths, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Plotting -import Logging +using Luna import FFTW import Hankel import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt gas = :Ar pres = 1.2 @@ -24,25 +19,7 @@ N = 512 grid = Grid.EnvGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12) q = Hankel.QDHT(R, N, dim=2) -energyfun, energyfun_ω = NonlinearRHS.energy_radial(grid, q) - -function prop(E, z) - Eω = FFTW.fft(E, 1) - Eωk = q * Eω - kzsq = @. (grid.ω/PhysData.c)^2 - (q.k^2)' - kzsq[kzsq .< 0] .= 0 - kz = sqrt.(kzsq) - @. Eωk *= exp(-1im * z * (kz - grid.ω/PhysData.c)) - Eω = q \ Eωk - E = FFTW.ifft(Eω, 1) - return E -end - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) .* Maths.gauss(q.r, w0/2)' - Et = @. sqrt(It) - return prop(Et, -0.3) -end +energyfun, energyfun_ω = Fields.energyfuncs(grid, q) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -51,19 +28,14 @@ responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) linop = LinearOps.make_const_linop(grid, q, PhysData.ref_index_fun(gas, pres)) -normfun = NonlinearRHS.norm_radial(grid.ω, q, PhysData.ref_index_fun(gas, pres)) +normfun = NonlinearRHS.const_norm_radial(grid, q, PhysData.ref_index_fun(gas, pres)) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0, propz=-0.3) -Eω, transform, FT = Luna.setup(grid, q, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, q, densityfun, normfun, responses, inputs) -# statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω), length(q.r))) -Ert = FT \ (q \ Eω) -println(energyfun(grid.t, Ert)) -println(energyfun_ω(grid.ω, Eω)) -error() +# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201) Luna.run(Eω, grid, linop, transform, FT, output) @@ -90,7 +62,7 @@ Ir = zeros(Float64, (length(q.r), length(zout))) energy = zeros(length(zout)) for ii = 1:size(Etout, 3) - energy[ii] = energyfun(t, Etout[:, :, ii]); + energy[ii] = energyfun(Etout[:, :, ii]); Ir[:, ii] = integrate(ω, Iωr[:, :, ii], SimpsonEven()); end @@ -103,11 +75,12 @@ idcs = [argmin(abs.(zout .- point)) for point in points] Epoints = [Hankel.symmetric(Etout[:, :, idxi], q) for idxi in idcs] rsym = Hankel.Rsymmetric(q); +import PyPlot:pygui, plt pygui(true) Iλ0 = Iωr[ω0idx, :, :] w1 = w0*sqrt(1+(L/2*λ0/(π*w0^2))^2) # w1 = w0 -Iλ0_analytic = Maths.gauss(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx) +Iλ0_analytic = Maths.gauss.(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx) plt.figure() plt.plot(q.r*1e3, Maths.normbymax(Iλ0[:, end])) plt.plot(q.r*1e3, Maths.normbymax(Iλ0_analytic), "--") diff --git a/examples/freespace/radial_hcf.jl b/examples/freespace/radial_hcf.jl new file mode 100644 index 00000000..14fe0800 --- /dev/null +++ b/examples/freespace/radial_hcf.jl @@ -0,0 +1,135 @@ +using Luna +import Luna.PhysData: wlfreq +import FFTW +import Hankel +import NumericalIntegration: integrate, SimpsonEven + +a = 15e-6 +gas = :Ar +pres = 5.0 + +τ = 30e-15 +λ0 = 800e-9 + +w0 = 0.64*a +energy = 1e-6 +L = 0.15 + +R = a +N = 32 + +grid = Grid.RealGrid(L, λ0, (200e-9, 3000e-9), 0.6e-12) +q = Hankel.QDHT(R, N, dim=2) + +energyfun, energyfun_ω = Fields.energyfuncs(grid, q) + +dens0 = PhysData.density(gas, pres) +densityfun(z) = dens0 + +ionpot = PhysData.ionisation_potential(gas) +ionrate = Ionisation.ionrate_fun!_PPTcached(gas, λ0) + +responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) + #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) + +linop = LinearOps.make_const_linop(grid, q, PhysData.ref_index_fun(gas, pres)) + +normfun = NonlinearRHS.const_norm_radial(grid, q, PhysData.ref_index_fun(gas, pres)) + +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0) + +Eω, transform, FT = Luna.setup(grid, q, densityfun, normfun, responses, inputs) + +# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +#output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω), length(q.r))) +output = Output.MemoryOutput(0, grid.zmax, 201) +Luna.run(Eω, grid, linop, transform, FT, output) + +ω = grid.ω +t = grid.t + +zout = output.data["z"] +Eout = output.data["Eω"] + +Erout = (q \ Eout) +Iωr = abs2.(Erout) +# Iω0 = Iωr[:, 1, :] +Er0 = dropdims(Hankel.onaxis(Eout, q), dims=2); +Iω0 = abs2.(Er0); +Iω0log = log10.(Maths.normbymax(Iω0)); +Etout = FFTW.irfft(Erout, length(grid.t), 1) + +Ilog = log10.(Maths.normbymax(abs2.(Eout))) + +It = PhysData.c * PhysData.ε_0/2 * abs2.(Maths.hilbert(Etout)); +Itlog = log10.(Maths.normbymax(It)) + +Ir = zeros(Float64, (length(q.r), length(zout))) + +Et = Maths.hilbert(Etout) +energy = zeros(length(zout)) +for ii = 1:size(Etout, 3) + energy[ii] = energyfun(Etout[:, :, ii]); + Ir[:, ii] = integrate(grid.ω, Iωr[:, :, ii], SimpsonEven()); +end + +ω0idx = argmin(abs.(grid.ω .- 2π*PhysData.c/λ0)) + +zr = π*w0^2/λ0 +points = L/2 .+ [-15, 3, 21].*zr +idcs = [argmin(abs.(zout .- point)) for point in points] + +Epoints = [Hankel.symmetric(Et[:, :, idxi], q) for idxi in idcs] +rsym = Hankel.Rsymmetric(q); + +import PyPlot:pygui, plt +pygui(true) + +plt.figure() +plt.pcolormesh(zout*1e2, q.r*1e3, Iωr[ω0idx, :, :]) +plt.colorbar() +plt.ylim(0, 4) +plt.xlabel("z (cm)") +plt.ylabel("r (mm)") +plt.title("I(r, ω=ω0)") + +plt.figure() +plt.pcolormesh(zout*1e2, q.r*1e3, It[length(grid.t)÷2, :, :]) +plt.colorbar() +plt.ylim(0, 4) +plt.xlabel("z (cm)") +plt.ylabel("r (mm)") +plt.title("I(r, t=0)") + +plt.figure() +plt.pcolormesh(zout*1e2, q.r*1e3, Ir) +# plt.pcolormesh(zout*1e2, q.r*1e3, log10.(Maths.normbymax(Ir))) +plt.colorbar() +plt.ylim(0, R*1e3) +# plt.clim(-6, 0) +plt.xlabel("z (cm)") +plt.ylabel("r (mm)") +plt.title("\$\\int I(r, \\omega) d\\omega\$") + +plt.figure() +plt.pcolormesh(zout*1e2, grid.ω*1e-15/2π, Iω0log) +plt.colorbar() +plt.clim(0, -6) +plt.xlabel("z (cm)") +plt.ylabel("f (PHz)") +plt.title("I(r=0, ω)") + +plt.figure() +plt.plot(zout.*1e2, energy.*1e6) +plt.xlabel("Distance [cm]") +plt.ylabel("Energy [μJ]") + +jw = Plotting.cmap_white("jet", 512, 10) +fig = plt.figure() +fig.set_size_inches(12, 4) +for ii in 1:3 + plt.subplot(1, 3, ii) + plt.pcolormesh(grid.t*1e15, rsym*1e3, abs2.(Epoints[ii]'), cmap=jw) + plt.xlim(-42, 42) + plt.ylim(-1.8, 1.8) +end diff --git a/examples/full_modal/basic_modal_full.jl b/examples/full_modal/basic_modal_full.jl index 0d8290a8..397ee1cb 100644 --- a/examples/full_modal/basic_modal_full.jl +++ b/examples/full_modal/basic_modal_full.jl @@ -1,41 +1,25 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy=1e-6 modes = ( - Modes.@delegated(Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0), - α=ω->0), - Modes.@delegated(Capillary.MarcatilliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0), - α=ω->0) + Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), + Capillary.MarcatilliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false), ) nmodes = length(modes) -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -45,18 +29,20 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=true) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = grid.ω t = grid.t diff --git a/examples/full_modal/basic_modal_full_bothpolarisations.jl b/examples/full_modal/basic_modal_full_bothpolarisations.jl index a835cb03..e9897eb8 100644 --- a/examples/full_modal/basic_modal_full_bothpolarisations.jl +++ b/examples/full_modal/basic_modal_full_bothpolarisations.jl @@ -1,41 +1,25 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 +flength = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy=1e-6 modes = ( - Modes.@delegated(Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0), - α=ω->0), - Modes.@delegated(Capillary.MarcatilliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0), - α=ω->0) + Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), + Capillary.MarcatilliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false), ) nmodes = length(modes) -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -45,18 +29,29 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :xy; full=true) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, + Stats.ω0(grid), + 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, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = grid.ω t = grid.t diff --git a/examples/gradients/gradient_modal.jl b/examples/gradients/gradient_modal.jl index 132c69a4..6ef6f06d 100644 --- a/examples/gradients/gradient_modal.jl +++ b/examples/gradients/gradient_modal.jl @@ -1,19 +1,12 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 L = 15e-2 @@ -24,35 +17,31 @@ modes = ( Capillary.MarcatilliMode(a, coren, n=1, m=2, kind=:HE, ϕ=0.0, loss=false) ) -grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - ionpot = PhysData.ionisation_potential(gas) ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=false) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = grid.ω t = grid.t diff --git a/examples/gradients/gradient_modal_env.jl b/examples/gradients/gradient_modal_env.jl index bff175d1..e80dc945 100644 --- a/examples/gradients/gradient_modal_env.jl +++ b/examples/gradients/gradient_modal_env.jl @@ -1,21 +1,13 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 +L = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 - -L = 15e-2 +energy = 1e-6 coren, densityfun = Capillary.gradient(gas, L, pres, pres); @@ -25,31 +17,27 @@ modes = ( ) nmodes = length(modes) -grid = Grid.EnvGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_env_modal() +energyfun = Fields.energyfuncs(grid)[1] normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end - responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=false) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = FFTW.fftshift(grid.ω) t = grid.t diff --git a/examples/gradients/gradient_modeAvg.jl b/examples/gradients/gradient_modeAvg.jl index c5008a96..39220e59 100644 --- a/examples/gradients/gradient_modeAvg.jl +++ b/examples/gradients/gradient_modeAvg.jl @@ -1,36 +1,22 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 L = 15e-2 -grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) coren, densityfun = Capillary.gradient(gas, L, pres, pres); m = Capillary.MarcatilliMode(a, coren, loss=false, model=:full); +aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_mode_avg(m) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - -# dens0 = PhysData.density(gas, pres) -# densityfun(z) = dens0 +energyfun, energyfunω = Fields.energyfuncs(grid) ionpot = PhysData.ionisation_potential(gas) ionrate = Ionisation.ionrate_fun!_ADK(ionpot) @@ -40,18 +26,20 @@ responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)), linop, βfun = LinearOps.make_linop(grid, m, λ0); -normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun) +normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = grid.ω t = grid.t @@ -71,7 +59,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end pygui(true) diff --git a/examples/gradients/gradient_modeAvg_env.jl b/examples/gradients/gradient_modeAvg_env.jl index 3085fc33..09c7d513 100644 --- a/examples/gradients/gradient_modeAvg_env.jl +++ b/examples/gradients/gradient_modeAvg_env.jl @@ -1,33 +1,22 @@ using Luna -import Logging -import FFTW -Logging.disable_logging(Logging.BelowMinLevel) a = 13e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 L = 15e-2 -grid = Grid.EnvGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 1e-12) coren, densityfun = Capillary.gradient(gas, L, pres, pres); m = Capillary.MarcatilliMode(a, coren, loss=false); aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end - -# dens0 = PhysData.density(gas, pres) -# densityfun(z) = dens0 +energyfun, energyfunω = Fields.energyfuncs(grid) ionpot = PhysData.ionisation_potential(gas) ionrate = Ionisation.ionrate_fun!_ADK(ionpot) @@ -38,10 +27,9 @@ linop, βfun = LinearOps.make_linop(grid, m, λ0); normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -49,10 +37,11 @@ statsfun = Stats.collect_stats(grid, Eω, Stats.peakpower(grid), Stats.fwhm_t(grid), Stats.density(densityfun)) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW import PyPlot:pygui, plt ω = grid.ω t = grid.t @@ -72,7 +61,7 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end pygui(true) diff --git a/examples/polarisation/elliptical.jl b/examples/polarisation/elliptical.jl index ef1f39c2..9cc7e8fd 100644 --- a/examples/polarisation/elliptical.jl +++ b/examples/polarisation/elliptical.jl @@ -22,17 +22,13 @@ nmodes = length(modes) grid = Grid.RealGrid(250e-2, λ0, (200e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end + Etlin = gausspulse(grid.t) -cenergy = energyfun(grid.t, Etlin) +cenergy = energyfun(Etlin) Etlin = sqrt(energy)/sqrt(cenergy) .* Etlin Et = [Etlin zero(grid.t)] @@ -60,16 +56,15 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :xy; full=false) Eω .= Ew -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) diff --git a/examples/polarisation/elliptical_env.jl b/examples/polarisation/elliptical_env.jl index 3f91889d..b256be83 100644 --- a/examples/polarisation/elliptical_env.jl +++ b/examples/polarisation/elliptical_env.jl @@ -1,18 +1,12 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps, Polarisation -import Logging -import FFTW +using Luna import StaticArrays: SMatrix, SVector import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt a = 50e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 400e-9 energy = 30e-6 @@ -22,7 +16,7 @@ nmodes = length(modes) grid = Grid.EnvGrid(10e-2, λ0, (160e-9, 3000e-9), 1e-12) -energyfun = NonlinearRHS.energy_env_modal() +energyfun = Fields.energyfuncs(grid)[1] normfun = NonlinearRHS.norm_modal(grid.ω) function gausspulse(t) @@ -30,8 +24,9 @@ function gausspulse(t) Et = @. sqrt(It) end +field = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) Etlin = gausspulse(grid.t) -cenergy = energyfun(grid.t, Etlin) +cenergy = energyfun(Etlin) Etlin = sqrt(energy)/sqrt(cenergy) .* Etlin Et = [Etlin zero(grid.t)] .+ 0im @@ -51,20 +46,22 @@ densityfun(z) = dens0 responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) # dummy, we don't use these -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :xy; full=false) Eω .= Ew -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW +import PyPlot:pygui, plt + ω = FFTW.fftshift(grid.ω) t = grid.t diff --git a/examples/rectangular/rectangular_modal.jl b/examples/rectangular/rectangular_modal.jl index 7542300a..ce5fe8c8 100644 --- a/examples/rectangular/rectangular_modal.jl +++ b/examples/rectangular/rectangular_modal.jl @@ -1,40 +1,23 @@ -import Luna -import Luna: Grid, Maths, RectModes, PhysData, Nonlinear, Ionisation, NonlinearRHS, RK45, Stats, Output, LinearOps -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -import LinearAlgebra: mul!, ldiv! -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - +using Luna a = 50e-6 b = 10e-6 gas = :Ar pres = 5 +L = 15e-2/10 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 5e-6 -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) -#modes = (RectModes.RectMode(a, b, gas, pres, :Ag, n=1, m=1), -# RectModes.RectMode(a, b, gas, pres, :Ag, n=2, m=1)) modes = collect(RectModes.RectMode(a, b, gas, pres, :Ag, n=n, m=m) for m in 1:3 for n in 1:6) nmodes = length(modes) -grid = Grid.RealGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) - -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -44,19 +27,20 @@ ionrate = Ionisation.ionrate_fun!_ADK(ionpot) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) -in1 = (func=gausspulse, energy=5e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :x; full=true) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) #output = Output.HDF5Output("test_full_modal_rect.h5", 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) +import FFTW + ω = grid.ω t = grid.t diff --git a/examples/rectangular/rectangular_modeAvg.jl b/examples/rectangular/rectangular_modeAvg.jl index 9ae087f5..1b012a43 100644 --- a/examples/rectangular/rectangular_modeAvg.jl +++ b/examples/rectangular/rectangular_modeAvg.jl @@ -1,13 +1,4 @@ -import Luna -import Luna: Grid, Maths, RectModes, PhysData, Modes, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - -import PyPlot:pygui, plt +using Luna a = 50e-6 b = 10e-6 @@ -15,22 +6,16 @@ gas = :Ar pres = 5 L = 15e-2 -energy = 5e-6 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 5e-6 grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) m = RectModes.RectMode(a, b, gas, pres, :Al) aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -45,16 +30,17 @@ linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t @@ -74,9 +60,11 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end +import PyPlot: pygui, plt + pygui(true) plt.figure() plt.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog)) diff --git a/examples/rectangular/rectangular_modeAvg_env.jl b/examples/rectangular/rectangular_modeAvg_env.jl index 6d845039..5c794c4f 100644 --- a/examples/rectangular/rectangular_modeAvg_env.jl +++ b/examples/rectangular/rectangular_modeAvg_env.jl @@ -1,33 +1,20 @@ -import Luna -import Luna: Grid, Maths, RectModes, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import DSP.Unwrap: unwrap - -import PyPlot:pygui, plt +using Luna a = 50e-6 b = 10e-6 gas = :Ar pres = 5 +L = 15e-2 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 energy = 5e-6 -grid = Grid.EnvGrid(15e-2, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 1e-12) m = RectModes.RectMode(a, b, gas, pres, :Al) aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - Et = @. sqrt(It) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -35,23 +22,24 @@ densityfun(z) = dens0 ionpot = PhysData.ionisation_potential(gas) ionrate = Ionisation.ionrate_fun!_ADK(ionpot) -responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) +responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t f = FFTW.fftshift(ω, 1)./2π.*1e-15 @@ -70,9 +58,10 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end +import PyPlot:pygui, plt pygui(true) plt.figure() plt.pcolormesh(f, zout, transpose(FFTW.fftshift(Ilog, 1))) diff --git a/examples/tapers/taper_modal.jl b/examples/tapers/taper_modal.jl index 88e4694c..012d65e1 100644 --- a/examples/tapers/taper_modal.jl +++ b/examples/tapers/taper_modal.jl @@ -1,22 +1,16 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 L = 15e-2 -grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) a0 = a aL = 3a/4 @@ -24,18 +18,12 @@ aL = 3a/4 afun = let a0=a0, aL=aL, L=L afun(z) = a0 + (aL-a0)*z/L end -# coren, dens = Capillary.gradient(gas, L, pres, pres); # dens is unused + modes = ( Capillary.MarcatilliMode(afun, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), Capillary.MarcatilliMode(afun, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false) ) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -47,20 +35,21 @@ responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) linop = LinearOps.make_linop(grid, modes, λ0); -energyfun = NonlinearRHS.energy_modal() +energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) -in1 = (func=gausspulse, energy=1e-6) -inputs = ((1,(in1,)),) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y, full=false) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),length(modes)), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t @@ -72,6 +61,7 @@ It = abs2.(Maths.hilbert(Etout)) Ilog = log10.(Maths.normbymax(abs2.(Eout))) +import PyPlot:pygui, plt pygui(true) for i = 1:length(modes) diff --git a/examples/tapers/taper_modeAvg.jl b/examples/tapers/taper_modeAvg.jl index e48528b8..42de9f11 100644 --- a/examples/tapers/taper_modeAvg.jl +++ b/examples/tapers/taper_modeAvg.jl @@ -1,22 +1,16 @@ -import Luna -import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes -import Logging -import FFTW -import NumericalIntegration: integrate, SimpsonEven -Logging.disable_logging(Logging.BelowMinLevel) - -import PyPlot:pygui, plt +using Luna a = 13e-6 gas = :Ar pres = 5 -τ = 30e-15 +τfwhm = 30e-15 λ0 = 800e-9 +energy = 1e-6 L = 15e-2 -grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) +grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) a0 = a aL = 3a/4 @@ -29,13 +23,7 @@ m = Capillary.MarcatilliMode(afun, gas, pres, loss=false, model=:full); aeff(z) = Modes.Aeff(m, z=z) -energyfun = NonlinearRHS.energy_modal() - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -50,16 +38,17 @@ linop, βfun = LinearOps.make_linop(grid, m, λ0); normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) -statsfun = Stats.collect_stats((Stats.ω0(grid), )) -output = Output.MemoryOutput(0, grid.zmax, 201, (length(grid.ω),), statsfun) +statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) +output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output) +import FFTW + ω = grid.ω t = grid.t @@ -79,10 +68,12 @@ zpeak = argmax(dropdims(maximum(It, dims=1), dims=1)) Et = Maths.hilbert(Etout) energy = zeros(length(zout)) for ii = 1:size(Etout, 2) - energy[ii] = energyfun(t, Etout[:, ii]) + energy[ii] = energyfun(Etout[:, ii]) end +import PyPlot:pygui, plt pygui(true) + plt.figure() plt.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog)) plt.clim(-6, 0) diff --git a/src/Fields.jl b/src/Fields.jl new file mode 100644 index 00000000..fbe90f37 --- /dev/null +++ b/src/Fields.jl @@ -0,0 +1,315 @@ +module Fields +import Luna +import Luna: Grid, Maths, PhysData +import NumericalIntegration: integrate, SimpsonEven +import Random: MersenneTwister, AbstractRNG +import Hankel + +abstract type AbstractField end +abstract type TimeField <: AbstractField end + +""" + PulseField(λ0, energy, ϕ, τ0, Itshape) + +Represents a temporal pulse with shape defined by `Itshape`. + +# Fields +- `λ0::Float64`: the central field wavelength +- `energy::Float64`: the pulse energy +- `ϕ::Float64`: the CEO phase +- `τ0::Float64`: the temproal shift from grid time 0 +- `Itshape`: a callable `f(t)` to get the shape of the intensity/power in the time domain +""" +struct PulseField{iT} <: TimeField + λ0::Float64 + energy::Float64 + ϕ::Float64 + τ0::Float64 + Itshape::iT +end + +""" + GaussField(;λ0, τfwhm, energy, ϕ=0.0, τ0=0.0, m=1) + +Construct a (super)Gaussian shaped pulse with intensity/power FWHM `τfwhm`, +superGaussian parameter `m=1` and other parameters as defined for [`PulseField`](@ref). +""" +function GaussField(;λ0, τfwhm, energy, ϕ=0.0, τ0=0.0, m=1) + PulseField(λ0, energy, ϕ, τ0, t -> Maths.gauss(t, fwhm=τfwhm, power=2*m)) +end + +""" + SechField(;λ0, energy, τw=nothing, τfwhm=nothing, ϕ=0.0, τ0=0.0) + +Construct a Sech^2(t/τw) shaped pulse, specifying either the +natural width `τw`, or the intensity/power FWHM `τfwhm`. +Other parameters are as defined for [`PulseField`](@ref). +""" +function SechField(;λ0, energy, τw=nothing, τfwhm=nothing, ϕ=0.0, τ0=0.0) + if τfwhm != nothing + if τw != nothing + error("only one of `τw` or `τfwhm` can be specified") + else + τw = τfwhm/(2*log(1 + sqrt(2))) + end + elseif τw == nothing + error("one of `τw` or `τfwhm` must be specified") + end + PulseField(λ0, energy, ϕ, τ0, t -> sech(t/τw)^2) +end + +""" + make_Et(p::PulseField, grid) + +Create electric field for `PulseField`, either the field (for `RealGrid`) or the envelope (for `EnvGrid`) +""" +function make_Et(p::PulseField, grid::Grid.RealGrid) + t = grid.t .- p.τ0 + ω0 = PhysData.wlfreq(p.λ0) + @. sqrt(p.Itshape(t))*cos(ω0*t + p.ϕ) +end + +function make_Et(p::PulseField, grid::Grid.EnvGrid) + t = grid.t .- p.τ0 + Δω = PhysData.wlfreq(p.λ0) - grid.ω0 + @. sqrt(p.Itshape(t))*exp(im*(p.ϕ + Δω*t)) +end + +""" + (p::PulseField)(Eω, grid, energy_t, FT) + +Add the field to `Eω` for the provided `grid`, `energy_t` function and Fourier transform `FT` +""" +function (p::PulseField)(Eω, grid, energy_t, FT) + Et = make_Et(p, grid) + Eω .+= FT * (sqrt(p.energy)/sqrt(energy_t(Et)) .* Et) +end + +""" + ShotNoise(seed=nothing) + +Creates one photon per mode quantum noise (shot noise) to add to an input field. +The random `seed` can optionally be provided. +""" +struct ShotNoise{rT<:AbstractRNG} <: TimeField + rng::rT +end + +function ShotNoise(seed=nothing) + ShotNoise(MersenneTwister(seed)) +end + +""" + (s::ShotNoise)(Eω, grid) + +Add shotnoise to `Eω` for the provided `grid`. The random `seed` can optionally be provided. +The optional parameters `energy_t` and `FT` are unused and are present for interface +compatibility with [`TimeField`](@ref). +""" +function (s::ShotNoise)(Eω, grid::Grid.RealGrid, energy_t=nothing, FT=nothing) + δω = grid.ω[2] - grid.ω[1] + δt = grid.t[2] - grid.t[1] + amp = @. sqrt(PhysData.ħ*grid.ω/δω) + rFFTamp = sqrt(2π)/2δt*amp + φ = 2π*rand(s.rng, size(Eω)...) + @. Eω += rFFTamp * exp(1im*φ) +end + +function (s::ShotNoise)(Eω, grid::Grid.EnvGrid, energy_t=nothing, FT=nothing) + δω = grid.ω[2] - grid.ω[1] + δt = grid.t[2] - grid.t[1] + amp = zero(grid.ω) + amp[grid.sidx] = @. sqrt(PhysData.ħ*grid.ω[grid.sidx]/δω) + FFTamp = sqrt(2π)/δt*amp + φ = 2π*rand(s.rng, size(Eω)...) + @. Eω += FFTamp * exp(1im*φ) +end + +""" + SpatioTemporalField(λ0, energy, ϕ, τ0, Ishape) + +Represents a spatiotemporal pulse with shape defined by `Ishape`. + +# Fields +- `λ0::Float64`: the central field wavelength +- `energy::Float64`: the pulse energy +- `ϕ::Float64`: the CEO phase +- `τ0::Float64`: the temproal shift from grid time 0 +- `Ishape`: a callable `f(t, xs)` to get the shape of the intensity/power in the time-space domain +""" +struct SpatioTemporalField{iT} <: AbstractField + λ0::Float64 + energy::Float64 + ϕ::Float64 + τ0::Float64 + Ishape::iT + propz::Float64 +end + +"Gaussian temporal-spatial field defined radially" +function GaussGauss(t, r::AbstractVector, fwhm, m, w0) + Maths.gauss.(t, fwhm=fwhm, power=2*m) .* Maths.gauss.(r, w0/2)' +end + +"Gaussian temporal-spatial field defined on x-y grid" +function GaussGauss(t, r::AbstractArray{T,3} where T, fwhm, m, w0) + Maths.gauss.(t, fwhm=fwhm, power=2*m) .* Maths.gauss.(r, w0/2) +end + +""" + GaussGaussField(;λ0, τfwhm, energy, w0, ϕ=0.0, τ0=0.0, m=1) + +Construct a (super)Gaussian shaped pulse with intensity/power FWHM `τfwhm`, +superGaussian parameter `m=1` and Gaussian shaped spatial profile with waist `w0`, +propagation distance from the waist of `propz`, +and other parameters as defined for [`TimeField`](@ref). +""" +function GaussGaussField(;λ0, τfwhm, energy, w0, ϕ=0.0, τ0=0.0, m=1, propz=0.0) + SpatioTemporalField(λ0, energy, ϕ, τ0, + (t, xs) -> GaussGauss(t, xs, τfwhm, m, w0), + propz) +end + +function make_Etr(s::SpatioTemporalField, grid::Grid.RealGrid, spacegrid) + t = grid.t .- s.τ0 + ω0 = PhysData.wlfreq(s.λ0) + sqrt.(s.Ishape(t, spacegrid.r)) .* cos.(ω0.*t .+ s.ϕ) +end + +function make_Etr(s::SpatioTemporalField, grid::Grid.EnvGrid, spacegrid) + t = grid.t .- s.τ0 + Δω = PhysData.wlfreq(s.λ0) - grid.ω0 + sqrt.(s.Ishape(t, spacegrid.r)) .* exp.(im .* (s.ϕ .+ Δω.*t)) +end + +transform(spacegrid::Hankel.QDHT, FT, Etr) = spacegrid * (FT * Etr) + +transform(spacegrid::Grid.FreeGrid, FT, Etr) = FT * Etr + +""" + (s::SpatioTemporalField)(Eωk, grid, spacegrid, energy_t, FT) + +Add the field to `Eωk` for the provided `grid`, `spacegrid` `energy_t` function +and Fourier transform `FT` +""" +function (s::SpatioTemporalField)(Eωk, grid, spacegrid, energy_t, FT) + Etr = make_Etr(s, grid, spacegrid) + Etr .*= sqrt(s.energy)/sqrt(energy_t(Etr)) + lEωk = transform(spacegrid, FT, Etr) + if s.propz != 0.0 + prop!(lEωk, s.propz, grid, spacegrid) + end + Eωk .+= lEωk +end + +# TODO: this for FreeGrid +function prop!(Eωk, z, grid, q::Hankel.QDHT) + kzsq = @. (grid.ω/PhysData.c)^2 - (q.k^2)' + kzsq[kzsq .< 0] .= 0 + kz = sqrt.(kzsq) + @. Eωk *= exp(-1im * z * (kz - grid.ω/PhysData.c)) +end + +"Calculate energy from modal field E(t)" +function energyfuncs(grid::Grid.RealGrid) + function energy_t(Et) + Eta = Maths.hilbert(Et) + return integrate(grid.t, abs2.(Eta), SimpsonEven()) + end + + prefac = 2π/(grid.ω[end]^2) + function energy_ω(Eω) + prefac*integrate(grid.ω, abs2.(Eω), SimpsonEven()) + end + return energy_t, energy_ω +end + +function energyfuncs(grid::Grid.EnvGrid) + function energy_t(Et) + return integrate(grid.t, abs2.(Et), SimpsonEven()) + end + + δω = grid.ω[2] - grid.ω[1] + Δω = length(grid.ω)*δω + prefac = 2π*δω/(Δω^2) + function energy_ω(Eω) + prefac*sum(abs2.(Eω)) + end + return energy_t, energy_ω +end + +function energyfuncs(grid::Grid.RealGrid, q::Hankel.QDHT) + function energy_t(Et) + Eta = Maths.hilbert(Et) + tintg = integrate(grid.t, abs2.(Eta), SimpsonEven()) + return 2π*PhysData.c*PhysData.ε_0/2 * Hankel.integrateR(tintg, q) + end + + prefac = 2π*PhysData.c*PhysData.ε_0/2 * 2π/(grid.ω[end]^2) + function energy_ω(Eω) + ωintg = integrate(grid.ω, abs2.(Eω), SimpsonEven()) + return prefac*Hankel.integrateK(ωintg, q) + end + return energy_t, energy_ω +end + +function energyfuncs(grid::Grid.EnvGrid, q::Hankel.QDHT) + function energy_t(Et) + tintg = integrate(grid.t, abs2.(Et), SimpsonEven()) + return 2π*PhysData.c*PhysData.ε_0/2 * Hankel.integrateR(tintg, q) + end + + δω = grid.ω[2] - grid.ω[1] + Δω = length(grid.ω)*δω + prefac = 2π*PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) + function energy_ω(Eω) + ωintg = dropdims(sum(abs2.(Eω); dims=1), dims=1) + return prefac*Hankel.integrateK(ωintg, q) + end + return energy_t, energy_ω +end + +function energyfuncs(grid::Grid.RealGrid, xygrid::Grid.FreeGrid) + δx = xygrid.x[2] - xygrid.x[1] + δy = xygrid.y[2] - xygrid.y[1] + δt = grid.t[2] - grid.t[1] + prefac_t = PhysData.c*PhysData.ε_0/2 * δx * δy * δt + function energy_t(Et) + Eta = Maths.hilbert(Et) + return prefac_t * sum(abs2.(Eta)) + end + + δω = grid.ω[2] - grid.ω[1] + Δω = grid.ω[end] + δkx = xygrid.kx[2] - xygrid.kx[1] + Δkx = length(xygrid.kx)*δkx + δky = xygrid.ky[2] - xygrid.ky[1] + Δky = length(xygrid.ky)*δky + prefac = PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) * 2π*δkx/(Δkx^2) * 2π*δky/(Δky^2) + energy_ω(Eω) = prefac * sum(abs2.(Eω)) + + return energy_t, energy_ω +end + +function energyfuncs(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid) + δx = xygrid.x[2] - xygrid.x[1] + δy = xygrid.y[2] - xygrid.y[1] + δt = grid.t[2] - grid.t[1] + prefac_t = PhysData.c*PhysData.ε_0/2 * δx * δy * δt + function energy_t(Et) + return prefac_t * sum(abs2.(Et)) + end + + δω = grid.ω[2] - grid.ω[1] + Δω = length(grid.ω)*δω + δkx = xygrid.kx[2] - xygrid.kx[1] + Δkx = length(xygrid.kx)*δkx + δky = xygrid.ky[2] - xygrid.ky[1] + Δky = length(xygrid.ky)*δky + prefac = PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) * 2π*δkx/(Δkx^2) * 2π*δky/(Δky^2) + energy_ω(Eω) = prefac * sum(abs2.(Eω)) + + return energy_t, energy_ω +end + +end diff --git a/src/Luna.jl b/src/Luna.jl index ef8fe402..e1cafe10 100644 --- a/src/Luna.jl +++ b/src/Luna.jl @@ -3,7 +3,7 @@ import FFTW import Hankel import Logging import LinearAlgebra: mul!, ldiv! -import Random: MersenneTwister +Logging.disable_logging(Logging.BelowMinLevel) """ HDF5LOCK @@ -78,50 +78,84 @@ include("Polarisation.jl") include("Tools.jl") include("Plotting.jl") include("Raman.jl") +include("Fields.jl") export Utils, Scans, Output, Maths, PhysData, Grid, RK45, Modes, Capillary, RectModes, Nonlinear, Ionisation, NonlinearRHS, LinearOps, Stats, Polarisation, - Tools, Plotting, Raman, Antiresonant + Tools, Plotting, Raman, Antiresonant, Fields -function setup(grid::Grid.RealGrid, energyfun, densityfun, normfun, responses, inputs, aeff) +# for a tuple of TimeFields we assume all inputs are for mode 1 +function doinput_sm(grid, inputs::Tuple{Vararg{T} where T <: Fields.TimeField}, FT) + out = fill(0.0 + 0.0im, length(grid.ω)) + energy_t = Fields.energyfuncs(grid)[1] + for input! in inputs + input!(out, grid, energy_t, FT) + end + return out +end + +# for a single Fields.TimeField we assume a single input for mode 1 +function doinput_sm(grid, inputs::Fields.TimeField, FT) + doinput_sm(grid, (inputs,), FT) +end + +function setup(grid::Grid.RealGrid, densityfun, normfun, responses, inputs, aeff) Utils.loadFFTwisdom() xo = Array{Float64}(undef, length(grid.to)) FTo = FFTW.plan_rfft(xo, 1, flags=settings["fftw_flag"]) transform = NonlinearRHS.TransModeAvg(grid, FTo, responses, densityfun, normfun, aeff) x = Array{Float64}(undef, length(grid.t)) FT = FFTW.plan_rfft(x, 1, flags=settings["fftw_flag"]) - Eω = make_init(grid, inputs, energyfun, FT) + Eω = doinput_sm(grid, inputs, FT) inv(FT) # create inverse FT plans now, so wisdom is saved inv(FTo) Utils.saveFFTwisdom() Eω, transform, FT end -function setup(grid::Grid.EnvGrid, energyfun, densityfun, normfun, responses, inputs, aeff) +function setup(grid::Grid.EnvGrid, densityfun, normfun, responses, inputs, aeff) Utils.loadFFTwisdom() x = Array{ComplexF64}(undef, length(grid.t)) FT = FFTW.plan_fft(x, 1, flags=settings["fftw_flag"]) xo = Array{ComplexF64}(undef, length(grid.to)) FTo = FFTW.plan_fft(xo, 1, flags=settings["fftw_flag"]) transform = NonlinearRHS.TransModeAvg(grid, FTo, responses, densityfun, normfun, aeff) - Eω = make_init(grid, inputs, energyfun, FT) + Eω = doinput_sm(grid, inputs, FT) inv(FT) # create inverse FT plans now, so wisdom is saved inv(FTo) Utils.saveFFTwisdom() Eω, transform, FT end -# for multimode setup, inputs is a tuple of ((mode_index, inputs), (mode_index, inputs), ..) -function setup(grid::Grid.RealGrid, energyfun, densityfun, normfun, responses, inputs, - modes, components; full=false) +# for a tuple of NamedTuple's with tuple fields we assume all is well +function doinput_mm!(Eω, grid, inputs::Tuple{Vararg{T} where T <: NamedTuple{<:Any, <:Tuple{Vararg{Any}}}}, FT) + energy_t = Fields.energyfuncs(grid)[1] + for input in inputs + out = @view Eω[:,input.mode] + for input! in input.fields + input!(out, grid, energy_t, FT) + end + end +end + +# for a tuple of TimeFields we assume all inputs are for mode 1 +function doinput_mm!(Eω, grid, inputs::Tuple{Vararg{T} where T <: Fields.TimeField}, FT) + doinput_mm!(Eω, grid, ((mode=1, fields=inputs),), FT) +end + +# for a single Fields.TimeField we assume a single input for mode 1 +function doinput_mm!(Eω, grid, inputs::Fields.TimeField, FT) + doinput_mm!(Eω, grid, ((mode=1, fields=(inputs,)),), FT) +end + +function setup(grid::Grid.RealGrid, densityfun, normfun, responses, inputs, + modes::Modes.ModeCollection, components; full=false) ts = Modes.ToSpace(modes, components=components) Utils.loadFFTwisdom() xt = Array{Float64}(undef, length(grid.t)) FTt = FFTW.plan_rfft(xt, 1, flags=settings["fftw_flag"]) Eω = zeros(ComplexF64, length(grid.ω), length(modes)) - for i in 1:length(inputs) - Eω[:,inputs[i][1]] .= make_init(grid, inputs[i][2], energyfun, FTt) - end + doinput_mm!(Eω, grid, inputs, FTt) x = Array{Float64}(undef, length(grid.t), length(modes)) FT = FFTW.plan_rfft(x, 1, flags=settings["fftw_flag"]) xo = Array{Float64}(undef, length(grid.to), ts.npol) @@ -135,17 +169,14 @@ function setup(grid::Grid.RealGrid, energyfun, densityfun, normfun, responses, i Eω, transform, FT end -# for multimode setup, inputs is a tuple of ((mode_index, inputs), (mode_index, inputs), ..) -function setup(grid::Grid.EnvGrid, energyfun, densityfun, normfun, responses, inputs, - modes, components; full=false) +function setup(grid::Grid.EnvGrid, densityfun, normfun, responses, inputs, + modes::Modes.ModeCollection, components; full=false) ts = Modes.ToSpace(modes, components=components) Utils.loadFFTwisdom() xt = Array{ComplexF64}(undef, length(grid.t)) FTt = FFTW.plan_fft(xt, 1, flags=settings["fftw_flag"]) Eω = zeros(ComplexF64, length(grid.ω), length(modes)) - for i in 1:length(inputs) - Eω[:,inputs[i][1]] .= make_init(grid, inputs[i][2], energyfun, FTt) - end + doinput_mm!(Eω, grid, inputs, FTt) x = Array{ComplexF64}(undef, length(grid.t), length(modes)) FT = FFTW.plan_fft(x, 1, flags=settings["fftw_flag"]) xo = Array{ComplexF64}(undef, length(grid.to), ts.npol) @@ -159,16 +190,27 @@ function setup(grid::Grid.EnvGrid, energyfun, densityfun, normfun, responses, in Eω, transform, FT end +function doinputs_fs!(Eωk, grid, spacegrid::Union{Hankel.QDHT,Grid.FreeGrid}, FT, + inputs::Tuple{Vararg{T} where T <: Fields.SpatioTemporalField}) + energy_t = Fields.energyfuncs(grid, spacegrid)[1] + for input! in inputs + input!(Eωk, grid, spacegrid, energy_t, FT) + end +end + +function doinputs_fs!(Eωk, grid, spacegrid::Union{Hankel.QDHT,Grid.FreeGrid}, FT, + inputs::Fields.SpatioTemporalField) + doinputs_fs!(Eωk, grid, spacegrid, FT, (inputs,)) +end + function setup(grid::Grid.RealGrid, q::Hankel.QDHT, - energyfun, densityfun, normfun, responses, inputs) + densityfun, normfun, responses, inputs) Utils.loadFFTwisdom() xt = zeros(Float64, length(grid.t), length(q.r)) FT = FFTW.plan_rfft(xt, 1, flags=settings["fftw_flag"]) Eω = zeros(ComplexF64, length(grid.ω), length(q.k)) - for input in inputs - Eω .+= scaled_input(grid, input, energyfun, FT) - end Eωk = q * Eω + doinputs_fs!(Eωk, grid, q, FT, inputs) xo = Array{Float64}(undef, length(grid.to), length(q.r)) FTo = FFTW.plan_rfft(xo, 1, flags=settings["fftw_flag"]) transform = NonlinearRHS.TransRadial(grid, q, FTo, responses, densityfun, normfun) @@ -179,15 +221,13 @@ function setup(grid::Grid.RealGrid, q::Hankel.QDHT, end function setup(grid::Grid.EnvGrid, q::Hankel.QDHT, - energyfun, densityfun, normfun, responses, inputs) + densityfun, normfun, responses, inputs) Utils.loadFFTwisdom() xt = zeros(ComplexF64, length(grid.t), length(q.r)) FT = FFTW.plan_fft(xt, 1, flags=settings["fftw_flag"]) Eω = zeros(ComplexF64, length(grid.ω), length(q.k)) - for input in inputs - Eω .+= scaled_input(grid, input, energyfun, FT) - end Eωk = q * Eω + doinputs_fs!(Eωk, grid, q, FT, inputs) xo = Array{ComplexF64}(undef, length(grid.to), length(q.r)) FTo = FFTW.plan_fft(xo, 1, flags=settings["fftw_flag"]) transform = NonlinearRHS.TransRadial(grid, q, FTo, responses, densityfun, normfun) @@ -198,16 +238,14 @@ function setup(grid::Grid.EnvGrid, q::Hankel.QDHT, end function setup(grid::Grid.RealGrid, xygrid::Grid.FreeGrid, - energyfun, densityfun, normfun, responses, inputs) + densityfun, normfun, responses, inputs) Utils.loadFFTwisdom() x = xygrid.x y = xygrid.y xr = Array{Float64}(undef, length(grid.t), length(y), length(x)) FT = FFTW.plan_rfft(xr, (1, 2, 3), flags=settings["fftw_flag"]) Eωk = zeros(ComplexF64, length(grid.ω), length(y), length(x)) - for input in inputs - Eωk .+= scaled_input(grid, input, energyfun, FT) - end + doinputs_fs!(Eωk, grid, xygrid, FT, inputs) xo = Array{Float64}(undef, length(grid.to), length(y), length(x)) FTo = FFTW.plan_rfft(xo, (1, 2, 3), flags=settings["fftw_flag"]) transform = NonlinearRHS.TransFree(grid, xygrid, FTo, @@ -219,16 +257,14 @@ function setup(grid::Grid.RealGrid, xygrid::Grid.FreeGrid, end function setup(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid, - energyfun, densityfun, normfun, responses, inputs) + densityfun, normfun, responses, inputs) Utils.loadFFTwisdom() x = xygrid.x y = xygrid.y xr = Array{ComplexF64}(undef, length(grid.t), length(y), length(x)) - FT = FFTW.plan_rfft(xr, (1, 2, 3), flags=settings["fftw_flag"]) + FT = FFTW.plan_fft(xr, (1, 2, 3), flags=settings["fftw_flag"]) Eωk = zeros(ComplexF64, length(grid.ω), length(y), length(x)) - for input in inputs - Eωk .+= scaled_input(grid, input, energyfun, FT) - end + doinputs_fs!(Eωk, grid, xygrid, FT, inputs) xo = Array{ComplexF64}(undef, length(grid.to), length(y), length(x)) FTo = FFTW.plan_fft(xo, (1, 2, 3), flags=settings["fftw_flag"]) transform = NonlinearRHS.TransFree(grid, xygrid, FTo, @@ -239,46 +275,6 @@ function setup(grid::Grid.EnvGrid, xygrid::Grid.FreeGrid, Eωk, transform, FT end -function make_init(grid, inputs, energyfun, FT) - out = fill(0.0 + 0.0im, length(grid.ω)) - for input in inputs - out .+= scaled_input(grid, input, energyfun, FT) - end - return out -end - -function scaled_input(grid, input, energyfun, FT) - Et = fixtype(grid, input.func(grid.t)) - energy = energyfun(grid.t, Et) - Et_sc = sqrt(input.energy)/sqrt(energy) .* Et - return FT * Et_sc -end - -# Make sure that envelope fields are complex to trigger correct dispatch -fixtype(grid::Grid.RealGrid, Et) = Et -fixtype(grid::Grid.EnvGrid, Et) = complex(Et) - -function shotnoise!(Eω, grid::Grid.RealGrid; seed=nothing) - rng = MersenneTwister(seed) - δω = grid.ω[2] - grid.ω[1] - δt = grid.t[2] - grid.t[1] - amp = @. sqrt(PhysData.ħ*grid.ω/δω) - rFFTamp = sqrt(2π)/2δt*amp - φ = 2π*rand(rng, size(Eω)...) - @. Eω += rFFTamp * exp(1im*φ) -end - -function shotnoise!(Eω, grid::Grid.EnvGrid; seed=nothing) - rng = MersenneTwister(seed) - δω = grid.ω[2] - grid.ω[1] - δt = grid.t[2] - grid.t[1] - amp = zero(grid.ω) - amp[grid.sidx] = @. sqrt(PhysData.ħ*grid.ω[grid.sidx]/δω) - FFTamp = sqrt(2π)/δt*amp - φ = 2π*rand(rng, size(Eω)...) - @. Eω += FFTamp * exp(1im*φ) -end - linoptype(l::AbstractArray) = "constant" linoptype(l) = "variable" diff --git a/src/Maths.jl b/src/Maths.jl index 18fb7722..0b8740f3 100644 --- a/src/Maths.jl +++ b/src/Maths.jl @@ -23,7 +23,7 @@ end "Gaussian or hypergaussian function (with std dev σ as input)" function gauss(x, σ; x0 = 0, power = 2) - return @. exp(-1//2 * ((x-x0)/σ)^power) + return exp(-1/2 * ((x-x0)/σ)^power) end "Gaussian or hypergaussian function (with FWHM as input)" diff --git a/src/Modes.jl b/src/Modes.jl index 3a2d9a0a..18e23d19 100644 --- a/src/Modes.jl +++ b/src/Modes.jl @@ -12,6 +12,9 @@ export dimlimits, neff, β, α, losslength, transmission, dB_per_m, dispersion, abstract type AbstractMode end +ModeCollection = Union{Tuple{Vararg{T} where T <: Modes.AbstractMode}, + AbstractArray{T} where T <: Modes.AbstractMode} + # make modes broadcast like a scalar Broadcast.broadcastable(m::AbstractMode) = Ref(m) diff --git a/src/NonlinearRHS.jl b/src/NonlinearRHS.jl index 6d2ebc1c..914ac13c 100644 --- a/src/NonlinearRHS.jl +++ b/src/NonlinearRHS.jl @@ -194,14 +194,14 @@ function TransModal(grid::Grid.EnvGrid, args...; kwargs...) TransModal(ComplexF64, grid, args...; kwargs...) end -@noinline function reset!(t::TransModal, Emω::Array{ComplexF64,2}, z::Float64) +function reset!(t::TransModal, Emω::Array{ComplexF64,2}, z::Float64) t.Emω .= Emω t.ncalls = 0 t.z = z t.dimlimits = Modes.dimlimits(t.ts.ms[1], z=z) end -@noinline function pointcalc!(fval, xs, t::TransModal) +function pointcalc!(fval, xs, t::TransModal) # TODO: parallelize this in Julia 1.3 for i in 1:size(xs, 2) x1 = xs[1, i] @@ -312,34 +312,6 @@ function (t::TransModeAvg)(nl, Eω, z) nl .*= t.grid.ωwin.*t.densityfun(z).*(-im.*t.grid.ω./2)./t.normfun(z) end -"Calculate energy from modal field E(t)" -function energy_modal(grid::Grid.RealGrid) - function energy_t(t, Et) - Eta = Maths.hilbert(Et) - return integrate(grid.t, abs2.(Eta), SimpsonEven()) - end - - prefac = 2π/(grid.ω[end]^2) - function energy_ω(ω, Eω) - prefac*integrate(ω, abs2.(Eω), SimpsonEven()) - end - return energy_t, energy_ω -end - -function energy_modal(grid::Grid.EnvGrid) - function energy_t(t, Et) - return integrate(grid.t, abs2.(Et), SimpsonEven()) - end - - δω = grid.ω[2] - grid.ω[1] - Δω = length(grid.ω)*δω - prefac = 2π*δω/(Δω^2) - function energy_ω(ω, Eω) - prefac*sum(abs2.(Eω)) - end - return energy_t, energy_ω -end - """ TransRadial @@ -466,37 +438,6 @@ function norm_radial(grid, q, nfun) return norm end -function energy_radial(grid::Grid.RealGrid, q) - function energy_t(t, Et) - Eta = Maths.hilbert(Et) - tintg = integrate(t, abs2.(Eta), SimpsonEven()) - return 2π*PhysData.c*PhysData.ε_0/2 * Hankel.integrateR(tintg, q) - end - - prefac = 2π*PhysData.c*PhysData.ε_0/2 * 2π/(grid.ω[end]^2) - function energy_ω(ω, Eω) - ωintg = integrate(ω, abs2.(Eω), SimpsonEven()) - return prefac*Hankel.integrateK(ωintg, q) - end - return energy_t, energy_ω -end - -function energy_radial(grid::Grid.EnvGrid, q) - function energy_t(t, Et) - tintg = integrate(t, abs2.(Et), SimpsonEven()) - return 2π*PhysData.c*PhysData.ε_0/2 * Hankel.integrateR(tintg, q) - end - - δω = grid.ω[2] - grid.ω[1] - Δω = length(grid.ω)*δω - prefac = 2π*PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) - function energy_ω(ω, Eω) - ωintg = dropdims(sum(abs2.(Eω); dims=1), dims=1) - return prefac*Hankel.integrateK(ωintg, q) - end - return energy_t, energy_ω -end - mutable struct TransFree{TT, FTT, nT, rT, gT, xygT, dT, iT} FT::FTT # 3D Fourier transform (space to k-space and time to frequency) normfun::nT # Function which returns normalisation factor @@ -629,47 +570,4 @@ function norm_free(grid, xygrid, nfun) end end -function energy_free(grid::Grid.RealGrid, xygrid) - δx = xygrid.x[2] - xygrid.x[1] - δy = xygrid.y[2] - xygrid.y[1] - δt = grid.t[2] - grid.t[1] - prefac_t = PhysData.c*PhysData.ε_0/2 * δx * δy * δt - function energy_t(t, Et) - Eta = Maths.hilbert(Et) - return prefac_t * sum(abs2.(Eta)) - end - - δω = grid.ω[2] - grid.ω[1] - Δω = grid.ω[end] - δkx = xygrid.kx[2] - xygrid.kx[1] - Δkx = length(xygrid.kx)*δkx - δky = xygrid.ky[2] - xygrid.ky[1] - Δky = length(xygrid.ky)*δky - prefac = PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) * 2π*δkx/(Δkx^2) * 2π*δky/(Δky^2) - energy_ω(ω, Eω) = prefac * sum(abs2.(Eω)) - - return energy_t, energy_ω -end - -function energy_free(grid::Grid.EnvGrid, xygrid) - δx = xygrid.x[2] - xygrid.x[1] - δy = xygrid.y[2] - xygrid.y[1] - δt = grid.t[2] - grid.t[1] - prefac_t = PhysData.c*PhysData.ε_0/2 * δx * δy * δt - function energy_t(t, Et) - return prefac_t * sum(abs2.(Et)) - end - - δω = grid.ω[2] - grid.ω[1] - Δω = length(grid.ω)*δω - δkx = xygrid.kx[2] - xygrid.kx[1] - Δkx = length(xygrid.kx)*δkx - δky = xygrid.ky[2] - xygrid.ky[1] - Δky = length(xygrid.ky)*δky - prefac = PhysData.c*PhysData.ε_0/2 * 2π*δω/(Δω^2) * 2π*δkx/(Δkx^2) * 2π*δky/(Δky^2) - energy_ω(ω, Eω) = prefac * sum(abs2.(Eω)) - - return energy_t, energy_ω -end - end \ No newline at end of file diff --git a/src/Stats.jl b/src/Stats.jl index b7073037..dd9df0ef 100644 --- a/src/Stats.jl +++ b/src/Stats.jl @@ -28,9 +28,9 @@ Create stats function to calculate the total energy. function energy(grid, energyfun_ω) function addstat!(d, Eω, Et, z, dz) if ndims(Eω) > 1 - d["energy"] = [energyfun_ω(grid.ω, Eω[:, i]) for i=1:size(Eω, 2)] + d["energy"] = [energyfun_ω(Eω[:, i]) for i=1:size(Eω, 2)] else - d["energy"] = energyfun_ω(grid.ω, Eω) + d["energy"] = energyfun_ω(Eω) end end return addstat! @@ -63,9 +63,9 @@ function energy_window(grid, energyfun_ω, window; label) key = "energy_$label" function addstat!(d, Eω, Et, z, dz) if ndims(Eω) > 1 - d[key] = [energyfun_ω(grid.ω, Eω[:, i].*window) for i=1:size(Eω, 2)] + d[key] = [energyfun_ω(Eω[:, i].*window) for i=1:size(Eω, 2)] else - d[key] = energyfun_ω(grid.ω, Eω.*window) + d[key] = energyfun_ω(Eω.*window) end end return addstat! diff --git a/test/runtests.jl b/test/runtests.jl index 9752b0a2..08731836 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -117,4 +117,9 @@ end include(joinpath(testdir, "test_antiresonant.jl")) end +@testset "Fields" begin + @info("================= test_fields.jl") + include(joinpath(testdir, "test_fields.jl")) +end + end \ No newline at end of file diff --git a/test/test_fields.jl b/test/test_fields.jl new file mode 100644 index 00000000..78c81d81 --- /dev/null +++ b/test/test_fields.jl @@ -0,0 +1,397 @@ +import Test: @test, @testset +import Luna: Fields, FFTW, Grid, Maths, PhysData + +# note that most of the Fields.jl code is tested in many other modules + +function getceo(t, Et, It, ω0) + Δt = t[argmax(It)] - t[argmax(Et)] + Δt*ω0 +end + +@testset "Wavelength" begin + # real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + λ0 = 320e-9 + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + λ0 = 800e-9 + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + λ0 = 320e-9 + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + λ0 = 320e-9 + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + λ0 = 800e-9 + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) + λ0 = 320e-9 + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4) +end + +@testset "Energy" begin + # real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + @test isapprox(energy_t(Et), energy, rtol=1e-14) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + @test isapprox(energy_t(Et), energy, rtol=1e-14) + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + @test isapprox(energy_t(Et), energy, rtol=1e-14) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + @test isapprox(energy_t(Et), energy, rtol=1e-14) +end + +@testset "Duration" begin + # real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=1e-5) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=1e-5) + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=2e-5) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=3e-5) +end + +@testset "Position" begin + # real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + # non zero + τ0 = -564e-15 + + #real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15) +end + +@testset "CEO" begin + # real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(getceo(grid.t, Et, It, PhysData.wlfreq(λ0)), ϕ, rtol=1e-15, atol=1e-15) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(getceo(grid.t, Et, It, PhysData.wlfreq(λ0)), ϕ, rtol=1e-15, atol=1e-15) + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + ϕ = 0.0 + τ0 = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(getceo(grid.t, real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0)), ϕ, rtol=1e-15, atol=1e-15) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox(getceo(grid.t, real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0)), ϕ, rtol=1e-15, atol=1e-15) + + # non zero + + #real + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + τ0 = 0.0 + grid = Grid.RealGrid(1.0, λ0, (100e-9, 3000e-9), 1e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{Float64}(undef, length(grid.t)) + FT = FFTW.plan_rfft(x, 1) + + # Make CEO exact multiple of one grid point to avoid issues with argmax() in getceo() + δt = grid.t[2] - grid.t[1] + for i = 1:10 + ϕ = i*δt*PhysData.wlfreq(λ0) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(getceo(grid.t, Et, It, PhysData.wlfreq(λ0)), ϕ, rtol=1e-10) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Maths.hilbert(Et)) + @test isapprox(getceo(grid.t, Et, It, PhysData.wlfreq(λ0)), ϕ, rtol=1e-10) + end + + # Envelope + τfwhm = 30e-15 + λ0 = 800e-9 + energy = 1e-6 + τ0 = 0.0 + grid = Grid.EnvGrid(1.0, λ0, (100e-9, 3000e-9), 1e-12) + energy_t = Fields.energyfuncs(grid)[1] + x = Array{ComplexF64}(undef, length(grid.t)) + FT = FFTW.plan_fft(x, 1) + + # Make CEO exact multiple of one grid point to avoid issues with argmax() in getceo() + δt = grid.t[2] - grid.t[1] + + for i = 1:10 + ϕ = i*δt*PhysData.wlfreq(λ0) + + input! = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox( + getceo(grid.t, real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0)), + ϕ, + rtol=1e-10) + + input! = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ, τ0=τ0) + Eω = fill(0.0 + 0.0im, length(grid.ω)) + input!(Eω, grid, energy_t, FT) + Et = FT \ Eω + It = abs2.(Et) + @test isapprox( + getceo(grid.t,real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0)), + ϕ, + rtol=1e-10) + end +end diff --git a/test/test_full_freespace.jl b/test/test_full_freespace.jl index 22d106b4..32e43638 100644 --- a/test/test_full_freespace.jl +++ b/test/test_full_freespace.jl @@ -1,9 +1,8 @@ import Test: @test, @testset @testset "Full 3D propagation" begin -import Luna -import Luna: Grid, Maths, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Hankel -import Logging +using Luna +Luna.set_fftw_mode(:estimate) import FFTW import Luna.PhysData: wlfreq import LinearAlgebra: norm @@ -26,13 +25,7 @@ xygrid = Grid.FreeGrid(R, N) x = xygrid.x y = xygrid.y -energyfun, energyfunω = NonlinearRHS.energy_free(grid, xygrid) - -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) .* Maths.gauss.(xygrid.r, w0/2) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end +energyfun, energyfunω = Fields.energyfuncs(grid, xygrid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 @@ -42,10 +35,9 @@ responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) linop = LinearOps.make_const_linop(grid, xygrid, PhysData.ref_index_fun(gas, pres)) normfun = NonlinearRHS.const_norm_free(grid, xygrid, PhysData.ref_index_fun(gas, pres)) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0) -Eω, transform, FT = Luna.setup(grid, xygrid, energyfun, densityfun, normfun, responses, inputs) +Eω, transform, FT = Luna.setup(grid, xygrid, densityfun, normfun, responses, inputs) output = Output.MemoryOutput(0, grid.zmax, 21) @@ -56,7 +48,7 @@ Eout = output.data["Eω"] # (ω, ky, kx, z) ω0idx = argmin(abs.(grid.ω .- wlfreq(λ0))) λ0 = 2π*PhysData.c/grid.ω[ω0idx] w1 = w0*sqrt(1+(L*λ0/(π*w0^2))^2) -Iω0_analytic = Maths.gauss(xygrid.x, w1/2) # analytical solution (in paraxial approx) +Iω0_analytic = Maths.gauss.(xygrid.x, w1/2) # analytical solution (in paraxial approx) Eω0yx = FFTW.ifft(Eout[ω0idx, :, :, end], (1, 2)) Iω0yx = abs2.(Eω0yx) diff --git a/test/test_gradient.jl b/test/test_gradient.jl index 42927f5a..0cede5be 100644 --- a/test/test_gradient.jl +++ b/test/test_gradient.jl @@ -12,14 +12,8 @@ pres = 5 L = 5e-2 # Common setup -grid = Grid.RealGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) +grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) # Constant @@ -27,11 +21,11 @@ dens0 = PhysData.density(gas, pres) dens(z) = dens0 m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) Eω, transform, FT = Luna.setup( - grid, energyfun, dens, normfun, responses, inputs, aeff) + grid, dens, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -42,11 +36,11 @@ Luna.run(Eω, grid, linop, transform, FT, output_const, status_period=10) coren, densityfun = Capillary.gradient(gas, L, pres, pres) m = Capillary.MarcatilliMode(a, coren, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) linop, βfun = LinearOps.make_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -57,11 +51,11 @@ Luna.run(Eω, grid, linop, transform, FT, output_grad, status_period=10) coren, densityfun = Capillary.gradient(gas, [0,L], [pres, pres]); m = Capillary.MarcatilliMode(a, coren, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) linop, βfun = LinearOps.make_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -81,14 +75,8 @@ pres = 5 L = 5e-2 # Common setup -grid = Grid.EnvGrid(L, 800e-9, (160e-9, 3000e-9), 1e-12) -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It) -end -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) +grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 1e-12) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),) # Constant @@ -96,10 +84,10 @@ dens0 = PhysData.density(gas, pres) dens(z) = dens0 m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid); +energyfun, energyfunω = Fields.energyfuncs(grid); linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0); normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -Eω, transform, FT = Luna.setup(grid, energyfun, dens, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, dens, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -110,11 +98,11 @@ Luna.run(Eω, grid, linop, transform, FT, output_const, status_period=10) coren, densityfun = Capillary.gradient(gas, L, pres, pres) m = Capillary.MarcatilliMode(a, coren, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) linop, βfun = LinearOps.make_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -125,11 +113,11 @@ Luna.run(Eω, grid, linop, transform, FT, output_grad, status_period=10) coren, densityfun = Capillary.gradient(gas, [0,L], [pres, pres]); m = Capillary.MarcatilliMode(a, coren, loss=false) aeff(z) = Modes.Aeff(m, z=z) -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) linop, βfun = LinearOps.make_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) diff --git a/test/test_maths.jl b/test/test_maths.jl index 6263066a..3882bb1c 100644 --- a/test/test_maths.jl +++ b/test/test_maths.jl @@ -23,7 +23,7 @@ end @testset "Moments" begin x = collect(range(-10, stop=10, length=513)) - y = Maths.gauss(x, 1, x0=1) + y = Maths.gauss.(x, 1, x0=1) @test Maths.moment(x, y) ≈ 1 @test Maths.moment(x, y, 2) ≈ 2 @test Maths.rms_width(x, y) ≈ 1 @@ -32,7 +32,7 @@ end σ = [0.1, 0.5, 1.0, 1.5, 1.5] y = zeros(length(x), length(x0)) for ii = 1:length(x0) - y[:, ii] = Maths.gauss(x, σ[ii], x0=x0[ii]) + y[:, ii] = Maths.gauss.(x, σ[ii], x0=x0[ii]) end @test_throws DomainError Maths.moment(x, y, dim=2) @test all(isapprox.(transpose(Maths.moment(x, y, dim=1)), x0, atol=1e-5)) @@ -46,7 +46,7 @@ end @testset "Fourier" begin t = collect(range(-10, stop=10, length=513)) - Et = Maths.gauss(t, fwhm=4).*cos.(4*t) + Et = Maths.gauss.(t, fwhm=4).*cos.(4*t) EtA = Maths.hilbert(Et) @test maximum(abs.(EtA)) ≈ 1 @test all(isapprox.(real(EtA), Et, atol=1e-9)) @@ -61,12 +61,12 @@ end @test all(out .≈ EtA) t = collect(range(-10, stop=10, length=512)) - Et = Maths.gauss(t, fwhm=4).*cos.(4*t) + Et = Maths.gauss.(t, fwhm=4).*cos.(4*t) to, Eto = Maths.oversample(t, Et, factor=4) @test 4*size(Et)[1] == size(Eto)[1] @test all(isapprox.(Eto[1:4:end], Et, rtol=1e-6)) - Etc = Maths.gauss(t, fwhm=4).*exp.(1im*4*t) + Etc = Maths.gauss.(t, fwhm=4).*exp.(1im*4*t) to, Etco = Maths.oversample(t, Etc, factor=4) @test 4*size(Etc)[1] == size(Etco)[1] @test all(isapprox.(Etco[1:4:end], Etc, rtol=1e-6)) @@ -120,7 +120,7 @@ end fw = 0.2 x = collect(range(-1, stop=1, length=2048)) δx = x[2] - x[1] - y = Maths.gauss(x, fwhm=fw) + y = Maths.gauss.(x, fwhm=fw) @test isapprox(Maths.fwhm(x, y), fw, rtol=1e-4) @test isapprox(Maths.fwhm(x, y, method=:spline), fw, rtol=1e-4) @test Maths.fwhm(x, y, method=:spline) == Maths.fwhm(x, y, method=:spline, minmax=:max) @@ -128,14 +128,14 @@ end fw = 0.1 sep = 0.5 - y = Maths.gauss(x, fwhm=fw, x0=-sep/2) + Maths.gauss(x, fwhm=fw, x0=sep/2) + y = Maths.gauss.(x, fwhm=fw, x0=-sep/2) + Maths.gauss.(x, fwhm=fw, x0=sep/2) @test isapprox(Maths.fwhm(x, y, minmax=:min), fw, rtol=1e-4) @test isapprox(Maths.fwhm(x, y, method=:spline, minmax=:min), fw, rtol=1e-4) @test isapprox(Maths.fwhm(x, y, minmax=:max), fw+sep, rtol=1e-4) @test isapprox(Maths.fwhm(x, y, method=:spline, minmax=:max), fw+sep, rtol=1e-4) hw = 1 - f(x) = Maths.gauss(x, fwhm=2*hw) + f(x) = Maths.gauss.(x, fwhm=2*hw) @test Maths.hwhm(f, 0) ≈ hw @test Maths.hwhm(f, 0; direction=:bwd) ≈ hw @@ -144,11 +144,11 @@ end y = zero(x) seed!(123) for i = 1:10000 - y .+= rand()*Maths.gauss(x, x0=rand(), fwhm=rand()/1000) + y .+= rand()*Maths.gauss.(x, x0=rand(), fwhm=rand()/1000) end m = 1.1*maximum(y) - y .+= m*Maths.gauss(x, x0=-0.05, fwhm=0.05) - y .+= m*Maths.gauss(x, x0=1.05, fwhm=0.05) + y .+= m*Maths.gauss.(x, x0=-0.05, fwhm=0.05) + y .+= m*Maths.gauss.(x, x0=1.05, fwhm=0.05) @test isapprox(Maths.fwhm(x, y, minmax=:max), 0.05 + 1.1, rtol=1e-4) @test isapprox(Maths.fwhm(x, y, minmax=:max, method=:spline), 0.05 + 1.1, rtol=1e-4) end diff --git a/test/test_modes.jl b/test/test_modes.jl index 952e3948..d55e7d25 100644 --- a/test/test_modes.jl +++ b/test/test_modes.jl @@ -64,7 +64,7 @@ r = collect(range(0, 4a, length=2^16)) m = Capillary.MarcatilliMode(a, :He, 1.0, model=:reduced, m=1) for i in eachindex(fac) w0 = fac[i]*a - Er = Maths.gauss(r, w0/sqrt(2)) + Er = Maths.gauss.(r, w0/sqrt(2)) ηn[i] = abs2.(Modes.overlap(m, r, Er, dim=1)[1]) end @test 0.63 < fac[argmax(ηn)] < 0.65 @@ -73,7 +73,7 @@ end total energy =# fac = 0.45 w0 = fac*a -Er = Maths.gauss(r, w0/sqrt(2)) +Er = Maths.gauss.(r, w0/sqrt(2)) s = 0 for mi = 1:10 m = Capillary.MarcatilliMode(a, :He, 1.0, model=:reduced, m=mi) @@ -93,7 +93,7 @@ a = 100e-6 # capillary radius q = Hankel.QDHT(2a, 512) grid = Grid.RealGrid(1, 800e-9, (200e-9, 2000e-9), 0.5e-12) # First pulse -It1 = Maths.gauss(grid.t, fwhm=30e-15) +It1 = Maths.gauss.(grid.t, fwhm=30e-15) Et1 = @. sqrt(It1)*cos(2π*PhysData.c/800e-9*grid.t) Eω1 = FFTW.rfft(Et1) # Spatial profile of the first pulse @@ -102,7 +102,7 @@ Er1 = besselj.(0, unm*q.r/a)' Er1[q.r .> a] .= 0 Etr1 = Et1 .* Er1 # create spatio-temporal pulse profile # Second pulse -It2 = 4*Maths.gauss(grid.t, fwhm=15e-15) +It2 = 4*Maths.gauss.(grid.t, fwhm=15e-15) Et2 = @. sqrt(It2)*cos(2π*PhysData.c/400e-9*grid.t) Eω2 = FFTW.rfft(Et2) # Spatial profile of the second pulse @@ -115,10 +115,10 @@ Etr2 = Et2 .* Er2 # create spatio-temporal pulse profile Etr = Etr1 .+ Etr2 Eωr = FFTW.rfft(Etr, 1) -ert, ekω = NonlinearRHS.energy_radial(grid, q) -energy1 = ert(grid.t, Etr1) -energy2 = ert(grid.t, Etr2) -@test ert(grid.t, Etr) ≈ energy1 + energy2 +ert, ekω = Fields.energyfuncs(grid, q) +energy1 = ert(Etr1) +energy2 = ert(Etr2) +@test ert(Etr) ≈ energy1 + energy2 m1 = Capillary.MarcatilliMode(a, :He, 1.0, model=:reduced, m=1) m2 = Capillary.MarcatilliMode(a, :He, 1.0, model=:reduced, m=2) @@ -128,10 +128,10 @@ Eωm2 = Modes.overlap(m2, q.r, Eωr; dim=2, norm=false) Etm1 = FFTW.irfft(Eωm1[:, 1], length(grid.t)) Etm2 = FFTW.irfft(Eωm2[:, 1], length(grid.t)) -et, eω = NonlinearRHS.energy_modal(grid) +et, eω = Fields.energyfuncs(grid) # check that the non-normalised overlap integral preserves the total energy -@test isapprox(et(grid.t, Etm1), energy1, rtol=1e-3) -@test isapprox(et(grid.t, Etm2), energy2, rtol=1e-3) +@test isapprox(et(Etm1), energy1, rtol=1e-3) +@test isapprox(et(Etm2), energy2, rtol=1e-3) # check that the non-normalised overlap integral preserves the spectral shape En1 = Eω1/norm(Eω1) diff --git a/test/test_multimode.jl b/test/test_multimode.jl index d565ea0e..defcf33b 100644 --- a/test/test_multimode.jl +++ b/test/test_multimode.jl @@ -2,9 +2,7 @@ import Test: @test, @testset, @test_throws @testset "Radial" begin # mode average and radial integral for single mode and only Kerr should be identical - - import Luna - import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes + using Luna import LinearAlgebra: norm a = 13e-6 gas = :Ar @@ -14,21 +12,16 @@ import Test: @test, @testset, @test_throws grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end + energyfun, energyfunω = Fields.energyfuncs(grid) + dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) - in1 = (func=gausspulse, energy=1e-6) - inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -38,10 +31,10 @@ import Test: @test, @testset, @test_throws modes = ( Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), ) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) + energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) - inputs = ((1,(in1,)),) - Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) + Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=false) outputr = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) @@ -54,9 +47,7 @@ end @testset "Full" begin # mode average and full integral for single mode and only Kerr should be identical - - import Luna - import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps, Modes + using Luna import LinearAlgebra: norm a = 13e-6 gas = :Ar @@ -66,21 +57,16 @@ end grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end + energyfun, energyfunω = Fields.energyfuncs(grid) + dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) - in1 = (func=gausspulse, energy=1e-6) - inputs = (in1, ) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) @@ -90,10 +76,10 @@ end modes = ( Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), ) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) + energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) - inputs = ((1,(in1,)),) - Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) + Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=true) outputf = Output.MemoryOutput(0, grid.zmax, 201, statsfun) linop = LinearOps.make_const_linop(grid, modes, λ0) @@ -103,3 +89,35 @@ end Iωf = abs2.(dropdims(outputf.data["Eω"], dims=2)) @test norm(Iω - Iωf)/norm(Iω) < 0.003 end + +@testset "FieldInputs" begin + using Luna + import LinearAlgebra: norm + a = 13e-6 + gas = :Ar + pres = 5 + τ = 30e-15 + λ0 = 800e-9 + grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12) + modes = ( + Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), + Capillary.MarcatilliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false) + ) + aeff(z) = Modes.Aeff(m, z=z) + energyfun, energyfunω = Fields.energyfuncs(grid) + dens0 = PhysData.density(gas, pres) + densityfun(z) = dens0 + responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) + normfun = NonlinearRHS.norm_modal(grid.ω) + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) + Eω_single, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, + modes, :y; full=false) + inputs = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6),) + Eω_tuple, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, + modes, :y; full=false) + inputs = ((mode=1, fields=(Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6),)),) + Eω_tuple_of_namedtuples, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, + modes, :y; full=false) + @test Eω_single ≈ Eω_tuple + @test Eω_single ≈ Eω_tuple_of_namedtuples +end diff --git a/test/test_output.jl b/test/test_output.jl index 72b9211b..a16b3ca7 100644 --- a/test/test_output.jl +++ b/test/test_output.jl @@ -111,21 +111,16 @@ fpath_comp = joinpath(homedir(), ".luna", "output_test", "test_comp.h5") grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12) m = Capillary.MarcatilliMode(a, gas, pres, loss=false) aeff(z) = Modes.Aeff(m, z=z) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) + energyfun, energyfunω = Fields.energyfuncs(grid) dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end - in1 = (func=gausspulse, energy=1e-6) - inputs = (in1, ) + + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) Eω, transform, FT = Luna.setup( - grid, energyfun, densityfun, normfun, responses, inputs, aeff) + grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), Stats.energy(grid, energyfunω)) diff --git a/test/test_polarisation_env.jl b/test/test_polarisation_env.jl index a462fdea..3a3de9ae 100644 --- a/test/test_polarisation_env.jl +++ b/test/test_polarisation_env.jl @@ -5,31 +5,26 @@ import Luna: Output # compare radial single pol, to radial linear pol at 45 degrees, # in capillary (non birefringent) these should be identical - import Luna - import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps + using Luna import LinearAlgebra: norm a = 13e-6 gas = :Ar pres = 5 τ = 30e-15 λ0 = 800e-9 + energy=1e-6 grid = Grid.EnvGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12) - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - Et = @. sqrt(It) - end dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) + energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) modes = ( Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), ) - in1 = (func=gausspulse, energy=1e-6) - inputs = ((1,(in1,)),) - Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy) + Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=true) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -44,10 +39,10 @@ import Luna: Output Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=π/2, loss=false) ) - in1 = (func=gausspulse, energy=1e-6/2.0) + inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),) # same field in each mode - inputs = ((1, (in1,)), (2, (in1,))) - Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, + inputs = ((mode=1, fields=inf), (mode=2, fields=inf)) + Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :xy; full=true) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), diff --git a/test/test_polarisation_field.jl b/test/test_polarisation_field.jl index 0b2c5397..4199e2b8 100644 --- a/test/test_polarisation_field.jl +++ b/test/test_polarisation_field.jl @@ -5,31 +5,26 @@ import Luna: Output # compare radial single pol, to radial linear pol at 45 degrees, # in capillary (non birefringent) these should be identical - import Luna - import Luna: Grid, Maths, Capillary, PhysData, Nonlinear, Ionisation, NonlinearRHS, Output, Stats, LinearOps + using Luna import LinearAlgebra: norm a = 13e-6 gas = :Ar pres = 5 τ = 30e-15 λ0 = 800e-9 + energy = 1e-6 grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12) - function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - end + dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) - energyfun, energyfunω = NonlinearRHS.energy_modal(grid) + energyfun, energyfunω = Fields.energyfuncs(grid) normfun = NonlinearRHS.norm_modal(grid.ω) modes = ( Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), ) - in1 = (func=gausspulse, energy=1e-6) - inputs = ((1,(in1,)),) - Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, + inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy) + Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y; full=true) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), @@ -44,10 +39,10 @@ import Luna: Output Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false), Capillary.MarcatilliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=π/2, loss=false) ) - in1 = (func=gausspulse, energy=1e-6/2.0) + inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),) # same field in each mode - inputs = ((1, (in1,)), (2, (in1,))) - Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, + inputs = ((mode=1, fields=inf), (mode=2, fields=inf)) + Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :xy; full=true) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid), diff --git a/test/test_radial.jl b/test/test_radial.jl index b057c828..aa1fbe48 100644 --- a/test/test_radial.jl +++ b/test/test_radial.jl @@ -21,7 +21,7 @@ N = 128 grid = Grid.RealGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12) q = Hankel.QDHT(R, N, dim=2) -energyfun, energyfun_ω = NonlinearRHS.energy_radial(grid, q) +energyfun, energyfun_ω = Fields.energyfuncs(grid, q) function prop(E, z) Eω = FFTW.rfft(E, 1) @@ -35,13 +35,6 @@ function prop(E, z) return E end -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) .* Maths.gauss(q.r, w0/2)' - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) - return prop(Et, -0.3) -end - dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 ionpot = PhysData.ionisation_potential(gas) @@ -50,9 +43,10 @@ responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)), Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)) linop = LinearOps.make_const_linop(grid, q, PhysData.ref_index_fun(gas, pres)) normfun = NonlinearRHS.const_norm_radial(grid, q, PhysData.ref_index_fun(gas, pres)) -in1 = (func=gausspulse, energy=energy) -inputs = (in1, ) -Eω, transform, FT = Luna.setup(grid, q, energyfun, densityfun, normfun, responses, inputs) + +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0, propz=-0.3) + +Eω, transform, FT = Luna.setup(grid, q, densityfun, normfun, responses, inputs) output = Output.MemoryOutput(0, grid.zmax, 201) Luna.run(Eω, grid, linop, transform, FT, output) @@ -64,7 +58,7 @@ Iωr = abs2.(Erout) Iλ0 = Iωr[ω0idx, :, :] λ0 = 2π*PhysData.c/grid.ω[ω0idx] w1 = w0*sqrt(1+(L/2*λ0/(π*w0^2))^2) -Iλ0_analytic = Maths.gauss(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx) +Iλ0_analytic = Maths.gauss.(q.r, w1/2)*(w0/w1)^2 # analytical solution (in paraxial approx) Ir = Maths.normbymax(Iλ0[:, end]) Ira = Maths.normbymax(Iλ0_analytic) @test maximum(abs.(Ir .- Ira)/norm(Ir)) < 1.5e-4 diff --git a/test/test_tapers.jl b/test/test_tapers.jl index 57c5501a..a2119dd3 100644 --- a/test/test_tapers.jl +++ b/test/test_tapers.jl @@ -19,14 +19,9 @@ afun = let a0=a0, aL=aL, L=L afun(z) = a0 + (aL-a0)*z/L end -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) ## mode-average @@ -35,9 +30,8 @@ m = Capillary.MarcatilliMode(afun, gas, pres, loss=false, model=:full); aeff(z) = Modes.Aeff(m, z=z) linop, βfun = LinearOps.make_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=600e-9) -inputs = (in1, ) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=600e-9) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output, status_period=10) @@ -51,9 +45,8 @@ modes = ( ) linop = LinearOps.make_linop(grid, modes, λ0) normfun = NonlinearRHS.norm_modal(grid.ω) -in1 = (func=gausspulse, energy=600e-9) -inputs = ((1,(in1,)),) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, +inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=600e-9) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, modes, :y, full=false) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) @@ -78,14 +71,9 @@ afun = let a=a z -> a end -function gausspulse(t) - It = Maths.gauss(t, fwhm=τ) - ω0 = 2π*PhysData.c/λ0 - Et = @. sqrt(It)*cos(ω0*t) -end dens0 = PhysData.density(gas, pres) densityfun(z) = dens0 -energyfun, energyfunω = NonlinearRHS.energy_modal(grid) +energyfun, energyfunω = Fields.energyfuncs(grid) responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),) ## mode-average @@ -94,9 +82,8 @@ m = Capillary.MarcatilliMode(afun, gas, pres, loss=false, model=:full); aeff(z) = Modes.Aeff(m, z=z) linop, βfun = LinearOps.make_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output, status_period=10) @@ -107,11 +94,10 @@ end Iωavg_c = let m = Capillary.MarcatilliMode(a, gas, pres, loss=false, model=:full); aeff(z) = Modes.Aeff(m, z=z) -in1 = (func=gausspulse, energy=1e-6) -inputs = (in1, ) +inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6) linop, βfun, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0) normfun = NonlinearRHS.norm_mode_average(grid.ω, βfun, aeff) -Eω, transform, FT = Luna.setup(grid, energyfun, densityfun, normfun, responses, inputs, aeff) +Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs, aeff) statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid)) output = Output.MemoryOutput(0, grid.zmax, 201, statsfun) Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)