diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0dee6a2..2590e60 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -15,7 +15,7 @@ jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} - timeout-minutes: 60 + timeout-minutes: 120 permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created actions: write contents: read @@ -37,6 +37,9 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 + - name: Monitor memory + if: always() + run: free -h - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v5 diff --git a/Project.toml b/Project.toml index 9a2867d..dc6ea89 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "0.1.0-DEV" authors = ["Sebastian Micluța-Câmpeanu and contributors"] [deps] +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" @@ -19,6 +20,8 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [compat] Aqua = "0.8" CSV = "0.10" +DataInterpolations = "8.9.0" +FFTW = "1" HypergeometricFunctions = "0.3.28" JET = "0.9, 0.10, 0.11" LaserTypes = "0.2" @@ -27,6 +30,7 @@ ModelingToolkit = "11" ModelingToolkitBase = "1.20.0" OrdinaryDiffEqNonlinearSolve = "1.10.0" OrdinaryDiffEqRosenbrock = "1.11.0" +OrdinaryDiffEqTsit5 = "1.9" OrdinaryDiffEqVerner = "1.2.0" PhysicalConstants = "0.2.3" Plots = "1" @@ -44,11 +48,13 @@ julia = "1.11" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LaserTypes = "e07c0bfa-524c-4f35-a151-c3dd916fa2f0" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -59,4 +65,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "CSV", "JET", "LaserTypes", "LinearAlgebra", "Test", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqVerner", "Plots", "Random", "SciMLBase", "StaticArrays", "Statistics", "SymbolicIndexingInterface"] +test = ["Aqua", "CSV", "FFTW", "JET", "LaserTypes", "LinearAlgebra", "Test", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Plots", "Random", "SciMLBase", "StaticArrays", "Statistics", "SymbolicIndexingInterface"] diff --git a/scripts/Project.toml b/scripts/Project.toml index 1b8513b..579ace0 100644 --- a/scripts/Project.toml +++ b/scripts/Project.toml @@ -1,10 +1,18 @@ [deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DiffEqGPU = "071ae1c0-96b5-11e9-1965-c90190d839ea" ElectronDynamicsModels = "acecdaf2-97b2-47e1-90eb-3efa7bb274e5" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2" +ProfileCanvas = "efd6af41-a80b-495e-886c-e51b0c7d77a3" +ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" diff --git a/scripts/thomson_scattering.jl b/scripts/thomson_scattering.jl new file mode 100644 index 0000000..0a231db --- /dev/null +++ b/scripts/thomson_scattering.jl @@ -0,0 +1,149 @@ +using ElectronDynamicsModels +using ModelingToolkit +using OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5 +using OrdinaryDiffEqNonlinearSolve +using SciMLBase +using StaticArrays +using SymbolicIndexingInterface +using LinearAlgebra +using FFTW +using CairoMakie + +# Atomic units +const c = 137.03599908330932 + +# Laser parameters +ω = 0.057 +τ = 150 / ω +λ = 2π * c / ω +w₀ = 75λ +Rmax = 3.25w₀ + +a₀ = 10.0 + +@named ref_frame = ProperFrame(:atomic) + +@named laser = LaguerreGaussLaser(; + wavelength = λ, + a0 = a₀, + beam_waist = w₀, + radial_index = 2, + azimuthal_index = -2, + ref_frame, + temporal_profile = :gaussian, + temporal_width = τ, + focus_position = 0.0, + polarization = :circular +) +@named elec = ClassicalElectron(; laser) +sys = mtkcompile(elec) + +# Time span +τi = -8τ +τf = 8τ +tspan = (τi, τf) + +# Single electron solve (for parameter access) +x⁰ = [τi * c, 0.0, 0.0, 0.0] +u⁰ = [c, 0.0, 0.0, 0.0] + +u0 = [ + sys.x => x⁰, + sys.u => u⁰, +] + +prob = ODEProblem{false, SciMLBase.FullSpecialize}( + sys, u0, tspan, u0_constructor = SVector{8}, fully_determined = true +) +sol0 = solve(prob, Vern9(), reltol = 1.0e-15, abstol = 1.0e-12) + +# Sunflower distribution for electron positions +const ϕ = (1 + √5) / 2 + +function radius(k, n, b) + if k > n - b + return 1.0 + else + return sqrt(k - 0.5) / sqrt(n - (b + 1) / 2) + end +end + +function sunflower(n, α) + points = Vector{Vector{Float64}}() + angle_stride = 2π / ϕ^2 + b = round(Int, α * sqrt(n)) + for k in 1:n + r = radius(k, n, b) + θ = k * angle_stride + push!(points, [r * cos(θ), r * sin(θ)]) + end + return points +end + +# Ensemble solve +N = 300 +R₀ = Rmax * sunflower(N, 2) +xμ = [[τi * c, r..., 0.0] for r in R₀] + +set_x = setsym_oop(prob, [Initial(sys.x); Initial(sys.u)]) + +function prob_func(prob, i, repeat) + x_new = SVector{4}(xμ[i]...) + u_new = SVector{4}(c, 0.0, 0.0, 0.0) + u0, p = set_x(prob, SVector{8}(x_new..., u_new...)) + return remake(prob; u0, p) +end + +function abserr(a₀) + amp = log10(a₀) + expo = -amp^2 / 27 + 32amp / 27 - 220 / 27 + return 10^expo +end + +ensemble = EnsembleProblem(prob; prob_func, safetycopy = false) +solution = solve( + ensemble, Vern9(), EnsembleThreads(); + reltol = 1.0e-12, abstol = abserr(a₀), trajectories = N +) + +# Radiation computation + +trajs = trajectory_interpolants(solution) + +# Screen parameters +const Z = 2.0e5λ +const δt = 2π / ω / 4 +const N_samples = floor(Int, (τf - τi) / δt) +const x⁰_start = c * τi + hypot(Z, 25w₀ + Rmax) + +Nx = 50 +Ny = 50 + +x⁰_samples = range(start = x⁰_start, step = c * δt, length = N_samples) + +screen = ObserverScreen( + LinRange(-25w₀, 25w₀, Nx), + LinRange(-25w₀, 25w₀, Ny), + Z, + x⁰_samples +) + +A_s = accumulate_potential(trajs, screen, Tsit5()) + +# GPU +# accumulate_potential(trajs, screen, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend())) +A_ω = rfft(A_s, 1) + +# Find fundamental frequency bin +freqs = rfftfreq(N_samples, 1 / δt) +idx_f1 = findmin(x -> abs(x - ω / 2π), freqs)[2] + +# Plot the y-component of the potential at the fundamental frequency +fig = Figure() +ax = Axis(fig[1, 1], aspect = 1, xlabel = "x", ylabel = "y", title = "Thomson scattering (ω₁)") +field = real.(A_ω[idx_f1, 3, :, :]) +heatmap!( + ax, collect(screen.x_grid), collect(screen.y_grid), field, + colorrange = maximum(abs, field) .* (-1, 1), colormap = :seismic +) +fig diff --git a/src/ElectronDynamicsModels.jl b/src/ElectronDynamicsModels.jl index 37b58ff..91dd43b 100644 --- a/src/ElectronDynamicsModels.jl +++ b/src/ElectronDynamicsModels.jl @@ -2,32 +2,32 @@ module ElectronDynamicsModels using ModelingToolkit using ModelingToolkitBase: AbstractSystem, SymbolicT, build_explicit_observed_function, get_systems -using SymbolicIndexingInterface: getname, setsym_oop +using SymbolicIndexingInterface: getname, setsym_oop, variable_index using PhysicalConstants, Unitful, UnitfulAtomic using PhysicalConstants.CODATA2018: c_0, e, m_e, ε_0 using LinearAlgebra using Symbolics using HypergeometricFunctions: HypergeometricFunctions, _₁F₁, pochhammer using StaticArrays +using SciMLBase +using DataInterpolations m_dot(x, y) = x[1] * y[1] - x[2] * y[2] - x[3] * y[3] - x[4] * y[4] export GaussLaser, LaguerreGaussLaser export ReferenceFrame, ProperFrame, LabFrame, - UniformField, - PlaneWave, + UniformField, PlaneWave, ParticleDynamics, - LandauLifshitzRadiation, - AbrahamLorentzRadiation, + LandauLifshitzRadiation, AbrahamLorentzRadiation, ChargedParticle, - ClassicalElectron, - RadiatingElectron, - LandauLifshitzElectron, - FieldEvaluator + ClassicalElectron, RadiatingElectron, LandauLifshitzElectron, + FieldEvaluator, + ObserverScreen, trajectory_interpolants, TrajectoryInterpolant, accumulate_potential include("base.jl") include("dynamics.jl") include("fields.jl") +include("radiation.jl") include("radiation_reaction.jl") include("external_fields.jl") include("systems.jl") diff --git a/src/external_fields.jl b/src/external_fields.jl index dd9e0da..4ed2750 100644 --- a/src/external_fields.jl +++ b/src/external_fields.jl @@ -13,9 +13,11 @@ Parameters: - t₀: pulse center time (= n_cycles × 2π/ω) - z₀: focus position along z-axis """ -@component function GaussLaser(; name, wavelength = nothing, frequency = nothing, +@component function GaussLaser(; + name, wavelength = nothing, frequency = nothing, a0 = 10.0, beam_waist = nothing, polarization = :linear, - n_cycles = 5, focus_position = nothing, ref_frame) + n_cycles = 5, focus_position = nothing, ref_frame + ) if wavelength === nothing && frequency === nothing wavelength = 1.0 end @@ -82,8 +84,7 @@ Parameters: E[1] ~ real(ξx * E_g) E[2] ~ real(ξy * E_g) E[3] ~ real( - 2im / (k * wz^2) * - (1 + im * (z / z_R)) * + 2im / (k * wz^2) * (1 + im * (z / z_R)) * (x[2] * (ξx * E_g) + x[3] * (ξy * E_g)), ) @@ -228,7 +229,7 @@ Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970) vars = nameof(iv) == :τ ? [x, t] : [x] - sys = System(eqs, iv, vars, [A, ω, k_dir, pol, λ]; name, systems = [ref_frame], initialization_eqs, bindings, initial_conditions) + sys = System(eqs, iv; name, systems = [ref_frame], initialization_eqs, bindings, initial_conditions) extend(sys, field_dynamics) end diff --git a/src/radiation.jl b/src/radiation.jl new file mode 100644 index 0000000..d403a15 --- /dev/null +++ b/src/radiation.jl @@ -0,0 +1,207 @@ +# Trajectory access (thread-safe interpolation wrapper) +struct TrajectoryInterpolant{I, R, U, T} + itp::I # DataInterpolations interpolant → SVector{8} + x_idxs::R # SVector{4, Int} indices for xμ = (x⁰, x¹, x², x³) + u_idxs::U # SVector{4, Int} indices for uμ = (u⁰, u¹, u², u³) + K::T # q_e / (4π ε₀ c) — Liénard-Wiechert prefactor +end + +function TrajectoryInterpolant(sol::SciMLBase.AbstractODESolution, x_syms, u_syms) + x_idxs = SVector{4, Int}(variable_index.((sol,), collect(x_syms))) + u_idxs = SVector{4, Int}(variable_index.((sol,), collect(u_syms))) + itp = CubicSpline(sol.u, sol.t; extrapolation = ExtrapolationType.Extension) + sys = sol.prob.f.sys + _ref_frame = _find_ref_frame(sys) + K = sol.ps[_ref_frame.q_e / (4π * _ref_frame.ε₀ * _ref_frame.c)] + return TrajectoryInterpolant(itp, x_idxs, u_idxs, K) +end + +function (t::TrajectoryInterpolant)(τ) + v = t.itp(τ) + xμ = v[t.x_idxs] + uμ = v[t.u_idxs] + return (xμ, uμ) +end + +""" + trajectory_interpolants(ensemble_sol::EnsembleSolution) + +Return a vector of `TrajectoryInterpolant` objects, which interpolate the trajectories +in a fast and thread safe manner. +""" +function trajectory_interpolants(ensemble_sol::EnsembleSolution) + sys = ensemble_sol.u[1].prob.f.sys + return [TrajectoryInterpolant(ensemble_sol.u[i], sys.x, sys.u) for i in axes(ensemble_sol.u, 1)] +end + +# Observer geometry + temporal sampling +struct ObserverScreen{G, T, R} + x_grid::G # e.g., LinRange for x + y_grid::G # e.g., LinRange for y + z::T # screen distance + x⁰_samples::R # uniform observer-time sampling grid +end + +function Base.show(io::IO, s::ObserverScreen) + Nx, Ny = length(s.x_grid), length(s.y_grid) + N = length(s.x⁰_samples) + δx⁰ = N > 1 ? step(s.x⁰_samples) : 0.0 + return print(io, "ObserverScreen($(Nx)×$(Ny) pixels, z=$(s.z), $(N) time samples, Δx⁰=$(δx⁰))") +end + +function Base.show(io::IO, ::MIME"text/plain", s::ObserverScreen) + Nx, Ny = length(s.x_grid), length(s.y_grid) + N = length(s.x⁰_samples) + δx⁰ = N > 1 ? step(s.x⁰_samples) : 0.0 + println(io, "ObserverScreen") + println(io, " pixels: $(Nx) × $(Ny)") + println(io, " x range: [$(first(s.x_grid)), $(last(s.x_grid))]") + println(io, " y range: [$(first(s.y_grid)), $(last(s.y_grid))]") + println(io, " z: $(s.z)") + println(io, " time samples: $(N)") + println(io, " x⁰ range: [$(first(s.x⁰_samples)), $(last(s.x⁰_samples))]") + return print(io, " Δx⁰: $(δx⁰)") +end + +# dτᵣ/dt = 1/(u⁰(τᵣ) - u⃗(τᵣ)·n̂(τᵣ, r_obs)) +function retarded_time_rhs(τᵣ, p, t) + traj, r_obs = p + xμ, uμ = traj(τᵣ) + xⁱ = xμ[SA[2, 3, 4]] + n̂ = (r_obs - xⁱ) * inv(norm(r_obs - xⁱ)) + return inv(uμ[1] - uμ[SA[2, 3, 4]] ⋅ n̂) +end + +function advanced_time(traj, τ, r_obs) + xμ, _ = traj(τ) + return xμ[1] + norm(r_obs - xμ[SA[2, 3, 4]]) +end + +function retarded_time_problem(traj::TrajectoryInterpolant, screen::ObserverScreen) + Nx, Ny = length(screen.x_grid), length(screen.y_grid) + CI = CartesianIndices((Nx, Ny)) + + r_obs_0 = SVector{3}(screen.x_grid[1], screen.y_grid[1], screen.z) + τi = first(traj.itp.t) + τf = last(traj.itp.t) + (x⁰_i_0, x⁰_f_0) = advanced_time(traj, τi, r_obs_0), advanced_time(traj, τf, r_obs_0) + + prob = ODEProblem{false, SciMLBase.FullSpecialize}( + retarded_time_rhs, + τi, + (x⁰_i_0, x⁰_f_0), + (traj, r_obs_0), + ) + + function set_pixel(prob, i, repeat) + ix, iy = Tuple(CI[i]) + r_obs = SVector{3}(screen.x_grid[ix], screen.y_grid[iy], screen.z) + (x⁰_i, x⁰_f) = advanced_time(traj, τi, r_obs), advanced_time(traj, τf, r_obs) + return remake(prob; p = (traj, r_obs), tspan = (x⁰_i, x⁰_f)) + end + + return EnsembleProblem(prob; prob_func = set_pixel, safetycopy = false) +end + +""" + accumulate_potential(trajs, screen, alg; solve_kwargs...) + accumulate_potential(trajs, screen, alg, ensemblealg; solve_kwargs...) + +Compute the Liénard-Wiechert 4-potential on `screen` from electron `trajs`. + +For each electron trajectory, solves the retarded-time ODE to map observer time to +proper time, then evaluates `Aμ = K uμ / (xr · u)` at uniform observer-time samples. +Returns `A[k, μ, ix, iy]` — the time-domain 4-potential ready for FFT. + +The two-argument `alg` form uses a `reinit!`-based integrator pool for efficient CPU +threading. The four-argument form with `ensemblealg` uses `EnsembleProblem` for +compatibility with GPU backends (e.g., `EnsembleGPUKernel`). + +# Arguments +- `trajs`: vector of `TrajectoryInterpolant` from [`trajectory_interpolants`](@ref) +- `screen`: `ObserverScreen` defining pixel grid and observer-time samples +- `alg`: ODE solver algorithm for the retarded-time problem (e.g., `Tsit5()`) +- `ensemblealg`: (optional) ensemble algorithm (e.g., `EnsembleGPUKernel(backend)`) +- `solve_kwargs...`: additional keyword arguments passed to the ODE solver +""" +function accumulate_potential(trajs::Vector{<:TrajectoryInterpolant}, screen::ObserverScreen, alg; solve_kwargs...) + x⁰_samples = screen.x⁰_samples + N_samples = length(x⁰_samples) + Nx, Ny = length(screen.x_grid), length(screen.y_grid) + + A = zeros(N_samples, 4, Nx, Ny) + + for (j, traj) in enumerate(trajs) + τi = first(traj.itp.t) + τf = last(traj.itp.t) + + r_obs_0 = SVector{3}(screen.x_grid[1], screen.y_grid[1], screen.z) + x⁰_i_0 = advanced_time(traj, τi, r_obs_0) + x⁰_f_0 = advanced_time(traj, τf, r_obs_0) + proto_prob = ODEProblem{false, SciMLBase.FullSpecialize}( + retarded_time_rhs, τi, (x⁰_i_0, x⁰_f_0), (traj, r_obs_0) + ) + + nworkers = Threads.nthreads() + integ_pool = Channel{Any}(nworkers) + for _ in 1:nworkers + put!(integ_pool, init(proto_prob, alg; saveat = x⁰_samples, solve_kwargs...)) + end + + Threads.@threads for ix in Base.OneTo(Nx) + integ = take!(integ_pool) + for iy in Base.OneTo(Ny) + r_obs = SVector{3}(screen.x_grid[ix], screen.y_grid[iy], screen.z) + x⁰_i = advanced_time(traj, τi, r_obs) + x⁰_f = advanced_time(traj, τf, r_obs) + + integ.p = (traj, r_obs) + reinit!(integ, τi; t0 = x⁰_i, tf = x⁰_f) + solve!(integ) + + _accumulate_pixel!(A, traj, screen, integ.sol.u, ix, iy) + end + put!(integ_pool, integ) + end + end + + return A +end + +# Ensemble path (for GPU / custom ensemble algorithms) +function accumulate_potential(trajs::Vector{<:TrajectoryInterpolant}, screen::ObserverScreen, alg, ensemblealg; solve_kwargs...) + x⁰_samples = screen.x⁰_samples + N_samples = length(x⁰_samples) + Nx, Ny = length(screen.x_grid), length(screen.y_grid) + + A = zeros(N_samples, 4, Nx, Ny) + + for (j, traj) in enumerate(trajs) + rt_prob = retarded_time_problem(traj, screen) + N_pixels = Nx * Ny + rt_sol = solve( + rt_prob, alg, ensemblealg; + trajectories = N_pixels, saveat = x⁰_samples, solve_kwargs... + ) + LI = LinearIndices((Nx, Ny)) + + Threads.@threads for ix in Base.OneTo(Nx) + for iy in Base.OneTo(Ny) + _accumulate_pixel!(A, traj, screen, rt_sol.u[LI[ix, iy]].u, ix, iy) + end + end + end + + return A +end + +function _accumulate_pixel!(A, traj, screen, τ_samples, ix, iy) + r_obs = SVector{3}(screen.x_grid[ix], screen.y_grid[iy], screen.z) + for (k, τ) in enumerate(τ_samples) + xμ, uμ = traj(τ) + disp = r_obs - xμ[SA[2, 3, 4]] + xr = SVector{4}(norm(disp), disp[1], disp[2], disp[3]) + @views A[k, :, ix, iy] .+= traj.K * uμ ./ m_dot(xr, uμ) + end + return +end diff --git a/src/systems.jl b/src/systems.jl index 0933471..5702cc2 100644 --- a/src/systems.jl +++ b/src/systems.jl @@ -1,18 +1,23 @@ """ - _find_ref_frame(external_field) + _find_ref_frame(sys) Find the reference frame subsystem within an external field by checking for the metric tensor `gμν` and electron mass `m_e` parameters. """ -function _find_ref_frame(external_field::AbstractSystem) - for s in get_systems(external_field) +function _find_ref_frame(sys::AbstractSystem) + if isempty(get_systems(sys)) && ModelingToolkitBase.has_parent(sys) + return _find_ref_frame(ModelingToolkitBase.get_parent(sys)) + end + for s in get_systems(sys) s isa AbstractSystem || continue names = getname.(parameters(s)) if :gμν ∈ names && :m_e ∈ names return s + elseif !isempty(get_systems(s)) + return _find_ref_frame(s) end end - error("No reference frame found in $(nameof(external_field))") + error("No reference frame found in $(nameof(sys))") end @component function ChargedParticle(; diff --git a/test/test_classical_electron.jl b/test/classical_electron.jl similarity index 100% rename from test/test_classical_electron.jl rename to test/classical_electron.jl diff --git a/test/test_gauss_laser.jl b/test/gauss_laser.jl similarity index 100% rename from test/test_gauss_laser.jl rename to test/gauss_laser.jl diff --git a/test/test_laguerre_gauss.jl b/test/laguerre_gauss.jl similarity index 100% rename from test/test_laguerre_gauss.jl rename to test/laguerre_gauss.jl diff --git a/test/test_radiation.jl b/test/radiation.jl similarity index 100% rename from test/test_radiation.jl rename to test/radiation.jl diff --git a/test/runtests.jl b/test/runtests.jl index b3d591e..e253ba3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,18 +8,25 @@ using JET Aqua.test_all(ElectronDynamicsModels) end @testset "Code linting (JET.jl)" begin - JET.test_package(ElectronDynamicsModels; target_modules = (ElectronDynamicsModels,)) + if VERSION < v"1.12" + JET.test_package(ElectronDynamicsModels; target_modules = (ElectronDynamicsModels,)) + else + @test_broken false # JET causes CI termination on 1.12 + end end @testset "Classical Electron" begin - include("test_classical_electron.jl") + include("classical_electron.jl") + end + @testset "Radiation Emission" begin + include("thomson_scattering.jl") end @testset "Radiation Reaction" begin - include("test_radiation.jl") + include("radiation.jl") end @testset "GaussLaser" begin - include("test_gauss_laser.jl") + include("gauss_laser.jl") end @testset "LaguerreGauss" begin - include("test_laguerre_gauss.jl") + include("laguerre_gauss.jl") end end diff --git a/test/thomson_scattering.jl b/test/thomson_scattering.jl new file mode 100644 index 0000000..1392c36 --- /dev/null +++ b/test/thomson_scattering.jl @@ -0,0 +1,102 @@ +using ElectronDynamicsModels +using ModelingToolkit +using OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5 +using StaticArrays +using Test +using SciMLBase +using LinearAlgebra +using FFTW + +@testset "Thomson Scattering — Dipole Pattern" begin + # Non-relativistic Thomson scattering from a single electron in a weak plane wave. + # For a₀ ≪ 1 (linear polarization along x), the scattered radiation at the + # fundamental frequency follows the dipole pattern: |A_ω|² ∝ sin²(θ) / R² + # where θ is the angle from the polarization (x) axis. + # Higher harmonics should be negligible. + + # Use natural units (c = 1) + @named ref_frame = ProperFrame(:natural) + + a₀ = 0.01 + @named laser = PlaneWave(; amplitude = a₀, frequency = 1.0, ref_frame) + @named elec = ClassicalElectron(; laser) + sys = mtkcompile(elec) + + ω = 1.0 + T = 2π / ω + c = 1.0 + + # Simulate for several periods to get a clean FFT + n_periods = 20 + τi = -n_periods * T / 2 + τf = n_periods * T / 2 + tspan = (τi, τf) + + x⁰ = [τi * c, 0.0, 0.0, 0.0] + u⁰ = [c, 0.0, 0.0, 0.0] # electron at rest + + u0 = [sys.x => x⁰, sys.u => u⁰] + prob = ODEProblem{false, SciMLBase.FullSpecialize}( + sys, u0, tspan, u0_constructor = SVector{8}, fully_determined = true + ) + sol = solve(prob, Vern9(), reltol = 1.0e-14, abstol = 1.0e-14) + + @test SciMLBase.successful_retcode(sol) + + # Build trajectory interpolant + trajs = [TrajectoryInterpolant(sol, sys.x, sys.u)] + + # Screen setup: observe from far field at distance Z along z + Z = 1.0e4 # far field: Z ≫ λ = 2π + δt = T / 8 # 8 samples per period (Nyquist up to 4th harmonic) + + Nx, Ny = 15, 15 + L = Z * 0.1 # small screen so all pixels have similar arrival times + + # Cover the full range: earliest signal (center pixel) to latest end (corner pixel) + x⁰_earliest = c * τi + Z + x⁰_latest = c * τf + hypot(Z, L * sqrt(2)) + N_samples = ceil(Int, (x⁰_latest - x⁰_earliest) / (c * δt)) + x⁰_samples = range(start = x⁰_earliest, step = c * δt, length = N_samples) + + screen = ObserverScreen( + LinRange(-L, L, Nx), + LinRange(-L, L, Ny), + Z, + x⁰_samples + ) + + A_s = accumulate_potential(trajs, screen, Tsit5()) + A_ω = rfft(A_s, 1) + + # Find fundamental frequency bin + freqs = rfftfreq(N_samples, 1 / δt) + idx_f1 = findmin(f -> abs(f - ω / 2π), freqs)[2] + + # Extract |A_ω|² at fundamental frequency (all 3 spatial components) + intensity = zeros(Nx, Ny) + A_ω_expected = zeros(Nx, Ny) + for ix in 1:Nx, iy in 1:Ny + intensity[ix, iy] = sum(abs2, A_ω[idx_f1, 2:4, ix, iy]) + r_obs = [screen.x_grid[ix], screen.y_grid[iy], Z] + R = norm(r_obs) + n̂ = r_obs / norm(r_obs) + A_ω_expected[ix, iy] = (1 - n̂[1]^2) / R^2 + end + + @testset "Dipole angular pattern" begin + for ix in 1:Nx, iy in 1:Ny + expected = A_ω_expected[ix, iy] / maximum(A_ω_expected) + @test intensity[ix, iy] / maximum(intensity) ≈ expected rtol = 0.05 + end + end + + @testset "Higher harmonics suppressed" begin + # For a₀ ≪ 1, the 2nd harmonic should be much weaker than fundamental + idx_f2 = findmin(f -> abs(f - 2ω / 2π), freqs)[2] + power_f1 = sum(abs2, A_ω[idx_f1, 2:4, :, :]) + power_f2 = sum(abs2, A_ω[idx_f2, 2:4, :, :]) + + @test power_f2 / power_f1 < 1.0e-4 + end +end