update to MTK@v11

.github/workflows/CI.yml

@@ -24,7 +24,7 @@ jobs:
       matrix:
         version:
           - '1.11'
-          - 'pre'
+          - '1'
         os:
           - ubuntu-latest
         arch:

Project.toml

@@ -15,18 +15,18 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 [compat]
 Aqua = "0.8"
 HypergeometricFunctions = "0.3.28"
-JET = "0.10"
+JET = "0.9, 0.10, 0.11"
 LaserTypes = "0.2"
 LinearAlgebra = "1.11.0"
-ModelingToolkit = "10.10"
+ModelingToolkit = "11"
 OrdinaryDiffEqNonlinearSolve = "1.10.0"
 OrdinaryDiffEqRosenbrock = "1.11.0"
 OrdinaryDiffEqVerner = "1.2.0"
 PhysicalConstants = "0.2.3"
 Plots = "1"
 SciMLBase = "2.102.1"
 Statistics = "1"
-Symbolics = "6.43.0"
+Symbolics = "7.10"
 Test = "1.11"
 Unitful = "1.21.0"
 UnitfulAtomic = "1.0.0"

src/ElectronDynamicsModels.jl

@@ -1,15 +1,12 @@
 module ElectronDynamicsModels
 
 using ModelingToolkit
-using Unitful, UnitfulAtomic, PhysicalConstants
-using PhysicalConstants.CODATA2018: e, m_e, c_0, ε_0
+using PhysicalConstants, Unitful, UnitfulAtomic
+using PhysicalConstants.CODATA2018: c_0, e, m_e, ε_0
 using LinearAlgebra
 using Symbolics
 using HypergeometricFunctions: HypergeometricFunctions, _₁F₁, pochhammer
 
-# Register hypergeometric function for symbolic use
-@register_symbolic HypergeometricFunctions._₁F₁(a, b, z)
-
 m_dot(x, y) = x[1] * y[1] - x[2] * y[2] - x[3] * y[3] - x[4] * y[4]
 
 export GaussLaser, LaguerreGaussLaser

src/base.jl

@@ -4,6 +4,16 @@ function ReferenceFrame(iv; name, c, ε₀, μ₀, m_e, q_e)
     System(Equation[], iv, [], GlobalScope.([constants..., gμν]); name)
 end
 
+# struct Constants{T <: Real}
+#     c::T
+#     ε₀::T
+#     μ₀::T
+#     m_e::T
+#     q_e::T
+# end
+
+# @symstruct Constants{T}
+
 function ReferenceFrame(iv, units::Symbol; name)
 
     if units == :SI
@@ -41,6 +51,11 @@ function ReferenceFrame(iv, units::Symbol; name)
         error("$units not supported!")
     end
 
+    # @parameters gμν[1:4, 1:4] = diagm([1, -1, -1, -1])
+    # @parameters constants = Constants(c, ε₀, μ₀, m_e, q_e)
+
+    # System(Equation[], iv, [], GlobalScope.([constants, gμν]); name)
+
     ReferenceFrame(iv; name, c, ε₀, μ₀, m_e, q_e)
 end
 

src/external_fields.jl

@@ -95,16 +95,34 @@ Parameters:
             exp(im * ω * t) *
             exp(-(((t - t₀) - (z - z₀) / c) / τ0)^2),
         )
+    ]
 
-        # Parameter relations
+    initialization_eqs = [
         ω ~ 2π * c / λ
         E₀ ~ a₀ * m_e * c * ω / abs(q_e)
         z_R ~ w₀^2 * k / 2
         k ~ 2π / λ
-        τ0 ~ 10 / ω
+    ]
+    bindings = [
+        ω => missing
+        λ => missing
+        k => missing
+        E₀ => missing
+        z_R => missing
+        w₀ => missing
+    ]
+
+    initial_conditions = [
+        τ0 => 10 / ω
     ]
 
