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)) diff --git a/src/Interface.jl b/src/Interface.jl index 997d4d6f..c0d17c0d 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -509,8 +509,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.PlasmaFourier(grid.ωo, Et, ionrate, ionpot, grid.to[2] - grid.to[1])) +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 @@ -523,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 4de039d6..d2b09753 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -83,9 +83,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 @@ -113,7 +115,7 @@ function PlasmaCumtrapz(t, E, ratefunc, ionpot) end "The plasma response for a scalar electric field" -function PlasmaScalar!(Plas::PlasmaCumtrapz, E) +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) @@ -155,8 +157,110 @@ function PlasmaVector!(Plas::PlasmaCumtrapz, E) 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 + 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 +end + +""" + make_fft(E::Vector) + +Plan a suitable FFT for the electric field type `E`. + +""" +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) + Utils.saveFFTwisdom() + FT +end + +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) + Utils.saveFFTwisdom() + 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 + rate = similar(ω) + fraction = similar(ω) + phase = 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.ω) +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.ω) +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))) @@ -217,11 +321,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 +343,8 @@ using response function `r`. """ function RamanPolarEnv(t, r) 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_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 6500b6bb..94c4b7e9 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)