Cos square pulse (#24)

In this PR, we add an implementation of the background field interface
for the cos-square pulse. Unit tests are adjusted accordingly.

Coincidently, a docstring in the interface was fixed as well.

---------

Co-authored-by: Uwe Hernandez Acosta <[email protected]>

src/QEDfields.jl

@@ -19,7 +19,10 @@ export reference_momentum, domain, phase_duration, envelope, amplitude, generic_
 
 export polarization_vector, oscillator
 
+export CosSquarePulse
+
 include("interfaces/background_field_interface.jl")
 include("polarization.jl")
+include("pulses/cos_square_pulse.jl")
 
 end

src/interfaces/background_field_interface.jl

@@ -178,8 +178,10 @@ Return the generic spectrum of the given field, for the given polarization direc
     The generic spectrum is defined as the Fourier transform of the respective amplitude function for the given polarization direction:
 
     ```math
-    x-\\mathrm{pol} \\to \\int_{-\\infty}^{\\infty} g(\\phi) \\cos(\\phi) \\exp{il\\phi}
-    y-\\mathrm{pol} \\to \\int_{-\\infty}^{\\infty} g(\\phi) \\sin(\\phi) \\exp{il\\phi}
+    \\begin{align*}
+        x-\\mathrm{pol} &\\to \\int_{-\\infty}^{\\infty} g(\\phi) \\cos(\\phi) \\exp(il\\phi)\\\\
+        y-\\mathrm{pol} &\\to \\int_{-\\infty}^{\\infty} g(\\phi) \\sin(\\phi) \\exp(il\\phi)
+    \\end{align*}
     ```
     where ``g(\\phi)`` is the [`envelope`](@ref) and ``l`` the photon number parameter.
 """

src/pulses/cos_square_pulse.jl

@@ -0,0 +1,70 @@
+#####################
+# Cosine square pulse
+#####################
+
+"""
+
+    CosSquarePulse(mom::M,pulse_width::T) where {M<:QEDbase.AbstractFourMomentum,T<:Real}
+
+Concrete implementation of an `AbstractPulsedPlaneWaveField` for cos-square pulses.
+
+!!! note "Pulse shape"
+
+    The pulse envelope of a cos-square pulse is defined as 
+
+    ```math
+    g(\\phi) = \\cos^2(\\frac{\\pi\\phi}{2\\Delta\\phi})
+    ```
+    for \$\\phi\\in (-\\Delta\\phi,\\Delta\\phi),\$where \$\\Delta\\phi`\$ denotes the `pulse_width`, and zero otherwise.
+
+"""
+struct CosSquarePulse{M<:QEDbase.AbstractFourMomentum,T<:Real} <:
+       AbstractPulsedPlaneWaveField
+    mom::M
+    pulse_width::T
+end
+
+@inline function _unsafe_cos_square_envelope(phi, dphi)
+    return cos(pi * phi / (2 * dphi))^2
+end
+
+function _cos_square_envelope(x, dphi)
+    return -dphi <= x <= dphi ? unsafe_envelope(x, dphi) : zero(x)
+end
+
+####
+# interface functions
+####
+
+reference_momentum(pulse::CosSquarePulse) = pulse.mom
+function domain(pulse::CosSquarePulse)
+    delta_phi = pulse.pulse_width
+    return Interval(-delta_phi, delta_phi)
+end
+
+phase_duration(pulse::CosSquarePulse) = pulse.pulse_width
+function _envelope(pulse::CosSquarePulse, phi::Real)
+    return _unsafe_cos_square_envelope(phi, pulse.pulse_width)
+end
+
+#######
+# Special implementation for the generic spectrum
+#######
+
+@inline function _gsinc(x::T) where {T<:Real}
+    return abs(x) == 1 ? one(x) / 2 : sinc(x) / (1 - x^2)
+end
+
+@inline function _generic_FT(l::Real, sig::Real)
+    return sig * _gsinc(sig * l / pi)
+end
+
+function generic_spectrum(field::CosSquarePulse, pol::PolX, pnum::Real)
+    dphi = field.pulse_width
+    return 0.5 * (_generic_FT(pnum + 1, dphi) + _generic_FT(pnum - 1, dphi))
+end
+
+function generic_spectrum(field::CosSquarePulse, pol::PolY, pnum::Real)
+    dphi = field.pulse_width
+    return -0.5im * (_generic_FT(pnum + 1, dphi) - _generic_FT(pnum - 1, dphi))
+end

test/Project.toml

@@ -2,6 +2,7 @@
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/pulses/cos_square_pulse.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 CosSquarePulseWrapper{C<:CosSquarePulse} <: AbstractPulsedPlaneWaveField
+    pulse::C
+end
+QEDfields.reference_momentum(p::CosSquarePulseWrapper) = reference_momentum(p.pulse)
+QEDfields.domain(p::CosSquarePulseWrapper) = domain(p.pulse)
+QEDfields.phase_duration(p::CosSquarePulseWrapper) = phase_duration(p.pulse)
+QEDfields._envelope(p::CosSquarePulseWrapper, x) = QEDfields._envelope(p.pulse, x)
+
+@testset "pulse interface" begin
+    @test hasmethod(reference_momentum, Tuple{CosSquarePulse})
+    @test hasmethod(domain, Tuple{CosSquarePulse})
+    @test hasmethod(phase_duration, Tuple{CosSquarePulse})
+    @test hasmethod(QEDfields._envelope, Tuple{CosSquarePulse,Real})
+end
+@testset "dphi: $dphi" for dphi in DPHIS
+    test_mom = rand(RNG, SFourMomentum)
+    test_pulse = CosSquarePulse(test_mom, dphi)
+
+    @testset "properties" begin
+        @test reference_momentum(test_pulse) == test_mom
+        @test domain(test_pulse) == Interval(-dphi, dphi) # cos_square specific
+        @test phase_duration(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 = CosSquarePulseWrapper(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

@@ -13,3 +13,7 @@ end
 @time @safetestset "polarization" begin
     include("polarization.jl")
 end
+
+@time @safetestset "pulses" begin
+    include("pulses/cos_square_pulse.jl")
+end