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,75 @@
+####################
+# The abstract phase integral interface
+#
+# In this file, the abstract interface for different copmutation methods of
+# phase integrals is defined.
+####################
+"""
+Abstract base type for defining a method for phase integral computation.
+"""
+abstract type Method end
+
+"""
+Analytic method for phase integral computation.
+
+Requires an existing implementation of analytic formulas for computing the
+entities in the phase integrals.
+"""
+struct Analytic <: Method end
+
+"""
+Fully numerical method for phase integral computation based on QuadGK.
+"""
+struct QuadGK <: Method 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{P<:AbstractPulsedPlaneWaveField, M<:Method}
+    pulse::P
+    method::M
+end
+
+# 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
+

src/phase_integrals/phase_integrals.jl

@@ -0,0 +1,119 @@
+###########################################
+# Definitions common to all phase integrals
+###########################################
+
+# TODO: We need to talk about the handling of field polarization!
+# Since some of the quantities required in phase integral computation are four-vectors, just as the bg-field itself,
+# we need to return these as four-vectors, or their components.
+# However, we do not even provide the field as a four-vector, yet.
+# We just return its "amplitude", which is only the oscillator times the envelope,
+# not even taking the correct maximum value a_0 of the field into account. (TODO!)
+#
+# I think both the polarization and the maximum value of the field should be a user defined quantity,
+# but then these should also be members of the field struct, shouldn't they?
+# Or are these separately set in the Process?
+#
+# All in all, I wonder how we should treat four-vectors in QEDfields.
+# The problem of missing vectorial information appears here in _field_integral(), _kinematic_vector_phase_factor(), and phase_integral_1().
+
+# TODO: The factors and integrals should not depend on the momenta of particles
+# but on the Process and Phase Space Point (the latter of which holds the momenta)
+# Question: Does the process hold a reference to the background field?
+#   If so, the explicit dependence on the background field should be removed to.
+
+# from QEDprocesses.jl/src/constants.jl
+# TODO: we might want to move the constants.jl file to QEDbase? See also TODO in QEDprocesses.jl/src/processes/one_photon_compton/perturbative/cross_section.jl
+const ALPHA = inv(137.035999074)
+const ELEMENTARY_CHARGE = sqrt(4 * pi * ALPHA)
+const ELEMENTARY_CHARGE_SQUARE = 4 * pi * ALPHA
+
+# definite integral of PPW field
+# TODO: Why does integration always start at 0?
+@inline function _field_integral(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    phi::T
+) where {T<:Real}
+    return quadgk(t -> amplitude(field, pol, t), 0.0, phi)[1]
+end
+
+# definite integral of squared PPW field
+# TODO: Why does integration always start at 0?
+@inline function _field_squared_integral(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    phi::T
+) where {T<:Real}
+    return quadgk(t -> amplitude(field, pol, t)^2, 0.0, phi)[1]
+end
+
+# kinematic vector factor alpha_1^mu appearing in Volkov phase
+# TODO: Is `ELEMENTARY_CHARGE` actually related to the scattering particle or just a factor
+@inline function _kinematic_vector_phase_factor(
+    field::AbstractPulsedPlaneWaveField,
+    p_in::MT,
+    p_out::MT,
+) where {MT<:AbstractFourMomentum}
+    k_mu = momentum(field)
+    return ELEMENTARY_CHARGE*(p_out/(k_mu*p_out) - p_in/(k_mu*p_in))
+end
+
+# kinematic scalar factor alpha_2 appearing in Volkov phase
+# TODO: Is `ELEMENTARY_CHARGE` actually related to the scattering particle or just a factor
+@inline function _kinematic_scalar_phase_factor(
+    field::AbstractPulsedPlaneWaveField,
+    p_in::MT,
+    p_out::MT,
+) where {MT<:AbstractFourMomentum}
+    k_mu = momentum(field)
+    return ELEMENTARY_CHARGE_SQUARE * (1/(k_mu*p_in) - 1/(k_mu*p_out))
+end
+
+# non-linear Volkov phase
+"""
+
+    _phase_function(field::AbstractPulsedPlaneWaveField, pol::AbstractPolarization, p_in::, p_out::, phi::)
+
+Return the phase function (or non-linear Volkov phase), for the given phase space point `p_in, p_out` and a given phase value `phi`.
+
+!!! note "Convention"
+
+    The non-linear Volkov phase is defined as:
+
+    ```math
+    \\begin{align*}
+        G(\\varphi,p, p^\\prime)& = \\alpha_1^\\mu \\int\\limits_0^\\varphi \\mathrm{d}\\varphi^\\prime A_\\mu(\\varphi^\\prime)
+            + \\alpha_2 \\int\\limits_0^\\varphi \\mathrm{d}\\varphi^\\prime A^2(\\varphi^\\prime) \\\\
+    \\end{align*}
+    ```
+    where ``A^\\mu(\\varphi)`` is the background field, ``\\alpha_1^\\mu`` is the [`kinematic vector phase factor`](@ref), and ``\\alpha_2`` is the [`kinematic scalar phase factor`](@ref).
+    Both ``\\alpha_1^\\mu`` and ``\\alpha_2`` depend on the given phase space point ``(p,p^\\prime)`` and the field's reference momentum ``k^\\mu`` the photon number parameter.
+"""
+@inline function _phase_function(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    p_in::MT,
+    p_out::MT,
+    phi::T
+) where {MT<:AbstractFourMomentum, T<:Real}
+    first = _kinematic_vector_phase_factor(field, p_in, p_out) * _field_integral(field, pol, phi)
+    second = _kinematic_scalar_phase_factor(field, p_in, p_out) * _field_squared_integral(field, pol, phi)
+    return first+second
+end
+
+# integrand shared by all phase integrals
+@inline function _shared_integrand(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    p_in::MT,
+    p_out::MT,
+    phi::T,
+    pnum::T
+) where {MT<:AbstractFourMomentum, T<:Real}
+    return exp(1im*(pnum*phi + _phase_function(field, pol, p_in, p_out, phi)))
+end
+
+
+include("zeroth_integral.jl")
+include("first_integral.jl")
+include("second_integral.jl")

src/phase_integrals/second_integral.jl

@@ -0,0 +1,26 @@
+#######################
+# Second phase integral
+#######################
+
+# Does it need to be the same T for all arguments?
+function phase_integral_2(
+    field::AbstractPulsedPlaneWaveField,
+    pol::AbstractPolarization,
+    p_in::T,
+    p_out::T,
+    pnum::T
+) where {T<:Real}
+    return quadgk(t -> amplitude(field, pol, t)^2 * _shared_integrand(field, pol, p_in, p_out, t, pnum), endpoints(domain(field))...)[1]
+end
+
+function phase_integral_2(
+    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_2(field, p_in, p_out, x), photon_number_parameter)
+end
+