Skip to content

Commit 73da542

Browse files
SebastianM-Cclaude
andcommitted
Add tests for Thomson scattering
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
1 parent d7320fc commit 73da542

4 files changed

Lines changed: 117 additions & 7 deletions

File tree

Project.toml

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
2121
Aqua = "0.8"
2222
CSV = "0.10"
2323
DataInterpolations = "8.9.0"
24+
FFTW = "1"
2425
HypergeometricFunctions = "0.3.28"
2526
JET = "0.9, 0.10, 0.11"
2627
LaserTypes = "0.2"
@@ -29,6 +30,7 @@ ModelingToolkit = "11"
2930
ModelingToolkitBase = "1.20.0"
3031
OrdinaryDiffEqNonlinearSolve = "1.10.0"
3132
OrdinaryDiffEqRosenbrock = "1.11.0"
33+
OrdinaryDiffEqTsit5 = "1.9"
3234
OrdinaryDiffEqVerner = "1.2.0"
3335
PhysicalConstants = "0.2.3"
3436
Plots = "1"
@@ -46,11 +48,13 @@ julia = "1.11"
4648
[extras]
4749
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
4850
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
51+
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
4952
JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
5053
LaserTypes = "e07c0bfa-524c-4f35-a151-c3dd916fa2f0"
5154
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
5255
OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
5356
OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
57+
OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
5458
OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
5559
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
5660
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -61,4 +65,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
6165
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
6266

