diff --git a/examples/simple_interface/PlasmaSSFBS.jl b/examples/simple_interface/PlasmaSSFBS.jl new file mode 100644 index 00000000..bd303cc9 --- /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) + +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)) diff --git a/src/Interface.jl b/src/Interface.jl index 6cf83cc1..472d9bba 100644 --- a/src/Interface.jl +++ b/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 diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index 4de039d6..0208e856 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -83,15 +83,18 @@ 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 @@ -99,74 +102,131 @@ 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ω) diff --git a/src/Stats.jl b/src/Stats.jl index a4048a21..40f548b2 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)