From f87276b356884aafbb048b6f00cee5d2df80c41e Mon Sep 17 00:00:00 2001 From: John Travers Date: Sun, 5 Sep 2021 23:55:38 +0100 Subject: [PATCH 01/18] Add plasma for envelopes to interface --- src/Interface.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index d26da035..2cbf5935 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -503,8 +503,20 @@ function makeplasma!(out, grid, gas, plasma::Symbol, pol) push!(out, Nonlinear.PlasmaCumtrapz(grid.to, Et, ionrate, ionpot)) end +function makeplasma!(out, grid::Grid.EnvGrid, gas, plasma::Symbol, pol) + ionpot = PhysData.ionisation_potential(gas) + if plasma == :ADK + error("ADK not supported for envelopes at present") + elseif plasma == :PPT + ionrate = Ionisation.ionrate_fun!_PPTcached(gas, grid.referenceλ, rcycle=false) + else + throw(DomainError(plasma, "Unknown ionisation rate $plasma.")) + end + Et = pol ? Array{Complex{Float64}}(undef, length(grid.to), 2) : Vector{Complex{Float64}}(undef, length(grid.to)) + push!(out, Nonlinear.PlasmaCumtrapz(grid.to, Et, ionrate, ionpot)) +end + function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol) - plasma && error("Plasma response for envelope fields has not been implemented yet.") isnothing(thg) && (thg = false) out = Any[] if kerr @@ -516,6 +528,7 @@ function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol) push!(out, Nonlinear.Kerr_env(PhysData.γ3_gas(gas))) end end + makeplasma!(out, grid, gas, plasma, pol) raman && push!(out, Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(gas))) Tuple(out) end From 2c3cb2eb936d6cc137ec4404dbdd9281f701447d Mon Sep 17 00:00:00 2001 From: John Travers Date: Sun, 5 Sep 2021 23:55:49 +0100 Subject: [PATCH 02/18] Soliton blue shift example --- examples/simple_interface/PlasmaSSFBS.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 examples/simple_interface/PlasmaSSFBS.jl diff --git a/examples/simple_interface/PlasmaSSFBS.jl b/examples/simple_interface/PlasmaSSFBS.jl new file mode 100644 index 00000000..c76a8985 --- /dev/null +++ b/examples/simple_interface/PlasmaSSFBS.jl @@ -0,0 +1,19 @@ +using Luna + +radius = 9e-6 # HCF core radius +flength = 25e-2 # HCF length + +gas = :Ar +pressure = 5.0 # gas pressure in bar + +λ0 = 1500e-9 # central wavelength of the pump pulse +τfwhm = 30e-15 # FWHM duration of the pump pulse +energy = 1.7e-6 # energy in the pump pulse + +pssfbs = prop_capillary(radius, flength, gas, pressure; λ0, τfwhm, energy, + trange=1.6e-12, λlims=(300e-9, 1.8e-6), loss=false, envelope=true, plasma=true) + +Plotting.prop_2D(pssfbs, :λ; modes=:sum, trange=(-700e-15, 100e-15), λrange=(300e-9, 1800e-9), + dBmin=-30) +Plotting.time_1D(pssfbs, range(0,flength,length=5), modes=:sum, trange=(-700e-15, 100e-15)) +Plotting.spec_1D(pssfbs, range(0,flength,length=5), modes=:sum, λrange=(300e-9, 1800e-9)) From f34daafae80e5883980fb6c5a417c8c114ccd58a Mon Sep 17 00:00:00 2001 From: John Travers Date: Sun, 5 Sep 2021 23:56:05 +0100 Subject: [PATCH 03/18] Roughly fix stats --- src/Stats.jl | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/src/Stats.jl b/src/Stats.jl index fca8935d..fd4b674d 100644 --- a/src/Stats.jl +++ b/src/Stats.jl @@ -216,7 +216,7 @@ If oversampling > 1, the field is oversampled before the calculation Oversampling can lead to a significant performance hit """ function electrondensity(grid::Grid.RealGrid, ionrate!, dfun, aeff; oversampling=1) - to, Eto = Maths.oversample(grid.t, complex(grid.t), factor=oversampling) + to, Eto = Maths.oversample(grid.t, grid.t, factor=oversampling) δt = to[2] - to[1] # ionfrac! stores the time-dependent ionisation fraction in out and returns the max # ionisation rate @@ -238,6 +238,29 @@ function electrondensity(grid::Grid.RealGrid, ionrate!, dfun, aeff; oversampling end end +function electrondensity(grid::Grid.EnvGrid, ionrate!, dfun, aeff; oversampling=1) + to, Eto = Maths.oversample(grid.t, complex(grid.t), factor=oversampling) + δt = to[2] - to[1] + # ionfrac! stores the time-dependent ionisation fraction in out and returns the max + # ionisation rate + function ionfrac!(out, Et) + ionrate!(out, Et) + ratemax = maximum(out) + Maths.cumtrapz!(out, δt) # in-place cumulative integration + @. out = 1 - exp(-out) + return ratemax + end + frac = similar(to) + function addstat!(d, Eω, Et, z, dz) + # note: oversampling returns its arguments without any work done if factor==1 + to, Eto = Maths.oversample(grid.t, Et, factor=oversampling) + @. Eto /= sqrt(ε_0*c*aeff(z)/2) + ratemax = ionfrac!(frac, abs.(Eto)) + d["electrondensity"] = frac[end]*dfun(z) + d["peak_ionisation_rate"] = ratemax + end +end + """ electrondensity(grid, ionrate, dfun, modes; oversampling=1) From 30351663e781dd70ba3d8200e9e09db3426ad4ab Mon Sep 17 00:00:00 2001 From: John Travers Date: Sun, 5 Sep 2021 23:56:16 +0100 Subject: [PATCH 04/18] Add envelope plasma --- src/Nonlinear.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index f5eb83bc..922c812b 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -140,6 +140,17 @@ function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E) Maths.cumtrapz!(out, Plas.J, Plas.δt) end +"The plasma response for a scalar electric field" +function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E::Vector{Complex{Float64}}) + Plas.ratefunc(Plas.rate, E) + Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) + @. Plas.fraction = 1-exp(-Plas.fraction) + @. Plas.phase = Plas.fraction * e_ratio * E + ω = 2π .* FFTW.fftfreq(length(E),1/Plas.δt) .+ 1.2557677115392355e15 + Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E + out .= FFTW.ifft(-FFTW.fft(Plas.phase) ./ ω.^2 .- 1im .* FFTW.fft(Jabs) ./ ω) +end + """ The plasma response for a vector electric field. From 1dea0ae4ffba343df3297f4f8e523b1c2493c1f9 Mon Sep 17 00:00:00 2001 From: John Travers Date: Mon, 6 Sep 2021 11:56:31 +0100 Subject: [PATCH 05/18] Add PlasmaFourier --- src/Nonlinear.jl | 46 +++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 922c812b..5cb206d7 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -98,9 +98,11 @@ function Kerr_env_thg(γ3, ω0, t) end end +abstract type PlasmaPolar end + "Response type for cumtrapz-based plasma polarisation, adapted from: M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec, Physical Review Letters 83, 2930 (1999)." -struct PlasmaCumtrapz{R, EType, tType} +struct PlasmaCumtrapz{R, EType, tType} <: PlasmaPolar ratefunc::R # the ionization rate function ionpot::Float64 # the ionization potential (for calculation of ionization loss) rate::tType # buffer to hold the rate @@ -146,9 +148,9 @@ function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E::Vector{Complex{Float64}}) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E - ω = 2π .* FFTW.fftfreq(length(E),1/Plas.δt) .+ 1.2557677115392355e15 - Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E - out .= FFTW.ifft(-FFTW.fft(Plas.phase) ./ ω.^2 .- 1im .* FFTW.fft(Jabs) ./ ω) + Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) + Plas.J .+= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E + Maths.cumtrapz!(out, Plas.J, Plas.δt) end """ @@ -179,8 +181,42 @@ function PlasmaVector!(out, Plas::PlasmaCumtrapz, E) Maths.cumtrapz!(out, Plas.J, Plas.δt) end +struct PlasmaFourier{R, EType, tType, FTType} <: PlasmaPolar + ratefunc::R # the ionization rate function + ionpot::Float64 # the ionization potential (for calculation of ionization loss) + rate::tType # buffer to hold the rate + fraction::tType # buffer to hold the ionization fraction + phase::EType # buffer to hold the plasma induced (mostly) phase modulation + J::EType # buffer to hold the plasma current + ω::tType # buffer to hold the frequency grid + FT::FTType # Fourier transform to use for integrations + δt::Float64 +end + +function PlasmaFourier(ω, E, ratefunc, ionpot, δt) + rate = similar(ω) + fraction = similar(ω) + phase = similar(E) + J = similar(E) + Utils.loadFFTwisdom() + FT = FFTW.plan_fft(E, 1, flags=Luna.settings["fftw_flag"]) + inv(FT) + Utils.saveFFTwisdom() + PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, ω, FT, δt) +end + +"The plasma response for a scalar electric field" +function PlasmaScalar!(out, Plas::PlasmaFourier, E::Vector{Complex{Float64}}) + Plas.ratefunc(Plas.rate, E) + Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) + @. Plas.fraction = 1-exp(-Plas.fraction) + @. Plas.phase = Plas.fraction * e_ratio * E + Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E # (E + conj.(E)) + out .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) +end + "Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`." -function (Plas::PlasmaCumtrapz)(out, Et) +function (Plas::PlasmaPolar)(out, Et) if ndims(Et) > 1 if size(Et, 2) == 1 # handle scalar case but within modal simulation PlasmaScalar!(out, Plas, reshape(Et, size(Et, 1))) From 0e5cb825074746d3de87f42e1740d15aec80d2d9 Mon Sep 17 00:00:00 2001 From: John Travers Date: Mon, 6 Sep 2021 11:56:40 +0100 Subject: [PATCH 06/18] Fix interface --- src/Interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index 2cbf5935..bc6bc578 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -513,7 +513,7 @@ function makeplasma!(out, grid::Grid.EnvGrid, gas, plasma::Symbol, pol) throw(DomainError(plasma, "Unknown ionisation rate $plasma.")) end Et = pol ? Array{Complex{Float64}}(undef, length(grid.to), 2) : Vector{Complex{Float64}}(undef, length(grid.to)) - push!(out, Nonlinear.PlasmaCumtrapz(grid.to, Et, ionrate, ionpot)) + push!(out, Nonlinear.PlasmaFourier(grid.ωo, Et, ionrate, ionpot, grid.to[2] - grid.to[1])) end function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol) From 01569769b91df6931b79167b3482382184db242b Mon Sep 17 00:00:00 2001 From: John Travers Date: Mon, 6 Sep 2021 12:37:46 +0100 Subject: [PATCH 07/18] comment on method --- src/Nonlinear.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 5cb206d7..f563b88d 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -142,6 +142,9 @@ function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E) Maths.cumtrapz!(out, Plas.J, Plas.δt) end +# This version, for envelope fields does not work. Not sure why. The exact same +# equations, with cumulative integration performed by FFT (see PlasmaFourier +# version below) does work. "The plasma response for a scalar electric field" function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E::Vector{Complex{Float64}}) Plas.ratefunc(Plas.rate, E) @@ -207,11 +210,20 @@ end "The plasma response for a scalar electric field" function PlasmaScalar!(out, Plas::PlasmaFourier, E::Vector{Complex{Float64}}) + # The equations used here for envelopes are very similar to the real field cases above. + # So far I have not been able to derive them directly. The modified loss term (Jabs) + # is taken from [1]. The plasma term is obtained by inspection. Essentially, you can start + # from the definition of the plasma refractive index, insert it into the general relation + # between polarisation and refractive index, and obtain this version (in the frequency domain) + # but with the electron density in the time domain(!). I also think this term could be obtained + # from Geissler's model [2], by making use of the Bedrosian identity for Hilbert transforms of + # products of functions, but I have not yet achieved it. Plas.ratefunc(Plas.rate, E) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E # (E + conj.(E)) + # here we perform cumulative integration via the Fourier transform. out .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) end From fbd259402aced85a383af0e7e50cf22ed87e8b55 Mon Sep 17 00:00:00 2001 From: John Travers Date: Tue, 7 Sep 2021 14:29:29 +0100 Subject: [PATCH 08/18] PlasmaFourier for real fields --- src/Nonlinear.jl | 68 +++++++++++++++++++++++++++++------------------- 1 file changed, 41 insertions(+), 27 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index f563b88d..6edd0e20 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -128,7 +128,7 @@ function PlasmaCumtrapz(t, E, ratefunc, ionpot) end "The plasma response for a scalar electric field" -function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E) +function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E::Vector{Float64}) Plas.ratefunc(Plas.rate, E) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @@ -142,20 +142,6 @@ function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E) Maths.cumtrapz!(out, Plas.J, Plas.δt) end -# This version, for envelope fields does not work. Not sure why. The exact same -# equations, with cumulative integration performed by FFT (see PlasmaFourier -# version below) does work. -"The plasma response for a scalar electric field" -function PlasmaScalar!(out, Plas::PlasmaCumtrapz, E::Vector{Complex{Float64}}) - Plas.ratefunc(Plas.rate, E) - Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = 1-exp(-Plas.fraction) - @. Plas.phase = Plas.fraction * e_ratio * E - Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) - Plas.J .+= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E - Maths.cumtrapz!(out, Plas.J, Plas.δt) -end - """ The plasma response for a vector electric field. @@ -196,15 +182,34 @@ struct PlasmaFourier{R, EType, tType, FTType} <: PlasmaPolar δt::Float64 end +""" + make_fft(E::Vector) + +Plan a suitable FFT for the electric field type `E`. + +""" +function make_fft(E::Vector{<:Real}) + Utils.loadFFTwisdom() + FT = FFTW.plan_rfft(E, 1, flags=Luna.settings["fftw_flag"]) + inv(FT) + Utils.saveFFTwisdom() + FT +end + +function make_fft(E::Vector{Complex{<:Real}}) + Utils.loadFFTwisdom() + FT = FFTW.plan_fft(E, 1, flags=Luna.settings["fftw_flag"]) + inv(FT) + Utils.saveFFTwisdom() + FT +end + function PlasmaFourier(ω, E, ratefunc, ionpot, δt) rate = similar(ω) fraction = similar(ω) phase = similar(E) J = similar(E) - Utils.loadFFTwisdom() - FT = FFTW.plan_fft(E, 1, flags=Luna.settings["fftw_flag"]) - inv(FT) - Utils.saveFFTwisdom() + FT = make_fft(E) PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, ω, FT, δt) end @@ -227,6 +232,21 @@ function PlasmaScalar!(out, Plas::PlasmaFourier, E::Vector{Complex{Float64}}) out .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) end +"The plasma response for a scalar electric field" +function PlasmaScalar!(out, Plas::PlasmaFourier, E::Vector{Float64}) + Plas.ratefunc(Plas.rate, E) + Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) + @. Plas.fraction = 1-exp(-Plas.fraction) + @. Plas.phase = Plas.fraction * e_ratio * E + for ii in eachindex(E) + if abs(E[ii]) > 0 + Plas.J[ii] += Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/E[ii] + end + end + # here we perform cumulative integration via the Fourier transform. + out .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) +end + "Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`." function (Plas::PlasmaPolar)(out, Et) if ndims(Et) > 1 @@ -309,10 +329,7 @@ harmonic generation component of the response. """ function RamanPolarField(t, ht; thg=true) h = zeros(length(t)*2) # note double grid size, see explanation below - Utils.loadFFTwisdom() - FT = FFTW.plan_rfft(h, 1, flags=Luna.settings["fftw_flag"]) - inv(FT) - Utils.saveFFTwisdom() + FT = make_fft(h) hω, Eω2, Pω = gethω!(h, t, ht, FT) E2 = similar(h) E2v = view(E2, 1:length(t)) @@ -331,10 +348,7 @@ using response function `ht`. """ function RamanPolarEnv(t, ht) h = zeros(length(t)*2) # note double grid size, see explanation below - Utils.loadFFTwisdom() - FT = FFTW.plan_fft(h, 1, flags=Luna.settings["fftw_flag"]) - inv(FT) - Utils.saveFFTwisdom() + FT = make_fft(h) hω, Eω2, Pω = gethω!(h, t, ht, FT) E2 = similar(hω) P = similar(hω) From 0d921e9c6b21405c03ad31887fbd76c2f6a42336 Mon Sep 17 00:00:00 2001 From: John Travers Date: Tue, 7 Sep 2021 14:38:05 +0100 Subject: [PATCH 09/18] Prototype of vector envelope plasma --- src/Nonlinear.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 6edd0e20..ff3a74d2 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -227,11 +227,23 @@ function PlasmaScalar!(out, Plas::PlasmaFourier, E::Vector{Complex{Float64}}) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E - Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E # (E + conj.(E)) + Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E # here we perform cumulative integration via the Fourier transform. out .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) end +function PlasmaVector!(out, Plas::PlasmaFourier, E::Array{Complex{Float64},2}) + Ex = E[:,1] + Ey = E[:,2] + Em = @. hypot.(Ex, Ey) + Plas.ratefunc(Plas.rate, Em) + Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) + @. Plas.fraction = 1-exp(-Plas.fraction) + @. Plas.phase = Plas.fraction * e_ratio * E + Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E + out .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) +end + "The plasma response for a scalar electric field" function PlasmaScalar!(out, Plas::PlasmaFourier, E::Vector{Float64}) Plas.ratefunc(Plas.rate, E) From 324fae0c0fbdf634e6087b12a6e91ac5f54e328e Mon Sep 17 00:00:00 2001 From: John Travers Date: Tue, 7 Sep 2021 14:57:24 +0100 Subject: [PATCH 10/18] fix make_fft --- src/Nonlinear.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index ff3a74d2..1a9ef059 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -188,7 +188,7 @@ end Plan a suitable FFT for the electric field type `E`. """ -function make_fft(E::Vector{<:Real}) +function make_fft(E::Array{T,N}) where T<:Real where N Utils.loadFFTwisdom() FT = FFTW.plan_rfft(E, 1, flags=Luna.settings["fftw_flag"]) inv(FT) @@ -196,7 +196,7 @@ function make_fft(E::Vector{<:Real}) FT end -function make_fft(E::Vector{Complex{<:Real}}) +function make_fft(E::Array{Complex{T},N}) where T<:Real where N Utils.loadFFTwisdom() FT = FFTW.plan_fft(E, 1, flags=Luna.settings["fftw_flag"]) inv(FT) From 81de6c74a9d9e61baa30b4e2323165f7786b767a Mon Sep 17 00:00:00 2001 From: John Travers Date: Fri, 10 Sep 2021 14:45:08 +0100 Subject: [PATCH 11/18] fix merge --- src/Interface.jl | 1 + src/Nonlinear.jl | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index 7f54f7d8..c0d17c0d 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -535,6 +535,7 @@ function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol) end end raman && push!(out, Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, gas))) + makeplasma!(out, grid, gas, plasma, pol) Tuple(out) end diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 0657cd12..2e13947b 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -164,6 +164,7 @@ struct PlasmaFourier{R, EType, tType, FTType} <: PlasmaPolar fraction::tType # buffer to hold the ionization fraction phase::EType # buffer to hold the plasma induced (mostly) phase modulation J::EType # buffer to hold the plasma current + P::EType # buffer to hold the plasma polarisation ω::tType # buffer to hold the frequency grid FT::FTType # Fourier transform to use for integrations δt::Float64 @@ -196,8 +197,9 @@ function PlasmaFourier(ω, E, ratefunc, ionpot, δt) fraction = similar(ω) phase = similar(E) J = similar(E) + P = similar(E) FT = make_fft(E) - PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, ω, FT, δt) + PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, P, ω, FT, δt) end "The plasma response for a scalar electric field" From 146f3b49b05ced9fa4a3a0c36033b907e58bf443 Mon Sep 17 00:00:00 2001 From: John Travers Date: Fri, 10 Sep 2021 15:25:08 +0100 Subject: [PATCH 12/18] vector env plasma --- src/Nonlinear.jl | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 2e13947b..d2b09753 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -192,7 +192,17 @@ function make_fft(E::Array{Complex{T},N}) where T<:Real where N FT end -function PlasmaFourier(ω, E, ratefunc, ionpot, δt) +function PlasmaFourier(ω, E::Array{T,N}, ratefunc, ionpot, δt) where T<:Real where N + rate = similar(E) + fraction = similar(E) + phase = similar(E) + J = similar(E) + P = similar(E) + FT = make_fft(E) + PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, P, ω, FT, δt) +end + +function PlasmaFourier(ω, E::Array{Complex{T},N}, ratefunc, ionpot, δt) where T<:Real where N rate = similar(ω) fraction = similar(ω) phase = similar(E) @@ -216,9 +226,9 @@ function PlasmaScalar!(Plas::PlasmaFourier, E::Vector{Complex{Float64}}) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E - Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E + Plas.J .= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E # here we perform cumulative integration via the Fourier transform. - Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) + Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω) end function PlasmaVector!(Plas::PlasmaFourier, E::Array{Complex{Float64},2}) @@ -229,8 +239,8 @@ function PlasmaVector!(Plas::PlasmaFourier, E::Array{Complex{Float64},2}) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E - Jabs = Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ Em .* E - Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) + Plas.J .= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ Em.^2 .* E + Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω) end "The plasma response for a scalar electric field" @@ -239,13 +249,14 @@ function PlasmaScalar!(Plas::PlasmaFourier, E::Vector{Float64}) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E + fill!(Plas.J, 0.0) for ii in eachindex(E) if abs(E[ii]) > 0 Plas.J[ii] += Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/E[ii] end end # here we perform cumulative integration via the Fourier transform. - Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Jabs) ./ Plas.ω) + Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω) end "Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`." From 628dca70f5de33e365ed7fb6c2531bbdd825f6c7 Mon Sep 17 00:00:00 2001 From: John Travers Date: Fri, 17 Sep 2021 21:29:45 +0100 Subject: [PATCH 13/18] make plasma default for envelopes too --- src/Interface.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index 30e181ce..c09d34f4 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -300,7 +300,6 @@ In this case, all keyword arguments except for `λ0` are ignored. - `:PPT` -- include plasma using the PPT ionisation rate. - `true` (default) -- same as `:PPT`. - `false` -- ignore plasma. - Note that plasma is only available for full-field simulations. - `thg::Bool`: Whether to include third-harmonic generation. Defaults to `true` for full-field simulations and to `false` for envelope simulations. If `raman` is `true`, then the following options apply: @@ -329,7 +328,7 @@ function prop_capillary(radius, flength, gas, pressure; pulses=nothing, shotnoise=true, modes=:HE11, model=:full, loss=true, - raman=false, kerr=true, plasma=nothing, + raman=false, kerr=true, plasma=true, rotation=true, vibration=true, saveN=201, filepath=nothing, scan=nothing, scanidx=nothing, filename=nothing, @@ -337,7 +336,6 @@ function prop_capillary(radius, flength, gas, pressure; pol = needpol(polarisation, pulses) || needpol_modes(modes) @info "X+Y polarisation "* (pol ? "required." : "not required.") - plasma = isnothing(plasma) ? !envelope : plasma thg = isnothing(thg) ? !envelope : thg gas = (gas == :He) ? :HeJ : gas From 26f44e98f2da5dacb63855fed832284cf8cae2d3 Mon Sep 17 00:00:00 2001 From: John Travers Date: Fri, 17 Sep 2021 21:31:17 +0100 Subject: [PATCH 14/18] fix example --- examples/simple_interface/PlasmaSSFBS.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/simple_interface/PlasmaSSFBS.jl b/examples/simple_interface/PlasmaSSFBS.jl index c76a8985..bd303cc9 100644 --- a/examples/simple_interface/PlasmaSSFBS.jl +++ b/examples/simple_interface/PlasmaSSFBS.jl @@ -11,7 +11,7 @@ pressure = 5.0 # gas pressure in bar energy = 1.7e-6 # energy in the pump pulse pssfbs = prop_capillary(radius, flength, gas, pressure; λ0, τfwhm, energy, - trange=1.6e-12, λlims=(300e-9, 1.8e-6), loss=false, envelope=true, plasma=true) + trange=1.6e-12, λlims=(300e-9, 1.8e-6), loss=false, envelope=true) Plotting.prop_2D(pssfbs, :λ; modes=:sum, trange=(-700e-15, 100e-15), λrange=(300e-9, 1800e-9), dBmin=-30) From 6596fe906c61b50f5586258c683290f901b1c6b9 Mon Sep 17 00:00:00 2001 From: John Travers Date: Sat, 18 Sep 2021 22:15:34 +0100 Subject: [PATCH 15/18] tidy --- src/Nonlinear.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index d2b09753..67ff05dc 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -141,7 +141,7 @@ A similar approach was used in: C Tailliez et al 2020 New J. Phys. 22 103038. function PlasmaVector!(Plas::PlasmaCumtrapz, E) Ex = E[:,1] Ey = E[:,2] - Em = @. hypot.(Ex, Ey) + Em = hypot.(Ex, Ey) Plas.ratefunc(Plas.rate, Em) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @@ -234,7 +234,7 @@ end function PlasmaVector!(Plas::PlasmaFourier, E::Array{Complex{Float64},2}) Ex = E[:,1] Ey = E[:,2] - Em = @. hypot.(Ex, Ey) + Em = hypot.(Ex, Ey) Plas.ratefunc(Plas.rate, Em) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) From e026681e9874ed7cd92e44b11d9286d4c14fe9e2 Mon Sep 17 00:00:00 2001 From: John Travers Date: Tue, 21 Sep 2021 21:19:41 +0100 Subject: [PATCH 16/18] clean up --- src/Interface.jl | 3 +- src/Nonlinear.jl | 160 +++++++++++++++++------------------------------ 2 files changed, 60 insertions(+), 103 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index c09d34f4..2368b0ed 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -530,7 +530,8 @@ function makeplasma!(out, grid::Grid.EnvGrid, gas, plasma::Symbol, pol) else throw(DomainError(plasma, "Unknown ionisation rate $plasma.")) end - Et = pol ? Array{Complex{Float64}}(undef, length(grid.to), 2) : Vector{Complex{Float64}}(undef, length(grid.to)) + pol && error("vector ionisation not suppported with envelopes at present") + Et = Vector{Complex{Float64}}(undef, length(grid.to)) push!(out, Nonlinear.PlasmaFourier(grid.ωo, Et, ionrate, ionpot, grid.to[2] - grid.to[1])) end diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 67ff05dc..84bb724e 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -93,7 +93,8 @@ struct PlasmaCumtrapz{R, EType, tType} <: PlasmaPolar rate::tType # buffer to hold the rate fraction::tType # buffer to hold the ionization fraction phase::EType # buffer to hold the plasma induced (mostly) phase modulation - J::EType # buffer to hold the plasma current + Jloss::EType # buffer to hold the ionization loss current + J::EType # buffer to hold the total plasma current P::EType # buffer to hold the plasma polarisation δt::Float64 # the time step end @@ -101,59 +102,61 @@ end """ PlasmaCumtrapz(t, E, ratefunc, ionpot) -Construct the Plasma polarisation response for a field on time grid `t` -with example electric field like `E`, an ionization rate callable -`ratefunc` and ionization potential `ionpot`. +Construct the Plasma polarisation response based on cumulative integration for a +real field on time grid `t` with example electric field like `E`, an ionization +rate callable `ratefunc` and ionization potential `ionpot`. """ -function PlasmaCumtrapz(t, E, ratefunc, ionpot) +function PlasmaCumtrapz(t, E::Array{T,N}, ratefunc, ionpot) where T<:Real where N rate = similar(t) fraction = similar(t) phase = similar(E) J = similar(E) + Jloss = similar(E) P = similar(E) - return PlasmaCumtrapz(ratefunc, ionpot, rate, fraction, phase, J, P, t[2]-t[1]) + return PlasmaCumtrapz(ratefunc, ionpot, rate, fraction, phase, Jloss, J, P, t[2]-t[1]) end -"The plasma response for a scalar electric field" -function PlasmaScalar!(Plas::PlasmaCumtrapz, E::Vector{Float64}) - Plas.ratefunc(Plas.rate, E) - Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = 1-exp(-Plas.fraction) - @. Plas.phase = Plas.fraction * e_ratio * E - Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) - for ii in eachindex(E) - if abs(E[ii]) > 0 - Plas.J[ii] += Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/E[ii] - end - end - Maths.cumtrapz!(Plas.P, Plas.J, Plas.δt) -end - -""" -The plasma response for a vector electric field. +"ionizing field" +ionE(E) = abs.(E) +ionE(E::Array{T,2}) where T = hypot.(E[:,1], E[:,2]) -We take the magnitude of the electric field to calculate the ionization -rate and fraction, and then solve the plasma polarisation component-wise -for the vector field. +"ionization loss current" +# for envelopes this follows Bergé et al. Rep. Prog. Phys. 70, 1633 (2007) +Jloss!(out, ip, rate, fraction, Em, E) = @. out = ifelse(Em > 0, ip*rate*(1 - fraction)/abs2(Em)*E, 0.0) +# minor optimisation for scalar real field, based on Geissler, PRL 83, 2930 (1999). +Jloss!(out, ip, rate, fraction, Em, E::Vector{Float64}) = @. out = ifelse(Em > 0, ip*rate*(1 - fraction)/E, 0.0) -A similar approach was used in: C Tailliez et al 2020 New J. Phys. 22 103038. """ -function PlasmaVector!(Plas::PlasmaCumtrapz, E) - Ex = E[:,1] - Ey = E[:,2] - Em = hypot.(Ex, Ey) +Calculate the ionization rate, fraction, loss and phase term. +""" +function precalc!(Plas::PlasmaPolar, E) + # This is basd on Geissler, PRL 83, 2930 (1999) for the real scalar field. + # In the vector case we take the magnitude of the electric field to calculate the ionization + # rate and fraction, and then solve the plasma polarisation component-wise for the vector field. + # A similar approach was used in: C Tailliez et al 2020 New J. Phys. 22 103038. + # The equations used here for envelopes are identical to the real field cases (apart from Jloss). + # So far I have not been able to derive them directly. The modified loss term + # is taken from Bergé et al. Rep. Prog. Phys. 70, 1633 (2007). The plasma phase term is obtained + # by inspection. Essentially, you can start from the definition of the plasma refractive index, + # insert it into the general relation between polarisation and refractive index, and obtain this + # version (in the frequency domain) but with the electron density in the time domain(!). + # I also think this term could be obtained from Geissler's model, by making use of the Bedrosian + # identity for Hilbert transforms of products of functions, but I have not yet achieved it. + Em = ionE(E) Plas.ratefunc(Plas.rate, Em) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) @. Plas.fraction = 1-exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E + Jloss!(Plas.Jloss, Plas.ionpot, Plas.rate, Plas.fraction, Em, E) +end + +""" +The plasma response for real electric field. +""" +function calc!(Plas::PlasmaCumtrapz, E::Array{Float64,N}) where N + precalc!(Plas, E) Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) - for ii in eachindex(Em) - if abs(Em[ii]) > 0 - pre = Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/Em[ii]^2 - Plas.J[ii,1] += pre*Ex[ii] - Plas.J[ii,2] += pre*Ey[ii] - end - end + Plas.J .+= Plas.Jloss Maths.cumtrapz!(Plas.P, Plas.J, Plas.δt) end @@ -163,7 +166,8 @@ struct PlasmaFourier{R, EType, tType, FTType} <: PlasmaPolar rate::tType # buffer to hold the rate fraction::tType # buffer to hold the ionization fraction phase::EType # buffer to hold the plasma induced (mostly) phase modulation - J::EType # buffer to hold the plasma current + Jloss::EType # buffer to hold the ionisation loss current + J::EType # buffer to hold the total plasma current P::EType # buffer to hold the plasma polarisation ω::tType # buffer to hold the frequency grid FT::FTType # Fourier transform to use for integrations @@ -176,7 +180,7 @@ end Plan a suitable FFT for the electric field type `E`. """ -function make_fft(E::Array{T,N}) where T<:Real where N +function make_fft(E::Array{T,N}) where {T<:Real, N} Utils.loadFFTwisdom() FT = FFTW.plan_rfft(E, 1, flags=Luna.settings["fftw_flag"]) inv(FT) @@ -184,7 +188,7 @@ function make_fft(E::Array{T,N}) where T<:Real where N FT end -function make_fft(E::Array{Complex{T},N}) where T<:Real where N +function make_fft(E::Array{Complex{T},N}) where {T<:Real, N} Utils.loadFFTwisdom() FT = FFTW.plan_fft(E, 1, flags=Luna.settings["fftw_flag"]) inv(FT) @@ -192,85 +196,37 @@ function make_fft(E::Array{Complex{T},N}) where T<:Real where N FT end -function PlasmaFourier(ω, E::Array{T,N}, ratefunc, ionpot, δt) where T<:Real where N - rate = similar(E) - fraction = similar(E) - phase = similar(E) - J = similar(E) - P = similar(E) - FT = make_fft(E) - PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, P, ω, FT, δt) -end - -function PlasmaFourier(ω, E::Array{Complex{T},N}, ratefunc, ionpot, δt) where T<:Real where N +function PlasmaFourier(ω, E, ratefunc, ionpot, δt) rate = similar(ω) fraction = similar(ω) phase = similar(E) + Jloss = similar(E) J = similar(E) P = similar(E) - FT = make_fft(E) - PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, J, P, ω, FT, δt) -end - -"The plasma response for a scalar electric field" -function PlasmaScalar!(Plas::PlasmaFourier, E::Vector{Complex{Float64}}) - # The equations used here for envelopes are very similar to the real field cases above. - # So far I have not been able to derive them directly. The modified loss term (Jabs) - # is taken from [1]. The plasma term is obtained by inspection. Essentially, you can start - # from the definition of the plasma refractive index, insert it into the general relation - # between polarisation and refractive index, and obtain this version (in the frequency domain) - # but with the electron density in the time domain(!). I also think this term could be obtained - # from Geissler's model [2], by making use of the Bedrosian identity for Hilbert transforms of - # products of functions, but I have not yet achieved it. - Plas.ratefunc(Plas.rate, E) - Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = 1-exp(-Plas.fraction) - @. Plas.phase = Plas.fraction * e_ratio * E - Plas.J .= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ abs2.(E) .* E - # here we perform cumulative integration via the Fourier transform. - Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω) -end - -function PlasmaVector!(Plas::PlasmaFourier, E::Array{Complex{Float64},2}) - Ex = E[:,1] - Ey = E[:,2] - Em = hypot.(Ex, Ey) - Plas.ratefunc(Plas.rate, Em) - Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = 1-exp(-Plas.fraction) - @. Plas.phase = Plas.fraction * e_ratio * E - Plas.J .= Plas.ionpot .* Plas.rate .* (1 .- Plas.fraction) ./ Em.^2 .* E - Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω) + PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, Jloss, J, P, ω, make_fft(E), δt) end -"The plasma response for a scalar electric field" -function PlasmaScalar!(Plas::PlasmaFourier, E::Vector{Float64}) - Plas.ratefunc(Plas.rate, E) - Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = 1-exp(-Plas.fraction) - @. Plas.phase = Plas.fraction * e_ratio * E - fill!(Plas.J, 0.0) - for ii in eachindex(E) - if abs(E[ii]) > 0 - Plas.J[ii] += Plas.ionpot * Plas.rate[ii] * (1-Plas.fraction[ii])/E[ii] - end - end - # here we perform cumulative integration via the Fourier transform. - Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.J) ./ Plas.ω) +""" +The plasma response for a scalar or vector, real field, or scalar envelope, based +on cumulative integration using Fourier transforms. +""" +function calc!(Plas::PlasmaFourier, E) + precalc!(Plas, E) + Plas.P .= Plas.FT \ (-(Plas.FT * Plas.phase) ./ Plas.ω.^2 .- 1im .* (Plas.FT * Plas.Jloss) ./ Plas.ω) end "Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`." function (Plas::PlasmaPolar)(out, Et, ρ) if ndims(Et) > 1 if size(Et, 2) == 1 # handle scalar case but within modal simulation - PlasmaScalar!(Plas, reshape(Et, size(Et,1))) + calc!(Plas, reshape(Et, size(Et,1))) out .+= ρ .* reshape(Plas.P, size(Et)) else - PlasmaVector!(Plas, Et) # vector case + calc!(Plas, Et) # vector case out .+= ρ .* Plas.P end else - PlasmaScalar!(Plas, Et) # straight scalar case + calc!(Plas, Et) # straight scalar case out .+= ρ .* Plas.P end end From 7ab9fe33dfb5d2e95146f5a5e4528ec4a2667587 Mon Sep 17 00:00:00 2001 From: John Travers Date: Tue, 21 Sep 2021 21:25:19 +0100 Subject: [PATCH 17/18] change defaults back --- src/Interface.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index 1a1799da..d09ff1e1 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -298,8 +298,9 @@ In this case, all keyword arguments except for `λ0` are ignored. - `plasma`: Can be one of - `:ADK` -- include plasma using the ADK ionisation rate. - `:PPT` -- include plasma using the PPT ionisation rate. - - `true` (default) -- same as `:PPT`. + - `true` -- same as `:PPT`. - `false` -- ignore plasma. + - `nothing` -- (default), implies `true` for field simulations and `false` for envelopes - `thg::Bool`: Whether to include third-harmonic generation. Defaults to `true` for full-field simulations and to `false` for envelope simulations. If `raman` is `true`, then the following options apply: @@ -343,13 +344,14 @@ function prop_capillary_args(radius, flength, gas, pressure; pulses=nothing, shotnoise=true, modes=:HE11, model=:full, loss=true, - raman=false, kerr=true, plasma=true, + raman=false, kerr=true, plasma=nothing, rotation=true, vibration=true, saveN=201, filepath=nothing, scan=nothing, scanidx=nothing, filename=nothing) pol = needpol(polarisation, pulses) || needpol_modes(modes) @info "X+Y polarisation "* (pol ? "required." : "not required.") + plasma = isnothing(plasma) ? !envelope : plasma thg = isnothing(thg) ? !envelope : thg gas = (gas == :He) ? :HeJ : gas From 6910a9058304556593bd7df65986a455620bb66b Mon Sep 17 00:00:00 2001 From: John Travers Date: Wed, 22 Sep 2021 10:59:22 +0100 Subject: [PATCH 18/18] try to fix tests --- src/Nonlinear.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 84bb724e..0208e856 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -153,7 +153,7 @@ end """ The plasma response for real electric field. """ -function calc!(Plas::PlasmaCumtrapz, E::Array{Float64,N}) where N +function calc!(Plas::PlasmaCumtrapz, E::AbstractArray{Float64,N}) where N precalc!(Plas, E) Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) Plas.J .+= Plas.Jloss @@ -299,7 +299,7 @@ using response function `r`. """ function RamanPolarEnv(t, r) h = zeros(length(t)*2) # note double grid size, see explanation below - FT = make_fft(h) + FT = make_fft(complex.(h)) ht = view(h, 1:length(t)) hω = FT * h Eω2 = similar(hω)