From f87276b356884aafbb048b6f00cee5d2df80c41e Mon Sep 17 00:00:00 2001
From: John Travers <jtravs@gmail.com>
Date: Sun, 5 Sep 2021 23:55:38 +0100
Subject: [PATCH 01/12] 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 <jtravs@gmail.com>
Date: Sun, 5 Sep 2021 23:55:49 +0100
Subject: [PATCH 02/12] 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 <jtravs@gmail.com>
Date: Sun, 5 Sep 2021 23:56:05 +0100
Subject: [PATCH 03/12] 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 <jtravs@gmail.com>
Date: Sun, 5 Sep 2021 23:56:16 +0100
Subject: [PATCH 04/12] 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 <jtravs@gmail.com>
Date: Mon, 6 Sep 2021 11:56:31 +0100
Subject: [PATCH 05/12] 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 <jtravs@gmail.com>
Date: Mon, 6 Sep 2021 11:56:40 +0100
Subject: [PATCH 06/12] 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 <jtravs@gmail.com>
Date: Mon, 6 Sep 2021 12:37:46 +0100
Subject: [PATCH 07/12] 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 <jtravs@gmail.com>
Date: Tue, 7 Sep 2021 14:29:29 +0100
Subject: [PATCH 08/12] 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 <jtravs@gmail.com>
Date: Tue, 7 Sep 2021 14:38:05 +0100
Subject: [PATCH 09/12] 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 <jtravs@gmail.com>
Date: Tue, 7 Sep 2021 14:57:24 +0100
Subject: [PATCH 10/12] 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 <jtravs@gmail.com>
Date: Fri, 10 Sep 2021 14:45:08 +0100
Subject: [PATCH 11/12] 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 <jtravs@gmail.com>
Date: Fri, 10 Sep 2021 15:25:08 +0100
Subject: [PATCH 12/12] 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`."