Add Gaussian Pulsed Plane Wave field

src/QEDfields.jl

@@ -11,10 +11,11 @@ export reference_momentum, domain, pulse_length, envelope, amplitude, generic_sp
 
 export polarization_vector, oscillator
 
-export CosSquarePulse
+export CosSquarePulse, GaussianPulse
 
 include("interfaces/background_field_interface.jl")
 include("polarization.jl")
 include("pulses/cos_square.jl")
+include("pulses/gaussian.jl")
 
 end

src/pulses/gaussian.jl

@@ -0,0 +1,59 @@
+############################
+# Gaussian pulsed plane wave
+############################
+
+"""
+    GaussianPulse(mom::M,pulse_length::T) where {M<:QEDbase.AbstractFourMomentum,T<:Real}
+
+Concrete implementation of an `AbstractPulsedPlaneWaveField` for Gaussian pulses.
+
+Propagates along the direction given by the space-like components of the reference momentum vector k_mu = (omega/speed_of_light, k_x, k_y, k_z), and omega = 2*pi*speed_of_light/wavelength.
+
+The time-like coordinate omega/speed_of_ligth of k_mu defines the field's oscillation frequency.
+
+In order to fulfill the vacuum dispersion relation, k_mu*k^mu=0 is required.
+
+!!! note "Pulse shape"
+
+    The longitudinal envelope of a Gaussian pulse is defined as
+
+    ```math
+    g(\\phi) = \\exp(-\\frac{\\phi^2}{2\\Delta\\phi^2})
+    ```
+    for \$\\phi\\in (-\\infty, \\infty) and \$\\Delta\\phi`\$ denotes the `pulse_length`.
+
+    There is no envelope in the transverse directions.
+
+"""
+struct GaussianPulse{M<:QEDbase.AbstractFourMomentum,T<:Real} <:
+       AbstractPulsedPlaneWaveField
+    mom::M
+    pulse_length::T
+end
+
+"""
+
+    _unsafe_gaussian_envelope(phi::Real, dphi::Real)
+
+The envelope of the Gaussian background field is defined according to the standard Gaussian distribution with `dphi` representing the distribution's standard deviation.
+"""
+@inline function _unsafe_gaussian_envelope(phi, dphi)
+    return exp(-0.5 * (phi / dphi)^2)
+end
+
+####
+# interface functions
+####
+
+reference_momentum(pulse::GaussianPulse) = pulse.mom
+
+function domain(pulse::GaussianPulse)
+    dphi = pulse.pulse_length
+    return Interval(-Inf, Inf)
+end
+
+pulse_length(pulse::GaussianPulse) = pulse.pulse_length
+
+function _envelope(pulse::GaussianPulse, phi::Real)
+    return _unsafe_gaussian_envelope(phi, pulse.pulse_length)
+end

test/pulses/gaussian.jl

@@ -0,0 +1,61 @@
+using Random
+using IntervalSets
+using QEDfields
+using QEDbase
+using QuadGK
+
+RNG = MersenneTwister(123456789)
+ATOL = 0.0
+RTOL = sqrt(eps())
+
+DPHIS = [rand(RNG), rand(RNG) * 10, rand(RNG) * 100, rand(RNG) * 1000, rand(RNG) * 10000]
+
+# wrapper implementation to test analytical solutions of the generic spectrum
+
+struct GaussianPulseWrapper{G<:GaussianPulse} <: AbstractPulsedPlaneWaveField
+    pulse::G
+end
+QEDfields.reference_momentum(p::GaussianPulseWrapper) = reference_momentum(p.pulse)
+QEDfields.domain(p::GaussianPulseWrapper) = domain(p.pulse)
+QEDfields.pulse_length(p::GaussianPulseWrapper) = pulse_length(p.pulse)
+QEDfields._envelope(p::GaussianPulseWrapper, x) = QEDfields._envelope(p.pulse, x)
+
+@testset "pulse interface" begin
+    @test hasmethod(reference_momentum, Tuple{GaussianPulse})
+    @test hasmethod(domain, Tuple{GaussianPulse})
+    @test hasmethod(pulse_length, Tuple{GaussianPulse})
+    @test hasmethod(QEDfields._envelope, Tuple{GaussianPulse,Real})
+end
+@testset "dphi: $dphi" for dphi in DPHIS
+    test_mom = rand(RNG, SFourMomentum)
+    test_pulse = GaussianPulse(test_mom, dphi)
+
+    @testset "properties" begin
+        @test reference_momentum(test_pulse) == test_mom
+        @test domain(test_pulse) == Interval(-Inf, Inf)
+        @test pulse_length(test_pulse) == dphi
+    end
+
+    @testset "envelope" begin
+        # unity at the origin
+        @test envelope(test_pulse, 0.0) == 1.0
+
+        # zero at the endpoints
+        @test isapprox(envelope(test_pulse, -Inf), 0.0, atol=ATOL, rtol=RTOL)
+        @test isapprox(envelope(test_pulse, Inf), 0.0, atol=ATOL, rtol=RTOL)
+    end
+    @testset "generic spectrum" begin
+        wrapper_pulse = GaussianPulseWrapper(test_pulse)
+        test_pnums = [1.0, -1.0, 1 + rand(RNG) * 0.1, -1 - rand(RNG) * 0.1]
+        @testset "pnum: $pnum" for pnum in test_pnums
+            test_val_xpol = generic_spectrum(test_pulse, PolX(), pnum)
+            test_val_ypol = generic_spectrum(test_pulse, PolY(), pnum)
+
+            groundtruth_xpol = generic_spectrum(wrapper_pulse, PolX(), pnum)
+            groundtruth_ypol = generic_spectrum(wrapper_pulse, PolY(), pnum)
+
+            @test isapprox(test_val_xpol, groundtruth_xpol, atol=ATOL, rtol=RTOL)
+            @test isapprox(test_val_ypol, groundtruth_ypol, atol=ATOL, rtol=RTOL)
+        end
+    end
+end

test/runtests.jl

@@ -17,4 +17,5 @@ end
 
 @time @safetestset "pulses" begin
     include("pulses/cos_square.jl")
+    include("pulses/gaussian.jl")
 end