Refac Field Infrastructure

Project.toml

@@ -10,16 +10,24 @@ authors = [
 version = "0.3.1"
 
 [deps]
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+GaussQuadrature = "d54b0c1a-921d-58e0-8e36-89d8069c0969"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
 QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
+FastGaussQuadrature = "1.0.2"
+GaussQuadrature = "0.5.8"
+IntervalSets = "0.7"
+LinearAlgebra = "1.10.0"
 QEDbase = "0.5"
 QEDcore = "0.4"
-IntervalSets = "0.7"
 QuadGK = "2"
+SpecialFunctions = "2.6.1"
 julia = "1.10"
 
 [extras]

src/QEDfields.jl

@@ -1,21 +1,75 @@
+# TODO: implement in-place versions:
+# - amplitude!
+# - compute!
+# - internal_integrals!
+# - volkov_phase!
+# - integrate!
+#
 module QEDfields
 
-using QEDbase
-using QEDcore
+# result types
+export InternalIntegralResult, InternalIntegrals, PhaseIntegralResult, PhaseIntegrals
+
+# background fields
+export AbstractBackgroundField
+
+# plane-wave fields
+export AbstractPlaneWaveField, polarization_type
+export oscillator
+export amplitude, reference_momentum, polarization, classical_nonlinearity_parameter
+export a0_parameter, internal_integrals, volkov_phase, phase_integrals
+export maximum_amplitude
+
+# pulsed fields
+export AbstractPulsedPlaneWaveField
+export pulse_profile
+export PulsedPlaneWaveField
+
+# pulse profiles
+export AbstractPulseProfile
+export domain, compact_domain, pulse_length, envelope
+export CosSquarePulse
+export GaussianPulse
+
+# integration methods
+export integrate
+export GaussKronrodQuadrature
+export GaussLegendreQuadrature
+export Analytical
 
-using IntervalSets
 using QuadGK
+using FastGaussQuadrature
+using GaussQuadrature
+using IntervalSets
+using LinearAlgebra
+using SpecialFunctions
+
+
+using QEDcore
+using QEDbase
+
+include("patch_QEDcore.jl")
+include("results.jl")
 
-export AbstractBackgroundField, AbstractPulsedPlaneWaveField
-export reference_momentum, domain, pulse_length, envelope, amplitude, generic_spectrum
+include("integration_methods/interface.jl")
+include("integration_methods/gauss_kronrod.jl")
+include("integration_methods/gauss_legendre.jl")
+include("integration_methods/analytical.jl")
 
-export polarization_vector, oscillator
+include("interface.jl")
+include("generic/amplitude.jl")
+include("generic/internal_integrals.jl")
+include("generic/volkov_phase.jl")
+include("generic/phase_integrals.jl")
 
-export CosSquarePulse, GaussianPulse
+include("pulse_profiles/interface.jl")
+include("pulse_profiles/generic.jl")
+include("pulse_profiles/cos_square/impl.jl")
+include("pulse_profiles/cos_square/internal_integrals.jl")
+include("pulse_profiles/gaussian_pulse/impl.jl")
 
-include("interfaces/background_field_interface.jl")
-include("polarization.jl")
-include("pulses/cos_square.jl")
-include("pulses/gaussian.jl")
+include("pulsed_fields/interface.jl")
+include("pulsed_fields/generic.jl")
+include("pulsed_fields/impl.jl")
 
 end

src/generic/amplitude.jl

@@ -0,0 +1,29 @@
+### Generic implementations for plane-wave fields
+
+# TODO: move oscillator to pulsed plane wave field. (or oscillating fields as a subtype?)
+# delegations
+oscillator(field::AbstractPlaneWaveField, phi::Real) = oscillator(polarization(field), phi)
+polarization_vector(field::AbstractPlaneWaveField) = polarization_vector(polarization(field), reference_momentum(field))
+
+"""
+
+    amplitude(field::AbstractPlaneWaveField, pol::AbstractDefinitePolarization, phi)
+
+Returns the value of the amplitude for a given polarization direction and phase variable `phi`.
+
+!!! 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.
+"""
+function amplitude(
+        field::AbstractPlaneWaveField, pol::AbstractDefinitePolarization, phi::Real
+    )
+    return phi in domain(field) ? _amplitude(field, pol, phi) : zero(phi)
+end
+
+function amplitude(
+        field::AbstractPlaneWaveField{P}, phi::Real,
+    ) where {P <: AbstractDefinitePolarization}
+    return phi in domain(field) ? _amplitude(field, phi) : zero(phi)
+end

src/generic/internal_integrals.jl

@@ -0,0 +1,31 @@
+### internal integrals
+
+@inline function _internal_integrals(field::AbstractPlaneWaveField{P}, method::AbstractIntegrationMethod, phi::T) where {T <: Real, P <: AbstractDefinitePolarization}
+    tmp_func1 = t -> _amplitude(field, t)
+    tmp_func2 = t -> _amplitude(field, t)^2
+
+    res1 = integrate(method, tmp_func1, zero(T), phi)
+    res2 = integrate(method, tmp_func2, zero(T), phi)
+    return InternalIntegrals(res1, res2)
+end
+
+# FIXME: This is wrong! We need to insert the xi-dependent terms
+@inline function _internal_integrals(field::AbstractPlaneWaveField{P}, method::AbstractIntegrationMethod, phi::T) where {T <: Real, P <: AbstractIndefinitePolarization}
+    tmp_func11 = t -> _amplitude(field, t, PolX())
+    tmp_func12 = t -> _amplitude(field, t, PolY())
+    tmp_func2 = t -> _amplitude(field, t, PolX())^2 + _amplitude(field, t, PolY())^2
+
+    res11 = integrate(method, tmp_func11, zero(T), phi)
+    res12 = integrate(method, tmp_func12, zero(T), phi)
+    res2 = integrate(method, tmp_func2, zero(T), phi)
+    return InternalIntegrals(res12, res12, res2)
+end
+
+function internal_integrals(field::AbstractPlaneWaveField{P}, method::AbstractIntegrationMethod, phi::T)::InternalIntegrals{T, 2} where {T <: Real, P <: AbstractDefinitePolarization}
+
+    dom = compact_domain(field)
+
+    phi_eval = phi <= infimum(dom) ? infimum(dom) : min(phi, supremum(dom))
+    return _internal_integrals(field, method, phi_eval)
+
+end

src/generic/phase_integrals.jl

@@ -0,0 +1,65 @@
+# TBW
+
+_tuple_pack(x::T) where {T <: Number} = (x,)
+_tuple_pack(x::Tuple) = x
+
+function phase_integrals(
+        field::AbstractPlaneWaveField{P},
+        internal_integral_method::AbstractIntegrationMethod,
+        phase_integral_method::AbstractIntegrationMethod,
+        pnum::T,
+        beta1::T,
+        beta2::T
+    )::PhaseIntegrals{Complex{T}, 2} where {
+        T <: Real,
+        P <: AbstractDefinitePolarization,
+    }
+
+    dom = domain(field)
+    max_amp = maximum_amplitude(field)
+
+    # TODO: make max_amp factors global
+    integrand1 = x -> max_amp * _amplitude(field, x) * exp(1im * pnum * x + 1im * volkov_phase(field, internal_integral_method, x, beta1, beta2))
+    integrand2 = x -> max_amp^2 * _amplitude(field, x)^2 * exp(1im * pnum * x + 1im * volkov_phase(field, internal_integral_method, x, beta1, beta2))
+
+    res1 = integrate(phase_integral_method, integrand1, dom)
+    res2 = integrate(phase_integral_method, integrand2, dom)
+
+    return PhaseIntegrals(
+        PhaseIntegralResult(_tuple_pack(res1)...),
+        PhaseIntegralResult(_tuple_pack(res2)...),
+    )
+end
+
+
+# FIXME: This is wrong! We need to insert the xi-dependent terms
+function phase_integrals(
+        field::AbstractPlaneWaveField{P},
+        internal_integral_method::AbstractIntegrationMethod,
+        phase_integral_method::AbstractIntegrationMethod,
+        pnum::Real,
+        beta11::Real,
+        beta12::Real,
+        beta2::Real
+    ) where {
+        P <: AbstractIndefinitePolarization,
+    }
+
+    dom = domain(field)
+    max_amp = maximum_amplitude(field)
+
+    # TODO: make max_amp factors global
+    integrand11 = x -> max_amp * _amplitude(field, PolX(), x) * exp(1im * pnum * x + 1im * volkov_phase(field, x, beta11, beta12, beta2))
+    integrand12 = x -> max_amp * _amplitude(field, PolY(), x) * exp(1im * pnum * x + 1im * volkov_phase(field, x, beta11, beta12, beta2))
+    integrand2 = x -> max_amp^2 * (_amplitude(field, PolX(), x)^2 + _amplitude(field, PolY(), x)^2) * exp(1im * pnum * x + 1im * volkov_phase(field, x, beta11, beta12, beta2))
+
+    res11 = integrate(method, integrand11, dom)
+    res12 = integrate(method, integrand12, dom)
+    res2 = integrate(method, integrand2, dom)
+
+    return PhaseIntegrals(
+        PhaseIntegralResult(_tuple_pack(res11)...),
+        PhaseIntegralResult(_tuple_pack(res12)...),
+        PhaseIntegralResult(_tuple_pack(res2)...),
+    )
+end

src/generic/volkov_phase.jl

@@ -0,0 +1,49 @@
+### non-linear Volkov phase
+
+function _volkov_phase(
+        field::AbstractBackgroundField{<:AbstractDefinitePolarization},
+        method::AbstractIntegrationMethod,
+        phi::Real,
+        beta1::Real,
+        beta2::Real
+    )
+
+    max_amp = maximum_amplitude(field)
+
+    ii = _internal_integrals(field, method, phi)
+
+    # "-" comes from eps_BG*eps_BG
+    return max_amp * beta1 * ii.I1.value - max_amp^2 * beta2 * ii.I2.value
+end
+
+function volkov_phase(
+        field::AbstractPlaneWaveField,
+        method::AbstractIntegrationMethod,
+        phi::T,
+        beta1::T,
+        beta2::T
+    )::T where {T <: Real}
+
+    dom = compact_domain(field)
+
+    if phi in dom
+        res = _volkov_phase(field, method, phi, beta1, beta2)
+    else
+        res = phi <= minimum(dom) ? _volkov_phase(field, method, minimum(dom), beta1, beta2) : _volkov_phase(field, method, maximum(dom), beta1, beta2)
+    end
+
+    return res
+end
+
+# WARN: Do not forget the xi-dependence!
+function _volkov_phase(
+        field::AbstractBackgroundField{AbstractIndefinitePolarization},
+        method::AbstractIntegrationMethod,
+        phi::Real,
+        beta11::Real,
+        beta12::Real,
+        beta2::Real
+    )
+    # add volkov phase for indefinite polarization
+
+end

src/integration_methods/analytical.jl

@@ -0,0 +1,7 @@
+"""
+
+    Analytical()
+
+General analytical solution of an integral.
+"""
+struct Analytical <: AbstractAnalyticalMethod end

src/integration_methods/gauss_kronrod.jl

@@ -0,0 +1,16 @@
+# TODO:
+# - pass kwargs to GaussKronrodQuadrature type
+# - consider extenting the interface for quadrature rules (e.g. order(::AbstractQuadratureMethod),...)
+
+"""
+
+    GaussKronrodQuadrature()
+
+Integration method using `QuadGK.quadgk`.
+"""
+struct GaussKronrodQuadrature <: AbstractQuadratureMethod end
+
+function integrate(meth::GaussKronrodQuadrature, func::Function, low::Real, high::Real)
+    res, _ = quadgk(func, low, 0, high)
+    return res
+end

src/integration_methods/gauss_legendre.jl

@@ -0,0 +1,44 @@
+# TODO:
+# - consider extenting the interface for quadrature rules (e.g. order(::AbstractQuadratureMethod),...)
+
+"""
+
+    GaussLegendreQuadrature(order)
+
+Integration method using Gauss-Legendre nodes provided by `FastGaussQuadrature.jl`
+"""
+struct GaussLegendreQuadrature{P, W} <: AbstractQuadratureMethod
+    order::Int
+    nodes::P
+    weights::W
+
+    function GaussLegendreQuadrature(order::Int; dtype = nothing)
+
+        x, w = gausslegendre(order)
+        #x, w = legendre(order)
+        if dtype != nothing
+            x = convert(Vector{dtype}, x)
+            w = convert(Vector{dtype}, w)
+        end
+        return new{typeof(x), typeof(w)}(order, x, w)
+    end
+end
+function Base.show(io::IO, method::GaussLegendreQuadrature)
+    return print(io, "GaussLegendreQuadrature(order = $(method.order))")
+end
+function quadrature_nodes(method::GaussLegendreQuadrature, low, high)
+    slope = (high - low) / 2
+    shift = (low + high) / 2
+    return @. slope * method.nodes + shift
+end
+
+function quadrature_weights(method::GaussLegendreQuadrature, low, high)
+    fac = (high - low) / 2
+    return @. fac * method.weights
+end
+
+function integrate(meth::GaussLegendreQuadrature, func::Function, low::Real, high::Real)
+    n = quadrature_nodes(meth, low, high)
+    w = quadrature_weights(meth, low, high)
+    return LinearAlgebra.dot(w, func.(n))
+end

src/integration_methods/interface.jl

@@ -0,0 +1,22 @@
+# TODO: consider using `Integrals.jl` here
+
+# integration methods
+
+abstract type AbstractIntegrationMethod end
+Base.broadcastable(method::AbstractIntegrationMethod) = Ref(method)
+abstract type AbstractAnalyticalMethod <: AbstractIntegrationMethod end
+abstract type AbstractNumericalIntegrationMethod <: AbstractIntegrationMethod end
+
+"""
+
+    integrate(::AbstractNumericalIntegrationMethod,func,low,up)
+
+Interface function for subtypes of `AbstractNumericalIntegrationMethod`.
+"""
+function integrate end
+
+@inline function integrate(method::AbstractNumericalIntegrationMethod, func::Function, domain::AbstractInterval)
+    return integrate(method, func, endpoints(domain)...)
+end
+
+abstract type AbstractQuadratureMethod <: AbstractNumericalIntegrationMethod end