Add phase integrals

src/QEDfields.jl

@@ -13,9 +13,12 @@ export polarization_vector, oscillator
 
 export CosSquarePulse, GaussianPulse
 
+export phase_integral_0, phase_integral_1, phase_integral_2
+
 include("interfaces/background_field_interface.jl")
 include("polarization.jl")
 include("pulses/cos_square.jl")
 include("pulses/gaussian.jl")
+include("phase_integrals/phase_integrals.jl")
 
 end

src/interfaces/background_field_interface.jl

@@ -2,7 +2,7 @@
 # The abstract background field interface
 #
 # In this file, the abstract interface for different types of background fields
-# is defined. 
+# is defined.
 ####################
 """
 Abstract base type for describing classical background fields.
@@ -45,13 +45,13 @@ function reference_momentum end
 
     domain(::AbstractPulsedPlaneWaveField)
 
-Interface function for [`AbstractPulsedPlaneWaveField`](@ref), which returns interval (as a `IntervalSets.Interval`) for the given background field. 
+Interface function for [`AbstractPulsedPlaneWaveField`](@ref), which returns interval (as a `IntervalSets.Interval`) for the given background field.
 
 """
 function domain end
 
 """
-    
+
     pulse_length(::AbstractPulsedPlaneWaveField)
 
 Interface function for [`AbstractPulsedPlaneWaveField`](@ref), which returns a dimensionless representative number for the duration of the background field,
@@ -68,35 +68,35 @@ Interface function for [`AbstractPulsedPlaneWaveField`](@ref), which returns the
 !!! note "Single point implementation"
 
     The interface function can be implemented for just one phase point as input. With that, evaluation on a vector of inputs is generically implemented by broadcasting.
-    However, if there is a better custom implementation for vectors in input values, consider implementing 
+    However, if there is a better custom implementation for vectors in input values, consider implementing
     ```Julia
-    
+
         _envelope(::AbstractPulsedPlaneWaveField, phi::AbstractVector{T<:Real})
 
     ```
 
 !!! note "unsafe implementation"
-    
-    This is the unsafe version of the phase envelope function, i.e. this should be implement without input checks like the domain check. 
-    In the safe version [`envelope`](@ref), a domain check is performed, i.e. it returns the value of `_envelope` if the passed in `phi` 
-    is in the `domain` of the field, and zero otherwise. 
+
+    This is the unsafe version of the phase envelope function, i.e. this should be implement without input checks like the domain check.
+    In the safe version [`envelope`](@ref), a domain check is performed, i.e. it returns the value of `_envelope` if the passed in `phi`
+    is in the `domain` of the field, and zero otherwise.
 
 """
 function _envelope end
 
 function _envelope(
     field::AbstractPulsedPlaneWaveField, phi::AbstractVector{T}
 ) where {T<:Real}
-    # TODO: maybe use broadcasting here 
+    # TODO: maybe use broadcasting here
     return map(x -> _envelope(field, x), phi)
 end
 
 """
-    
+
     envelope(pulsed_field::AbstractPulsedPlaneWaveField, phi::Real)
-    
-Return the value of the phase envelope function (also referred to as pulse envelope or pulse shape) 
-for given `pulsed_field` and phase `phi`. Performs domain check on `phi` before calling [`_envelope`](@ref); 
+
+Return the value of the phase envelope function (also referred to as pulse envelope or pulse shape)
+for given `pulsed_field` and phase `phi`. Performs domain check on `phi` before calling [`_envelope`](@ref);
 returns zero if `phi` is not in the domain returned by `[domain](@ref)`.
 """
 function envelope(field::AbstractPulsedPlaneWaveField, phi::Real)
@@ -106,7 +106,7 @@ end
 function envelope(
     field::AbstractPulsedPlaneWaveField, phi::AbstractVector{T}
 ) where {T<:Real}
-    # TODO: maybe use broadcasting here 
+    # TODO: maybe use broadcasting here
     return map(x -> envelope(field, x), phi)
 end
 
@@ -123,15 +123,15 @@ function _amplitude(
     pol::AbstractDefinitePolarization,
     phi::AbstractVector{T},
 ) where {T<:Real}
-    # TODO: maybe use broadcasting here 
+    # TODO: maybe use broadcasting here
     return map(x -> _amplitude(field, pol, x), phi)
 end
 
 """
 
     amplitude(field::AbstractPulsedPlaneWaveField, pol::AbstractDefinitePolarization, phi)
 
-Returns the value of the amplitude for a given polarization direction and phase variable `phi`. 
+Returns the value of the amplitude for a given polarization direction and phase variable `phi`.
 
 !!! note "Conventions"
 
@@ -143,7 +143,7 @@ Returns the value of the amplitude for a given polarization direction and phase
     ```
 
 !!! note "Safe implementation"
-    
+
     In this function, a domain check is performed, i.e. if `phi` is in the domain of the field,
     the value of the amplitude is returned, and zero otherwise.
 """
@@ -158,7 +158,7 @@ function amplitude(
     pol::AbstractDefinitePolarization,
     phi::AbstractVector{T},
 ) where {T<:Real}
-    # TODO: maybe use broadcasting here 
+    # TODO: maybe use broadcasting here
     return map(x -> amplitude(field, pol, x), phi)
 end
 
@@ -197,6 +197,6 @@ function generic_spectrum(
     pol::AbstractDefinitePolarization,
     photon_number_parameter::AbstractVector{T},
 ) where {T<:Real}
-    # TODO: maybe use broadcasting here 
+    # TODO: maybe use broadcasting here
     return map(x -> generic_spectrum(field, pol, x), photon_number_parameter)
 end

src/interfaces/phase_integrals.jl

@@ -0,0 +1,135 @@
+####################
+# The abstract phase integral interface
+#
+# In this file, the abstract interface for different computation methods of
+# phase integrals is defined.
+####################
+"""
+Abstract base type for defining a method for phase integral computation.
+"""
+abstract type ComputeMethod end
+
+"""
+Analytic method for phase integral computation.
+
+Requires an existing implementation of analytic formulas for computing the
+entities in the phase integrals.
+"""
+struct Analytic <: ComputeMethod end
+
+"""
+Fully numerical method for phase integral computation.
+"""
+struct Numeric <: ComputeMethod end
+
+"""
+Struct holding setup specific information to compute phase integrals.
+
+ToDo: We mix physical and numerical aspects in this class.
+This does not seem ideal to me (Klaus).
+"""
+struct PhaseIntegral{F<:AbstractBackgroundField, C<:ComputeMethod}
+    bgfield::F
+    method::C
+
+    function PhaseIntegral(bgfield::F, method::C) where {F<:AbstractBackgroundField, C<:ComputeMethod}
+        _assert_compatibility(bgfield, method)
+        new(bgfield, method)
+end
+
+
+"""
+
+    _assert_compatibility(::AbstractBackgroundField, ::ComputeMethod)
+
+Interface function enforcing essential error check whether the respective functions exist that compute the phase integrals of the provided background field by the required method.
+This function represents the default case stating that compute method and background field are not compatible.
+
+It must be specialized for compatible types as
+    _assert_compatibility(::CompatibleField, ::CompatibleComputeMethod) = nothing
+"""
+function _assert_compatibility(<:AbstractBackgroundField, <:ComputeMethod)
+    throw(ArgumentError("phase integral of provided background field cannot be computed with required method!"))
+end
+
+
+"""
+
+    compute(pi::PhaseIntegral, vi::VertexInput)
+
+Interface function for [`PhaseIntegral`](@ref), which returns the result of the computation.
+
+Need to be specialized? I need to check on what the computation of phase integrals actually depends.
+"""
+function compute end
+
+
+"""
+Interface function returning background field.
+All interface functions of background field, e.g. [`reference_momentum()`](@ref) can then be called on it.
+"""
+background_field(pi::PhaseIntegral) = pi.bgfield
+
+
+"""
+Composite type providing all (kinematic) information, required for dressed vertex computation and physe integral computation further down the line, via interface functions
+
+TODO: Should pnum become an explicit dependence of compute(pi, vi) in order to not construct a new VertexInput for every l, which is probably very slow. Because in the end we want to integrate the dressed vertex function Gamma^\mu over pnum.
+"""
+struct VertexInput{M<:QEDbase.AbstractFourMomentum,T<:Real}
+    pin::M
+    pout::M
+    pnum::T
+end
+
+
+"""
+Interface functions of VertexInput
+"""
+in_momentum(vi::VertexInput) = vi.pin
+out_momentum(vi::VertexInput) = vi.pout
+photon_number_parameter(vi::VertexInput) = vi.pnum
+
+
+
+# TODO: The follwing is old stuff from 2024 and should be removed if not required in the implementation started in May 2025
+# phase integrals B_i(l)
+#
+# TODO: Up to now, all of this is just copy paste and needs to be adapted to phase integrals
+"""
+
+    computeB1(ph_int_stp::PhaseIntegral, pol::AbstractPolarization, a0, pnum, alpha1x, alpha1y, alpha2)
+
+Return the first phase integral for the given setup `ph_Int_stp`, background field strength `a0`, photon number parameter `pnum`, components of the kinematic vector factor ``\\alpha_1^\\mu``, and kinematic scalar factor ``\\alpha_2``.
+
+!!! note "Convention"
+
+    The first phase integral is defined as:
+
+    ```math
+    \\begin{align*}
+        B_1^\\mu(l, p, p^\\prime)& = \\int \\mathrm{d}\\varphi A^\\mu(\\varphi)\\exp[\\imath l \\varphi + \\imath G(\\varphi)] \\\\
+    \\end{align*}
+    ```
+    where ``A^\\mu(\\varphi)`` is the background field, ``G(\\varphi,p, p^\\prime)`` is the [`phase function`](@ref), ``(p,p^\\prime)`` the given phase space point, and ``l`` the photon number parameter.
+"""
+function computeB1 end
+
+# TODO: REWORK THE FOLLOWING DOCUMENTATION
+"""
+
+    computeB2()
+
+Return the second phase integral.
+
+!!! note "Convention"
+
+    The second phase integral is defined as:
+
+    ```math
+    \\begin{align*}
+        B_2(l, p, p^\\prime)& = \\int \\mathrm{d}\\varphi A(\\varphi)^2 \\exp[\\imath l \\varphi + \\imath G(\\varphi)] \\\\
+    \\end{align*}
+    ```
+"""
+function computeB2 end

src/phase_integrals/first_integral.jl

@@ -0,0 +1,26 @@
+######################
+# First phase integral
+######################
+
+# Does it need to be the same T for all arguments?
+function phase_integral_1(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    p_in::T,
+    p_out::T,
+    pnum::T
+) where {T<:Real}
+    return quadgk(t -> amplitude(field, pol, t)*_shared_integrand(field, pol, p_in, p_out, t, pnum), endpoints(domain(field))...)[1]
+end
+
+function phase_integral_1(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    p_in::T,
+    p_out::T,
+    photon_number_parameter::AbstractVector{T}
+) where {T<:Real}
+    # TODO: maybe use broadcasting here
+    return map(x -> phase_integral_1(field, pol, p_in, p_out, x), photon_number_parameter)
+end
+