6367
[targets]
64-
test = ["Aqua", "CSV", "JET", "LaserTypes", "LinearAlgebra", "Test", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqVerner", "Plots", "Random", "SciMLBase", "StaticArrays", "Statistics", "SymbolicIndexingInterface"]
68+
test = ["Aqua", "CSV", "FFTW", "JET", "LaserTypes", "LinearAlgebra", "Test", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Plots", "Random", "SciMLBase", "StaticArrays", "Statistics", "SymbolicIndexingInterface"]

scripts/thomson_scattering.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -128,9 +128,10 @@ screen = ObserverScreen(
128128
x⁰_samples
129129
)
130130

131-
K = sol0.ps[ref_frame.q_e / (4π * ref_frame.ε₀ * c)]
131+
A_s = accumulate_potential(trajs, screen, Tsit5())
132132

133-
A_s = accumulate_potential(trajs, screen, K, alg = Tsit5())
133+
# GPU
134+
# accumulate_potential(trajs, screen, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()))
134135
A_ω = rfft(A_s, 1)
135136

136137
# Find fundamental frequency bin

test/runtests.jl

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -11,15 +11,18 @@ using JET
1111
JET.test_package(ElectronDynamicsModels; target_modules = (ElectronDynamicsModels,))
1212
end
1313
@testset "Classical Electron" begin
14-
include("test_classical_electron.jl")
14+
include("classical_electron.jl")
15+
end
16+
@testset "Radiation Emission" begin
17+
include("thomson_scattering.jl")
1518
end
1619
@testset "Radiation Reaction" begin
17-
include("test_radiation.jl")
20+
include("radiation.jl")
1821
end
1922
@testset "GaussLaser" begin
20-
include("test_gauss_laser.jl")
23+
include("gauss_laser.jl")
2124
end
2225
@testset "LaguerreGauss" begin
23-
include("test_laguerre_gauss.jl")
26+
include("laguerre_gauss.jl")
2427
end
2528
end

test/thomson_scattering.jl

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
using ElectronDynamicsModels
2+
using ModelingToolkit
3+
using OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5
4+
using StaticArrays
5+
using Test
6+
using SciMLBase
7+
using LinearAlgebra
8+
using FFTW
9+
10+
@testset "Thomson Scattering — Dipole Pattern" begin
11+
# Non-relativistic Thomson scattering from a single electron in a weak plane wave.
12+
# For a₀ ≪ 1 (linear polarization along x), the scattered radiation at the
13+
# fundamental frequency follows the dipole pattern: |A_ω|² ∝ sin²(θ) / R²
14+
# where θ is the angle from the polarization (x) axis.
15+
# Higher harmonics should be negligible.
16+
17+
# Use natural units (c = 1)
18+
@named ref_frame = ProperFrame(:natural)
19+
20+
a₀ = 0.01
21+
@named laser = PlaneWave(; amplitude = a₀, frequency = 1.0, ref_frame)
22+
@named elec = ClassicalElectron(; laser)
23+
sys = mtkcompile(elec)
24+
25+
ω = 1.0
26+
T = 2π / ω
27+
c = 1.0
28+
29+
# Simulate for several periods to get a clean FFT
30+
n_periods = 20
31+
τi = -n_periods * T / 2
32+
τf = n_periods * T / 2
33+
tspan = (τi, τf)
34+
35+
x⁰ = [τi * c, 0.0, 0.0, 0.0]
36+
u⁰ = [c, 0.0, 0.0, 0.0] # electron at rest
37+
38+
u0 = [sys.x => x⁰, sys.u => u⁰]
39+
prob = ODEProblem{false, SciMLBase.FullSpecialize}(
40+
sys, u0, tspan, u0_constructor = SVector{8}, fully_determined = true
41+
)
42+
sol = solve(prob, Vern9(), reltol = 1.0e-14, abstol = 1.0e-14)
43+
44+
@test SciMLBase.successful_retcode(sol)
45+
46+
# Build trajectory interpolant
47+
trajs = [TrajectoryInterpolant(sol, sys.x, sys.u)]
48+
49+
# Screen setup: observe from far field at distance Z along z
50+
Z = 1.0e4 # far field: Z ≫ λ = 2π
51+
δt = T / 8 # 8 samples per period (Nyquist up to 4th harmonic)
52+
53+
Nx, Ny = 15, 15
54+
L = Z * 0.1 # small screen so all pixels have similar arrival times
55+
56+
# Cover the full range: earliest signal (center pixel) to latest end (corner pixel)
57+
x⁰_earliest = c * τi + Z
58+
x⁰_latest = c * τf + hypot(Z, L * sqrt(2))
59+
N_samples = ceil(Int, (x⁰_latest - x⁰_earliest) / (c * δt))
60+
x⁰_samples = range(start = x⁰_earliest, step = c * δt, length = N_samples)
61+
62+
screen = ObserverScreen(
63+
LinRange(-L, L, Nx),
64+
LinRange(-L, L, Ny),
65+
Z,
66+
x⁰_samples
67+
)
68+
69+
A_s = accumulate_potential(trajs, screen, Tsit5())
70+
A_ω = rfft(A_s, 1)
71+
72+
# Find fundamental frequency bin
73+
freqs = rfftfreq(N_samples, 1 / δt)
74+
idx_f1 = findmin(f -> abs(f - ω / 2π), freqs)[2]
75+
76+
# Extract |A_ω|² at fundamental frequency (all 3 spatial components)
77+
intensity = zeros(Nx, Ny)
78+
A_ω_expected = zeros(Nx, Ny)
79+
for ix in 1:Nx, iy in 1:Ny
80+
intensity[ix, iy] = sum(abs2, A_ω[idx_f1, 2:4, ix, iy])
81+
r_obs = [screen.x_grid[ix], screen.y_grid[iy], Z]
82+
R = norm(r_obs)
83+
= r_obs / norm(r_obs)
84+
A_ω_expected[ix, iy] = (1 - n̂[1]^2) / R^2
85+
end
86+
87+
@testset "Dipole angular pattern" begin
88+
for ix in 1:Nx, iy in 1:Ny
89+
expected = A_ω_expected[ix, iy] / maximum(A_ω_expected)
90+
@test intensity[ix, iy] / maximum(intensity) expected rtol = 0.05
91+
end
92+
end
93+
94+
@testset "Higher harmonics suppressed" begin
95+
# For a₀ ≪ 1, the 2nd harmonic should be much weaker than fundamental
96+
idx_f2 = findmin(f -> abs(f - 2ω / 2π), freqs)[2]
97+
power_f1 = sum(abs2, A_ω[idx_f1, 2:4, :, :])
98+
power_f2 = sum(abs2, A_ω[idx_f2, 2:4, :, :])
99+
100+
@test power_f2 / power_f1 < 1.0e-4
101+
end
102+
end

0 commit comments

Comments
 (0)