-    sys = System(eqs, iv, [x, t, wz, z, r], [λ, a₀, ω, k, E₀, w₀, z_R, T0, τ0]; name, systems=[ref_frame])
+    sys = System(eqs, iv, [x, t, wz, z, r], [λ, a₀, ω, k, E₀, w₀, z_R, T0, τ0];
+        name,
+        systems=[ref_frame],
+        initial_conditions,
+        initialization_eqs,
+        bindings
+    )
 
     extend(sys, field_dynamics)
 end
@@ -118,7 +136,7 @@ electrons can exhibit figure-8 motion when a₀ ~ 1.
 
 Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970)
 """
-@component function PlaneWave(; name, amplitude=1.0, frequency=1.0, k_vector=[0,0,1], ref_frame)
+@component function PlaneWave(; name, amplitude=1.0, frequency=1.0, k_direction=[0,0,1], polarization=[1,0,0], ref_frame)
     # New interface with spacetime
     @named field_dynamics = EMFieldDynamics(; ref_frame)
 
@@ -127,29 +145,58 @@ Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970)
     iv = ModelingToolkit.get_iv(ref_frame)
 
     # Create local position and time variables
-    @variables x(iv)[1:4] t(iv)
+    @variables x(iv)[1:4]
+    if nameof(iv) == :τ
+        @variables t(iv)
+    else
+        t = iv
+    end
 
     @unpack E, B = field_dynamics
 
-    @parameters A=amplitude ω=frequency k[1:3]=k_vector λ
-    @variables x⃗(iv)[1:3]
+    @parameters A=amplitude ω=frequency k_dir[1:3]=k_direction pol[1:3]=polarization λ
+
+    # Normalize k direction to get unit vector k̂
+    k_norm = sqrt(k_dir[1]^2 + k_dir[2]^2 + k_dir[3]^2)
+    k̂ = [k_dir[1] / k_norm, k_dir[2] / k_norm, k_dir[3] / k_norm]
+
+    # Normalize polarization vector (user must ensure it's perpendicular to k)
+    pol_norm = sqrt(pol[1]^2 + pol[2]^2 + pol[3]^2)
+    ê = [pol[1] / pol_norm, pol[2] / pol_norm, pol[3] / pol_norm]
+
+    # B-field direction: k̂ × ê
+    b̂ = [
+        k̂[2] * ê[3] - k̂[3] * ê[2],
+        k̂[3] * ê[1] - k̂[1] * ê[3],
+        k̂[1] * ê[2] - k̂[2] * ê[1]
+    ]
+
+    # Spatial position from 4-position (x = [ct, x, y, z] or [t, x, y, z])
+    x⃗ = [x[2], x[3], x[4]]
+
+    # Phase: k·r - ωt where |k| = ω/c (dispersion relation)
+    phase = (ω / c) * (k̂[1] * x⃗[1] + k̂[2] * x⃗[2] + k̂[3] * x⃗[3]) - ω * t
 
     eqs = [
-        # Define spatial position from 4-position
-        x⃗[1] ~ x[2]
-        x⃗[2] ~ x[3]
-        x⃗[3] ~ x[4]
-        E[1] ~ A * cos(dot(k, x⃗) - ω * t)
-        E[2] ~ 0
-        E[3] ~ 0
-        B[1] ~ 0
-        B[2] ~ A/c * cos(dot(k, x⃗) - ω * t)
-        B[3] ~ 0
-        # parameters
+        E[1] ~ A * ê[1] * cos(phase)
+        E[2] ~ A * ê[2] * cos(phase)
+        E[3] ~ A * ê[3] * cos(phase)
+        B[1] ~ (A / c) * b̂[1] * cos(phase)
+        B[2] ~ (A / c) * b̂[2] * cos(phase)
+        B[3] ~ (A / c) * b̂[3] * cos(phase)
+    ]
+
+    initialization_eqs = [
         λ ~ (2π * c) / ω
     ]
+    bindings = [
+        λ => missing
+        ω => missing
+    ]
+
+    vars = nameof(iv) == :τ ? [x, t] : [x]
 
-    sys = System(eqs, iv, [x, t, x⃗], [A, ω, k, λ]; name, systems=[ref_frame])
+    sys = System(eqs, iv, vars, [A, ω, k_dir, pol, λ]; name, systems=[ref_frame], initialization_eqs, bindings)
 
     extend(sys, field_dynamics)
 end
@@ -376,13 +423,21 @@ Reference: Allen et al., Phys. Rev. A 45, 8185 (1992)
                 )
             )
         )
-
-        # Parameter relations
+    ]
+    initialization_eqs = [
         ω ~ 2π * c / λ
         k ~ 2π / λ
         z_R ~ π * w₀^2 / λ
         E₀ ~ a₀ * m_e * c * ω / abs(q_e)
     ]
+    bindings = [
+        ω => missing
+        λ => missing
+        k => missing
+        E₀ => missing
+        z_R => missing
+        w₀ => missing
+    ]
 
     vars = if nameof(iv) == :τ
         [x, t, z, r, θ, wz, σ, rwz, env]
@@ -391,6 +446,7 @@ Reference: Allen et al., Phys. Rev. A 45, 8185 (1992)
         [x, τ, z, r, θ, wz, σ, rwz, env]
     end
 
-    sys = System(eqs, iv, vars, params; name, systems=[ref_frame])
+    sys = System(eqs, iv, vars, params;
+        name, systems=[ref_frame], initialization_eqs, bindings)
     extend(sys, field_dynamics)
 end

src/radiation_reaction.jl

@@ -91,10 +91,10 @@ This is valid when ω << m*c²/ℏ (classical regime).
     iv = ModelingToolkit.get_iv(ref_frame)
     @named field = ElectromagneticSystem(iv)
 
-    @unpack c, gμν = ref_frame
+    @unpack c, gμν, ε₀ = ref_frame
     @unpack u = particle
 
-    @parameters m=1.0 ε₀=1.0 q=charge
+    @parameters m=1.0 q=charge
     @variables begin
         F_rad(iv)[1:4]     # 4-force from radiation reaction
         P_rad(iv)          # Radiated power (invariant)

test/runtests.jl

@@ -8,7 +8,7 @@ using JET
         Aqua.test_all(ElectronDynamicsModels)
     end
     @testset "Code linting (JET.jl)" begin
-        JET.test_package(ElectronDynamicsModels; target_defined_modules = (ElectronDynamicsModels,))
+        JET.test_package(ElectronDynamicsModels; target_modules = (ElectronDynamicsModels,))
     end
     @testset "Classical Electron" begin
         include("test_classical_electron.jl")

test/test_classical_electron.jl

@@ -26,7 +26,7 @@ using LinearAlgebra
         c = 1
         @independent_variables τ
         @named ref_frame = ReferenceFrame(τ; c, ε₀=1, μ₀=1, m_e=1, q_e=1)
-        external_field = PlaneWave(; ref_frame, name = :plane_wave, k_vector = [0, 0, k])
+        external_field = PlaneWave(; ref_frame, name = :plane_wave, k_direction = [0, 0, k])
         @named electron = ChargedParticle(; ref_frame, external_field)
         sys = mtkcompile(electron)
 

test/test_radiation.jl

@@ -255,10 +255,11 @@ using Statistics
             v_max = maximum(sqrt.(sum(u_spatial .^ 2, dims = 2)) ./ u_0)
             @test v_max < 0.1  # v < 0.1c
 
-            # In natural units (q=1, m=1, ε₀=1, c=1):
-            # P_expected = (2/3) * (1/(6π)) * E₀²
-            # where τ₀ = 1/(6π) is the classical electron radius time scale
-            τ₀ = 1.0 / (6π)
+            # In natural units (q=1, m=1, c=1), with ε₀ = 1/(4πα):
+            # τ₀ = q²/(6πε₀mc³) = 2α/3 where α ≈ 1/137
+            # P_expected = (2/3) * τ₀ * E₀²
+            # Get τ₀ directly from the system to ensure consistency
+            τ₀ = sol[sys_nr.radiation.τ₀, end]
             P_expected = (2 / 3) * τ₀ * E₀^2
 
             # Check power at end of simulation