From ad72264843076937a2081cf7d6d336d2132077ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 30 Jan 2026 20:26:53 +0200 Subject: [PATCH 1/9] update to MTK@v11 --- Project.toml | 6 ++--- src/ElectronDynamicsModels.jl | 3 --- src/base.jl | 15 ++++++++++++ src/external_fields.jl | 46 +++++++++++++++++++++++++++++------ 4 files changed, 57 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 90ebdc5..4fe24d5 100644 --- a/Project.toml +++ b/Project.toml @@ -15,10 +15,10 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [compat] Aqua = "0.8" HypergeometricFunctions = "0.3.28" -JET = "0.10" +JET = "0.10, 0.11" LaserTypes = "0.2" LinearAlgebra = "1.11.0" -ModelingToolkit = "10.10" +ModelingToolkit = "10.10, 11" OrdinaryDiffEqNonlinearSolve = "1.10.0" OrdinaryDiffEqRosenbrock = "1.11.0" OrdinaryDiffEqVerner = "1.2.0" @@ -26,7 +26,7 @@ PhysicalConstants = "0.2.3" Plots = "1" SciMLBase = "2.102.1" Statistics = "1" -Symbolics = "6.43.0" +Symbolics = "6.43.0, 7" Test = "1.11" Unitful = "1.21.0" UnitfulAtomic = "1.0.0" diff --git a/src/ElectronDynamicsModels.jl b/src/ElectronDynamicsModels.jl index 3b04412..8c922fe 100644 --- a/src/ElectronDynamicsModels.jl +++ b/src/ElectronDynamicsModels.jl @@ -7,9 +7,6 @@ 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 diff --git a/src/base.jl b/src/base.jl index 7530f3a..870e16f 100644 --- a/src/base.jl +++ b/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 diff --git a/src/external_fields.jl b/src/external_fields.jl index 575574e..5806d02 100644 --- a/src/external_fields.jl +++ b/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 ] - sys = System(eqs, iv, [x, t, wz, z, r], [λ, a₀, ω, k, E₀, w₀, z_R, T0, τ0]; name, systems=[ref_frame]) + 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], + initial_conditions, + initialization_eqs, + bindings + ) extend(sys, field_dynamics) end @@ -145,9 +163,14 @@ Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970) B[1] ~ 0 B[2] ~ A/c * cos(dot(k, x⃗) - ω * t) B[3] ~ 0 - # parameters + ] + initialization_eqs = [ λ ~ (2π * c) / ω ] + bindings = [ + λ => missing + ω => missing + ] sys = System(eqs, iv, [x, t, x⃗], [A, ω, k, λ]; name, systems=[ref_frame]) @@ -376,13 +399,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 +422,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 From 60e975edcbaa0e425295c3eeaeab556fd9638493 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 2 Feb 2026 18:35:52 +0200 Subject: [PATCH 2/9] fix PlaneWave implementation to take k into account correctly Co-Authored-By: Claude Opus 4.5 --- src/external_fields.jl | 54 ++++++++++++++++++++++++++++++------------ 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/external_fields.jl b/src/external_fields.jl index 5806d02..fb4b7a8 100644 --- a/src/external_fields.jl +++ b/src/external_fields.jl @@ -136,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) @@ -145,25 +145,47 @@ 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 + 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) / ω ] @@ -172,7 +194,9 @@ Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970) ω => missing ] - sys = System(eqs, iv, [x, t, x⃗], [A, ω, k, λ]; name, systems=[ref_frame]) + vars = nameof(iv) == :τ ? [x, t] : [x] + + sys = System(eqs, iv, vars, [A, ω, k_dir, pol, λ]; name, systems=[ref_frame], initialization_eqs, bindings) extend(sys, field_dynamics) end From 7a5fd34e2e391b5ed769f96cc1802afc965f386d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 2 Feb 2026 18:36:13 +0200 Subject: [PATCH 3/9] fix JET testing --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 715469f..85789ec 100644 --- a/test/runtests.jl +++ b/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") From 02159efd384411432c7175961dba87fd768526e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 2 Feb 2026 18:36:36 +0200 Subject: [PATCH 4/9] update PlaneWave test --- test/test_classical_electron.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_classical_electron.jl b/test/test_classical_electron.jl index 59165a3..cccb78f 100644 --- a/test/test_classical_electron.jl +++ b/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) From 438f4699b53fa2d6bedab4d2d71e847d79401ac2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 2 Feb 2026 18:36:49 +0200 Subject: [PATCH 5/9] sort imports --- src/ElectronDynamicsModels.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ElectronDynamicsModels.jl b/src/ElectronDynamicsModels.jl index 8c922fe..07421bb 100644 --- a/src/ElectronDynamicsModels.jl +++ b/src/ElectronDynamicsModels.jl @@ -1,8 +1,8 @@ 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 From 180c3ee1d0f1eeadce3e2fb6f00073ff490e3ce8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 2 Feb 2026 18:37:06 +0200 Subject: [PATCH 6/9] fix parameter propagation in AbrahamLorentzRadiation --- src/radiation_reaction.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/radiation_reaction.jl b/src/radiation_reaction.jl index 597fe7e..27c18a6 100644 --- a/src/radiation_reaction.jl +++ b/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) From 7d3e8b85239ed95b0e5ee41abdedd9bee02a67ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Mon, 2 Feb 2026 19:07:53 +0200 Subject: [PATCH 7/9] =?UTF-8?q?fix=20=CE=B5=E2=82=80=20definition?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Claude --- test/test_radiation.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/test_radiation.jl b/test/test_radiation.jl index c166d02..aa93bc9 100644 --- a/test/test_radiation.jl +++ b/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 From 6cd0a04f8ecff5b5184c41bfe59d3c7560b4e31b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 8 Feb 2026 22:41:01 +0200 Subject: [PATCH 8/9] update compats --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 4fe24d5..657b837 100644 --- a/Project.toml +++ b/Project.toml @@ -15,10 +15,10 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [compat] Aqua = "0.8" HypergeometricFunctions = "0.3.28" -JET = "0.10, 0.11" +JET = "0.9, 0.10, 0.11" LaserTypes = "0.2" LinearAlgebra = "1.11.0" -ModelingToolkit = "10.10, 11" +ModelingToolkit = "11" OrdinaryDiffEqNonlinearSolve = "1.10.0" OrdinaryDiffEqRosenbrock = "1.11.0" OrdinaryDiffEqVerner = "1.2.0" @@ -26,7 +26,7 @@ PhysicalConstants = "0.2.3" Plots = "1" SciMLBase = "2.102.1" Statistics = "1" -Symbolics = "6.43.0, 7" +Symbolics = "7.10" Test = "1.11" Unitful = "1.21.0" UnitfulAtomic = "1.0.0" From 7ad112e18d0f452acf882efa75ac975895a05764 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Sun, 8 Feb 2026 22:50:47 +0200 Subject: [PATCH 9/9] run ci on 1.12 instead of pre JET not compatible with pre --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 6f534cf..0dee6a2 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,7 +24,7 @@ jobs: matrix: version: - '1.11' - - 'pre' + - '1' os: - ubuntu-latest arch: