Fix indefinite pols

src/QEDfields.jl

@@ -50,6 +50,7 @@ using QEDbase
 
 include("patch_QEDcore.jl")
 include("results.jl")
+include("polarization.jl")
 
 include("integration_methods/interface.jl")
 include("integration_methods/gauss_kronrod.jl")
@@ -63,6 +64,8 @@ include("generic/internal_integrals.jl")
 include("generic/volkov_phase.jl")
 include("generic/phase_integrals.jl")
 
+include("kinematic_factors.jl")
+
 include("pulse_profiles/interface.jl")
 include("pulse_profiles/generic.jl")
 include("pulse_profiles/cos_square/impl.jl")

src/generic/amplitude.jl

@@ -4,6 +4,7 @@
 # delegations
 oscillator(field::AbstractPlaneWaveField, phi::Real) = oscillator(polarization(field), phi)
 polarization_vector(field::AbstractPlaneWaveField) = polarization_vector(polarization(field), reference_momentum(field))
+polarization_vector(field::AbstractPlaneWaveField, pol::AbstractDefinitePolarization) = polarization_vector(polarization(pol), reference_momentum(field))
 
 """
 

src/generic/phase_integrals.jl

@@ -32,8 +32,6 @@ function phase_integrals(
     )
 end
 
-
-# FIXME: This is wrong! We need to insert the xi-dependent terms
 function phase_integrals(
         field::AbstractPlaneWaveField{P},
         internal_integral_method::AbstractIntegrationMethod,
@@ -46,7 +44,7 @@ function phase_integrals(
         P <: AbstractIndefinitePolarization,
     }
 
-    dom = domain(field)
+    dom = compact_domain(field)
     max_amp = maximum_amplitude(field)
 
     # TODO: make max_amp factors global

src/generic/volkov_phase.jl

@@ -35,7 +35,6 @@ function volkov_phase(
     return res
 end
 
-# WARN: Do not forget the xi-dependence!
 function _volkov_phase(
         field::AbstractBackgroundField{AbstractIndefinitePolarization},
         method::AbstractIntegrationMethod,
@@ -44,6 +43,11 @@ function _volkov_phase(
         beta12::Number,
         beta2::Number
     )
-    # add volkov phase for indefinite polarization
+    max_amp = maximum_amplitude(field)
+
+    ii = _internal_integrals(field, method, phi)
+
+    # "-" comes from eps_BG*eps_BG
+    return max_amp * (beta11 * ii.I11.value + beta12 * ii.I12.value) - max_amp^2 * beta2 * ii.I2.value
 
 end

src/interface.jl

@@ -58,6 +58,8 @@ Interface function for background fields. Returns the value of the amplitude at
 function _amplitude end
 
 #TODO: consider moving to generics
+# - OR consider making this the actual interface function, because the three-args version
+# is never called.
 @inline function _amplitude(field::AbstractPlaneWaveField{P}, phi::Real) where {P <: AbstractDefinitePolarization}
     return _amplitude(field, polarization(field), phi)
 end

src/kinematic_factors.jl

@@ -1,51 +1,127 @@
-# kinematic factors
+# TODO: implement kin factors for indefinite polarization (see offline note)
 
+# kinematic factors
 """
-    kinematic_factor1(field::AbstractBackground, p, p_prime) -> Float64
+    kinematic_factor1(field::AbstractBackgroundField, pol, p, p_prime) -> Real
+    kinematic_factor1(field::AbstractBackgroundField, p, p_prime) -> Real
+
+Return the first order kinematic factor
+\$\\beta_1(p, p' \\mid k, \\varepsilon)\$ for an electron interacting with a
+background electromagnetic field.
 
-Compute the kinematic factor \$\\beta_1(p, p'|k, a)\$ for strong-field QED processes.
+The function supports both *definite* and *indefinite* polarization models via
+multiple dispatch.
 
 # Arguments
-- `field::AbstractBackground`: Background electromagnetic field configuration containing
-  the reference momentum \$K\$ and polarization vector \$\\varepsilon^\\mu\$
-- `p`: Initial four-momentum of the particle
-- `p_prime`: Final four-momentum of the particle
+
+- `field::AbstractBackgroundField`: Background electromagnetic field configuration,
+providing the reference momentum \$k\$ and polarization data.
+- `pol::AbstractDefinitePolarization` (optional): Definite polarization state \$\\varepsilon\$.
+- `p`: Incoming electron four-momentum.
+- `p_prime`: Outgoing electron four-momentum.
 
 # Returns
-  - Kinematic factor: \$e((p'\\varepsilon)/(p'K) - (p\\varepsilon)/(pK))\$ where \$e\$ is the elementary charge
 
-# See also
-- [`kinematic_factor2`](@ref): The quadratic kinematic factor \$\\beta_2\$
+- `Real`: The kinematic factor
+    ```math
+    \\begin{align*}
+
+    \\beta_1(p, p' \\mid k, \\varepsilon) =  e\\left(
+      \\frac{p' \\varepsilon}{p' k}
+      -
+      \\frac{p \\varepsilon}{p k}
+    \\right),
+    \\end{align*}
+  ```
+  where \$e\$ is the elementary charge, \$\\varepsilon\$ is the polarization
+  vector, and \$k\$ is the reference momentum of the background field. For
+  mismatched definite polarizations, the return value is zero.
+
+# Calling modes
+
+## Explicit polarization
+
+```julia
+kinematic_factor1(field, pol, p, p_prime)
+```
+If the background field carries a definite polarization, the kinematic factor is evaluated
+only when the explicitly provided polarization `pol` matches the polarization type
+`field`. Otherwise, the function returns zero. In contrast, if the field carries an
+indefinite polarization, the kinematic factor is always evaluated for the given definite
+polarization pol.
+
+## Implicite polarization
+
+```Julia
+kinematic_factor1(field, p, p_prime)
+```
+
+This in only avaliable for background fields with definite polarization, such that the
+polarization of `field` is used implicitly.
+
+- [`kinematic_factor2`](@ref): The second order kinematic factor \$\\beta_1\$
 """
-function kinematic_factor1(field::AbstractBackground, p, p_prime)
-    K = reference_momentum(field)
-    Eps = polarization_vector(field)
+kinematic_factor1
+
+### definite polarization
+# return zero, if the polarization passed in does not match the field polarization
+kinematic_factor1(field::AbstractBackgroundField{<:AbstractDefinitePolarization}, pol::AbstractDefinitePolarization, p, p_prime) = zero(eltype(p))
+
+# if given pol and field pol match, proceed
+kinematic_factor1(field::AbstractBackgroundField{P}, pol::P, p, p_prime) where {P <: AbstractDefinitePolarization} = _kinematic_factor1(field, pol, p, p_prime)
+
+# default: take field pol
+kinematic_factor1(field::AbstractBackgroundField{P}, p, p_prime) where {P <: AbstractDefinitePolarization} = kinematic_factor1(field, polarization(field), p, p_prime)
 
-    term1 = (p_prime * Eps) / (p_prime * K)
-    term2 = (p * Eps) / (p * K)
+### indefinite polarization
+
+kinematic_factor1(field::AbstractBackgroundField{P}, pol::AbstractDefinitePolarization, p, p_prime) where {P <: AbstractIndefinitePolarization} = _kinematic_factor1(field, pol, p, p_prime)
+
+# generic implementation
+function _kinematic_factor1(field::AbstractBackgroundField, pol::AbstractDefinitePolarization, p, p_prime)
+    k = reference_momentum(field)
+    Eps = polarization_vector(field, pol)
+
+    term1 = (p_prime * Eps) / (p_prime * k)
+    term2 = (p * Eps) / (p * k)
     return ELEMENTARY_CHARGE * (term1 - term2)
 end
 
 """
-    kinematic_factor2(field::AbstractBackground, p, p_prime) -> Float64
+    kinematic_factor2(field::AbstractBackgroundField, p, p_prime) -> Float64
 
-Compute the kinematic factor \$\\beta_2(p, p'|k, a)\$ for strong-field QED processes.
+Return the sectond order kinematic factor
+\$\\beta_2(p, p' \\mid k, \\varepsilon)\$ for an electron interacting with a
+background electromagnetic field.
 
 # Arguments
-- `field::AbstractBackground`: Background electromagnetic field configuration containing
-  the reference momentum \$K\$ and polarization vector \$\\varepsilon^\\mu\$
+
+- `field::AbstractBackgroundField`: Background electromagnetic field configuration containing
+  the reference momentum \$k\$ and polarization vector \$\\varepsilon^\\mu\$
 - `p`: Initial four-momentum of the particle
 - `p_prime`: Final four-momentum of the particle
 
 # Returns
-  - Kinematic factor: \$-e^2(1/(p'K) - 1/(pK))\$ where \$e\$ is the elementary charge
+
+- `Real`: The kinematic factor
+  ```math
+    \\begin{align*}
+
+    \\beta_2(p, p' \\mid k, \\varepsilon) =  -e^2\\left(
+      \\frac{1}{p' k}
+      -
+      \\frac{1}{p k}
+    \\right),
+    \\end{align*}
+  ```
+   where \$e\$ is the elementary charge
 
 # See also
-- [`kinematic_factor1`](@ref): The quadratic kinematic factor \$\\beta_1\$
+- [`kinematic_factor1`](@ref): The first order kinematic factor \$\\beta_1\$
 """
-function kinematic_factor2(field::AbstractBackground, p, p_prime)
-    K = reference_momentum(field)
-    term1 = inv(p_prime * K)
-    term2 = inv(p * K)
+function kinematic_factor2(field::AbstractBackgroundField, p, p_prime)
+    k = reference_momentum(field)
+    term1 = inv(p_prime * k)
+    term2 = inv(p * k)
     return ELEMENTARY_CHARGE^2 * (term2 - term1)
 end

src/polarization.jl

@@ -0,0 +1,50 @@
+struct EllipticPolarization{T} <: AbstractIndefinitePolarization
+    xi::T
+end
+
+
+"""
+
+    polarization_vector(pol::AbstractPolarization, mom::QEDbase.AbstractFourMomentum)
+
+Return the polarization vector for a given polarization and four-momentum `mom`.
+For a definite polarization, the respective `LorentzVector` is returned,
+where as for an indefinite polarization, a tuple of polarization vectors is returned.
+
+!!! note "Convention"
+
+    In the current implementation, we use the `base_state` function for `Photon` provided by `QEDcore.jl`.
+
+"""
+@inline function polarization_vector(pol::AbstractDefinitePolarization, mom)
+    return base_state(Photon(), Incoming(), mom, pol)
+end
+
+"""
+
+    oscillator(pol::AbstractPolarizaion, phi::Real)
+
+Return the value of the base oscillator associated with a given polarization `pol` at a given point `phi`.
+
+!!! note "Convention"
+
+    The current default implementation are
+
+    ```Julia
+    PolX() -> cos(phi)
+    PolY() -> sin(phi)
+    AllPol() -> (cos(phi), sin(phi))
+    ```
+
+"""
+function oscillator end
+
+@inline oscillator(::PolX, phi) = cos(phi)
+@inline oscillator(::PolY, phi) = sin(phi)
+
+# TODO: remove this
+@inline function oscillator(::AbstractIndefinitePolarization, phi)
+    sincos_res = sincos(phi)
+    @inbounds cossin_res = (sincos_res[2], sincos_res[1])
+    return cossin_res
+end

src/pulse_profiles/cos_square/internal_integrals.jl

@@ -1,12 +1,14 @@
 @inline _sinc(x) = sinc(x / pi)
 
+# solved integral_0^phi dphi' g(phi') * cos(phi)
 @inline function _internal_integral1_cos_square(::PolX, phi, dphi)
 
     fac_plus = pi / dphi + 1
     fac_minus = pi / dphi - 1
     return sin(phi) / 2 + phi / 4 * (_sinc(fac_plus * phi) + _sinc(fac_minus * phi))
 end
 
+# solved integral_0^phi dphi' g(phi') * sin(phi)
 @inline function _internal_integral1_cos_square(::PolY, phi, dphi)
 
     fac_plus = 1 + pi / dphi
@@ -17,6 +19,7 @@ end
     )
 end
 
+# solved integral_0^phi dphi' (g(phi') * cos(phi))^2
 @inline function _internal_integral2_cos_square(::PolX, phi, dphi)
     k0 = pi / dphi
     k1p = 1 + k0
@@ -35,6 +38,7 @@ end
     return res / 16
 end
 
+# solved integral_0^phi dphi' (g(phi') * sin(phi))^2
 @inline function _internal_integral2_cos_square(::PolY, phi, dphi)
     k0 = pi / dphi
     k1p = 1 + k0
@@ -53,10 +57,9 @@ end
     return res / 16
 end
 
-# FIXME: This is wrong! We need to insert the xi-dependent terms
-@inline function _internal_integral2_cos_square(phi, dphi)
-    value_I21 = _internal_integral2_cos_square(pol, phi, pulse_len)
-    value_I22 = _internal_integral2_cos_square(pol, phi, pulse_len)
+@inline function _internal_integral2_cos_square_indef(phi, dphi, xi)
+    value_I21 = oscillator(PolX(), pol.xi)^2 * _internal_integral2_cos_square(pol, phi, pulse_len)
+    value_I22 = oscillator(PolY(), pol.xi)^2 * _internal_integral2_cos_square(pol, phi, pulse_len)
 
     return value_I21 + value_I22
 end
@@ -75,14 +78,12 @@ end
     )
 end
 
-# FIXME: This is wrong! We need to insert the xi-dependent terms
-@inline function _internal_integrals(pulse::CosSquarePulse, pol::P, method::Analytical, phi::T) where {T <: Number, P <: AbstractIndefinitePolarization}
+@inline function _internal_integrals(pulse::CosSquarePulse, pol::EllipticPolarization, method::Analytical, phi::T) where {T <: Number}
     pulse_len = pulse_length(pulse)
-    value_I11 = _internal_integral1_cos_square(PolX(), phi, pulse_len)
-    value_I12 = _internal_integral2_cos_square(PolY(), phi, pulse_len)
+    value_I11 = oscillator(PolX(), pol.xi) * _internal_integral1_cos_square(PolX(), phi, pulse_len)
+    value_I12 = oscillator(PolY(), pol.xi) * _internal_integral2_cos_square(PolY(), phi, pulse_len)
 
-    # TODO: update with xi!
-    value_I2 = _internal_integral2_cos_square(phi, pulse_len)
+    value_I2 = _internal_integral2_cos_square_indef(phi, pulse_len, pol.xi)
 
     return InternalIntegrals(res12, res12, res2)
 end

src/pulse_profiles/generic.jl

@@ -16,13 +16,16 @@ function _internal_integrals(pulse::AbstractPulseProfile, pol::P, method::Abstra
     )
 end
 
-# FIXME: This is wrong! We need to insert the xi-dependent terms
 @inline function _internal_integrals(pulse::AbstractPulseProfile, pol::P, method::AbstractNumericalIntegrationMethod, phi::T) where {T <: Number, P <: AbstractIndefinitePolarization}
-    tmp_func11 = t -> oscillator(PolX(), phi) * _envelope(pulse, phi)
-    tmp_func12 = t -> oscillator(PolY(), phi) * _envelope(pulse, phi)
+    pol_fac_X = oscillator(PolX(), pol.xi)
+    pol_fac_Y = oscillator(PolY(), pol.xi)
 
-    # FIXME: this one misses xi
-    tmp_func2 = t -> _envelope(pulse, t)^2 * (oscillator(PolX(), t)^2 + oscillator(PolY(), t)^2)
+    tmp_func11 = t -> pol_fac_X * oscillator(PolX(), phi) * _envelope(pulse, phi)
+    tmp_func12 = t -> pol_fac_Y * oscillator(PolY(), phi) * _envelope(pulse, phi)
+
+    tmp_func2 = t -> _envelope(pulse, t)^2 * (
+        (pol_fac_X * oscillator(PolX(), t))^2 + (pol_fac_Y * oscillator(PolY(), t))^2
+    )
 
     res11 = integrate(method, tmp_func11, zero(T), phi)
     res12 = integrate(method, tmp_func12, zero(T), phi)