Envelope plasma take 2

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)
+
+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))

src/Interface.jl

@@ -298,9 +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.
-    Note that plasma is only available for full-field simulations.
+    - `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:
@@ -550,9 +550,22 @@ 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
+    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
+
 function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol,
                       rotation, vibration)
-    plasma && error("Plasma response for envelope fields has not been implemented yet.")
     isnothing(thg) && (thg = false) 
     out = Any[]
     if kerr
@@ -572,6 +585,7 @@ function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol,
         rr = Raman.raman_response(grid.to, gas, rotation=rotation, vibration=vibration)
         push!(out, Nonlinear.RamanPolarEnv(grid.to, rr))
     end
+    makeplasma!(out, grid, gas, plasma, pol)
     Tuple(out)
 end
 

src/Nonlinear.jl

@@ -83,90 +83,150 @@ 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
     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
 
 """
     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)
-    Plas.ratefunc(Plas.rate, E)
+"ionizing field"
+ionE(E) = abs.(E)
+ionE(E::Array{T,2}) where T = hypot.(E[:,1], E[:,2])
+
+"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)
+
+"""
+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::AbstractArray{Float64,N}) where N
+    precalc!(Plas, 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
+    Plas.J .+= Plas.Jloss
     Maths.cumtrapz!(Plas.P, 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
+    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
+    δt::Float64
+end
+
 """
-The plasma response for a vector electric field.
+    make_fft(E::Vector)
 
-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.
+Plan a suitable FFT for the electric field type `E`.
 
-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)
-    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
-    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
-    Maths.cumtrapz!(Plas.P, Plas.J, Plas.δt)
+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)
+    Utils.saveFFTwisdom()
+    FT
+end
+
+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)
+    Utils.saveFFTwisdom()
+    FT
+end
+
+function PlasmaFourier(ω, E, ratefunc, ionpot, δt)
+    rate = similar(ω)
+    fraction = similar(ω)
+    phase = similar(E)
+    Jloss = similar(E)
+    J = similar(E)
+    P = similar(E)
+    PlasmaFourier(ratefunc, ionpot, rate, fraction, phase, Jloss, J, P, ω, make_fft(E), δt)
+end
+
+"""
+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::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!(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
@@ -217,11 +277,8 @@ harmonic generation component of the response.
 """
 function RamanPolarField(t, r; thg=true)
     h = zeros(length(t)*2) # note double grid size, see explanation below
+    FT = make_fft(h)
     ht = view(h, 1:length(t))
-    Utils.loadFFTwisdom()
-    FT = FFTW.plan_rfft(h, 1, flags=Luna.settings["fftw_flag"])
-    inv(FT)
-    Utils.saveFFTwisdom()
     hω = FT * h
     Eω2 = similar(hω)
     Pω = similar(hω)
@@ -242,11 +299,8 @@ using response function `r`.
 """
 function RamanPolarEnv(t, r)
     h = zeros(length(t)*2) # note double grid size, see explanation below
+    FT = make_fft(complex.(h))
     ht = view(h, 1:length(t))
-    Utils.loadFFTwisdom()
-    FT = FFTW.plan_fft(h, 1, flags=Luna.settings["fftw_flag"])
-    inv(FT)
-    Utils.saveFFTwisdom()
     hω = FT * h
     Eω2 = similar(hω)
     Pω = similar(hω)

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)