ElectronDynamicsModels

FieldAccumulator(screen, backend; mode = Val(:split))

The device-resident accumulation buffers of accumulate_field, held as an object so that they survive one call and several calls can sum into the same buffers.

By default (buffers = nothing, finish = true) accumulate_field allocates a set per call and downloads it at the end — the single-call API, unchanged. With finish = false it instead returns the FieldAccumulator holding that call's electrons — after draining the call's launches, so the buffers are safe to hand to another task — and passing it back as buffers on the next call continues the accumulation in place; the cube is downloaded once, by the call that runs with finish = true (or by finish_field).

This is the memory model behind the solver scripts' EDM_ELECTRON_BATCH knob: the host then holds one batch of trajectory splines instead of every electron's for the whole field phase, while the device keeps the same single buffer set as an unbatched run. The kernels, the per-electron uploads and the launch order are untouched, so on one device a batched run adds the same per-electron contributions in the same order as the single call (bit-identical); across devices the shard composition changes with the batching, which moves the sum by the last bits.

Fields: E1, B1 (far field in :split, total in :total), E2, B2 (near field; aliases of E1, B1 in :total, where the kernel sums far + near before the write and never touches them), the mode, the device the buffers live on, and the running n_electrons count.

acc = accumulate_field(trajs[1:500], screen, alg, backend; finish = false)
acc = accumulate_field(trajs[501:1000], screen, alg, backend; buffers = acc, finish = false)
fld = accumulate_field(trajs[1001:1500], screen, alg, backend; buffers = acc)   # downloads once

FieldEvaluator(laser)

Evaluate electromagnetic fields E(x,y,z,t) and B(x,y,z,t) at arbitrary spacetime points without running an ODE solver. Internally wraps the laser in a ClassicalElectron to produce a compilable system, then uses build_explicit_observed_function to get a fast compiled function for the field observables.

Usage

The input is [t, x, y, z] (bare time, not c*t). The speed of light c is obtained from the reference frame and applied internally.

@named world = Worldline(:τ, :atomic)
@named laser = GaussLaser(; wavelength=1.0, a0=1.0, world)
fe = FieldEvaluator(laser)
result = fe([t, x, y, z])  # (E = [...], B = [...])

GPUCubicSpline{V, M}

A GPU-compatible cubic spline interpolant. Stores precomputed coefficients in flat arrays (works with both Vector and CuArray). The coefficient matrices are column-major over the knots (N × D: for one component the knots are contiguous), because the kernels' reuse of a cache line is along the knot index — consecutive samples of a pixel and the lanes of a wave sit on the same or the neighbouring knot. The knot-major alternative (D × N) issues fewer, wider loads but turns a cache line over every two knots and measured 20 % slower on the MI300X; see the manual page "The device spline".

Evaluates the standard natural cubic spline formula: S(t) = (z[i]dt2³ + z[i+1]dt1³) / (6h[i+1]) + c1[i]dt1 + c2[i]*dt2 where dt1 = t - t[i], dt2 = t[i+1] - t.

Fields

• t: knot positions (length N)
• h: interval widths, h[i] = t[i] - t[i-1] (length N, h[1] = 0)
• z: second derivatives at knots (N × D matrix, D = number of components)
• c1: precomputed linear coefficients per interval ((N-1) × D matrix)
• c2: precomputed linear coefficients per interval ((N-1) × D matrix)

(spline::GPUCubicSpline)(τ)
(spline::GPUCubicSpline)(τ, guess) -> (value, idx)

Evaluate the spline at time τ, returning an SVector of interpolated values. The two-argument form starts the knot search at interval guess (see _searchsorted_left) and also returns the interval it used, to be fed back as the next call's guess; its value is bit-identical to the one-argument form.

GPUCubicSpline(itp::DataInterpolations.CubicSpline)

Construct a GPUCubicSpline by extracting precomputed coefficients from an existing DataInterpolations.CubicSpline.

GPUKernelNewton

Sentinel solver type selecting the Newton light-cone GPU kernel path: the retarded proper time at each saveat slot is obtained by solving the light-cone condition x⁰_target − x⁰(τ) − |r_obs − x⃗(τ)| = 0 with n_iters warm-started Newton corrections, instead of marching the retarded-time ODE as GPUKernelRK4 does. The residual is evaluated in the equivalent screen-relative light-front spelling f = tₖ − ψ(τ) − ρ²/(R + d³) (see the file header), which keeps Float64 rounding at the interaction scale instead of the ε·Z floor of the absolute coordinates.

Per-slot errors are independent across observer-time slots (no accumulation). Prefer this kernel for boosted forward-scattering geometries (γ ≫ 1), where n_iters = 1–2 beats matched-accuracy RK4 by ~4× in cost; prefer GPUKernelRK4 for rest-electron field maps and deep-floor harmonic B-map studies, where its smoother error spectrum preserves near-floor ring structure. Measured campaign numbers in the header of gpu/kernel_newton.jl.

GPUKernelRK4

Sentinel solver type that selects the unified GPU kernel path.

IntervalCoefs{D, T}

Everything a cubic-spline evaluation needs on one knot interval: its end knots, inv(6h) and the four D-vectors z[i], z[i+1], c1[i], c2[i] (3 + 4D scalars). _fetch_interval loads it once; _eval_poly evaluates the cubic at any τ from it. The per-slot solvers can hold it in registers across the Newton corrections (or RK4 stages) of a slot, which leave the interval only at a knot crossing, instead of reloading the coefficients for every evaluation (coef_reuse = Val(true) in the accumulate functions: bit-identical values, fewer loads, more live registers).

ObserverScreen(x_grid, y_grid, z, x⁰_samples; c)

Keyword form requiring the speed of light c in the working units (e.g. getdefault(world.c)). c has no default on purpose: the x⁰ = c·t axis is meaningless without the unit system, and a wrong default would silently corrupt δt, FFT frequencies, and recommended_n_substeps.

ReduceStats()

Wall-clock accounting of the multi-device reduce, filled in by accumulate_field_sharded when passed as reduce_stats: fold_s is the time the finishing tasks spent folding their partials into the accumulator (under the lock; device-to-device adds for reduce = :device, host download-permute-add for :host), n_folds how many partials were folded, and download_s the final download and permute into the host cube (:device only; the :host path downloads as it folds). The solver scripts record them as [timing] reduce_fold / reduce_download, next to field and kernel, so a sharded cell's non-kernel time is visible.

ShardedFieldAccumulator(devices)

Per-device FieldAccumulators of a batched accumulate_field_sharded run: one buffer set per device, each holding that device's shard of every batch accumulated so far. Created by the first finish = false call and passed back as buffers; the device reduce runs once, on the call with finish = true. accs[i] is nothing until device i has drawn its first shard.

Gaussian laser pulse electromagnetic field.

Represents a focused Gaussian beam with a temporal envelope. The beam propagates along the z-direction with waist w₀ at focus.

Parameters:

• λ: wavelength
• a₀: normalized vector potential (a0 kwarg)
• w₀: beam waist (defaults to 75λ)
• n_cycles: number of optical cycles to pulse center (determines t₀)
• τ0: temporal envelope half-width
• t₀: pulse center time (= n_cycles × 2π/ω)
• z₀: focus position along z-axis

Laguerre-Gauss laser beam electromagnetic field.

Represents a focused beam with orbital angular momentum (OAM). The beam is characterized by radial index p and azimuthal index m.

Parameters:

• λ: wavelength
• a₀: normalized vector potential (a0 kwarg)
• ϕ₀: initial carrier phase (initial_phase kwarg, defaults to 0.0)
• w₀: beam waist (defaults to 75λ)
• radial_index (p): radial mode number, p ≥ 0
• azimuthal_index (m): azimuthal mode number (orbital angular momentum)
• temporal_profile: :gaussian (pulsed) or :constant (CW)
• temporal_width: pulse width for gaussian profile (defaults to 100.0)
• focus_position: focal position along z-axis (defaults to 0.0)
• k_direction: propagation direction unit vector. Currently restricted to [0, 0, 1] (default, +ẑ) or [0, 0, -1] (−ẑ).

Reference: Allen et al., Phys. Rev. A 45, 8185 (1992)

Landau-Lifshitz Radiation Reaction

Implements the Landau-Lifshitz formulation of radiation reaction for a charged particle. This avoids the runaway solutions of Abraham-Lorentz by eliminating second derivatives.

References:

• Landau, L.D. & Lifshitz, E.M. "The Classical Theory of Fields" §76
• Niel et. al. 2018, 10.1103/PhysRevE.97.043209, eq. 2 and 3

The Landau-Lifshitz equation in covariant form: dp^μ/dτ = q F^μν uν + (2τₑ/3)[q ∂ν F^μλ u^ν uλ + q²/m F^μν Fνλ u^λ + q²/mc² (F^νλ u^λ) (F_νγ u^γ) u^μ]

Plane wave electromagnetic field.

For a plane wave with normalized vector potential a₀ = eA/(mc²), electrons can exhibit figure-8 motion when a₀ ~ 1.

Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970)

Uniform electromagnetic field component.

In crossed E and B fields with E⊥B and |E| < |B|c, particles drift with velocity v_drift = E×B/B²

Reference: Jackson, "Classical Electrodynamics", Section 12.4

Worldline(parameter, units; name)

Construct a worldline parameterization context for relativistic dynamics. The returned MTK system carries the Minkowski metric, physical constants, and the chosen integration parameter (independent variable).

parameter is :τ (proper time) or :t (lab time). Both choices describe the same inertial (lab) frame; they differ only in which scalar parameterizes the worldline that the integrator steps along. They are not different Lorentz frames — the components of x^μ, u^μ, F^{μν} are always resolved in the frame where the external fields are defined.

units is :SI, :atomic, or :natural.

_device_add!(backend, dev, dst, src) -> dst

Fold src, resident on the calling task's current device, into dst, resident on device dev (1-based vendor id), without a host copy of either. Same device (or the CPU backend): a chunked in-place broadcast. Different devices: each 1/16 chunk is copied device-to-device into a staging buffer on dev and added there; the vendor runtime routes the copy over the peer link when the two GPUs can address each other and through host memory otherwise. Should the direct copy be refused by the runtime, the chunk goes through a host slab instead. The calling task's device is restored on return.

_download_permuted(buf; backend = nothing, dev = 0, workers = 1) -> Array{T,4}

Download an accumulation buffer laid out [ix, iy, μ, k] (pixel-fastest for coalesced device writes; observer slot k slowest) into the cube layout [k, μ, ix, iy] without a full-size intermediate. The old permutedims(Array(buf), (4, 3, 1, 2)) held cube + copy on the host — a 2× peak that capped carrier-resolved cube designs at ~half the node's RAM (the permute lives on the host because the >2³¹-element GPU permutedims overflows its 32-bit linear indices, cudaError 700). Downloading contiguous k-chunks through the linear copyto! DMA path and permutedims!-ing each into a view of the preallocated cube keeps the host peak at ~(1 + 1/16)× cube; device memory and kernels are untouched, and every transfer stays far below 2³¹ elements. Works unchanged on the CPU backend (buf::Array).

workers > 1 spreads the 16 chunks over that many tasks, each with its own staging slab (host peak ~(1 + workers/16)× cube): the download DMA and the single-threaded permutedims! of the chunks then overlap. The tasks pin themselves to dev through gpu_device!(backend, dev), so backend and dev are required for that path.

_download_permuted_add!(out, buf) -> out

In-place sibling of _download_permuted: download buf ([ix, iy, μ, k]) chunk by chunk and ADD it into the preallocated host cube out ([k, μ, ix, iy]). Two staging buffers of 1/16 cube (the linear download slab and its permuted copy) are the only transient, so folding a device partial into a running host sum costs ~(1/8)× cube of extra host memory — the primitive behind the streamed multi-device reduce, which keeps ONE host cube resident instead of one per device.

_find_world(sys)

Find the Worldline subsystem within an external field by checking for the metric tensor gμν and electron mass m_e parameters.

_kz_sign(k_direction) -> ±1

Validate that k_direction is along ±ẑ and return its sign. Currently only ±ẑ propagation is supported for LaguerreGaussLaser; oblique directions require defining a transverse basis convention that has not yet been implemented.

_searchsorted_left(t, x, guess)

Warm-started variant: returns exactly what _searchsorted_left(t, x) returns (the largest i ∈ [1, length(t)-1] with t[i] ≤ x, or 1), but starts from guess. If x lies in [t[guess], t[guess+1]) the answer is guess after two knot reads; otherwise the search gallops from guess towards x (steps 1, 2, 4, …) and finishes with a binary search inside the bracket it found — about 2·log₂(distance) reads instead of log₂(N). The per-pixel kernels feed the interval of their previous evaluation back in: consecutive slots move the retarded time by a fraction of a knot, so the fast path is the common case. The result is identical by construction (same predicate, same bracket invariants), so kernel output does not change.

_searchsorted_left(t, x)

Binary search for the interval index: find largest i such that t[i] ≤ x. Clamps to [1, length(t)-1] for evaluation safety. GPU-compatible: no allocations, no dynamic dispatch.

a0_from_peak(a_peak; mode) -> a₀

EDM a0 for the LG mode mode = (p, m) whose focal-plane peak field equals that of a Gaussian beam with a0 = a_peak (same polarization): a_peak / lg_peak_factor(p, m).

Use it to match a measured peak intensity, e.g. a_peak = 0.855·λ[μm]·√(I / 10¹⁸ W cm⁻²); use a0_from_pulse_energy to match a pulse energy instead. At m = 7 the two targets and the bare a0 differ by orders of magnitude.

a0_from_pulse_energy(W, w₀, τ₀, ω; world, mode = (p = 0, m = 2)) -> a₀

Compute the dimensionless vector potential a₀ for a Laguerre-Gauss pulse of total energy W, host-beam waist w₀, field-envelope half-width τ₀, and angular frequency ω. All inputs in the unit system of world (atomic, SI, or natural). mode = (p, m) selects the LG mode; the formula is p-independent (in EDM's normalization), so only |m| matters in practice.

Formula: a₀² = 2W·|qe|² / (ε₀ c³ A(p,m) (me ω)² w₀² τ₀ √(π/2)) with A(p,m) = (π/2)·(|m|!)².

Use this script-side when u0_constructor (e.g. SVector{8}) precludes LaguerreGaussLaser's built-in pulse_energy initialization, since the init sub-problem inherits the constructor with the wrong dimension. Pass the result to LaguerreGaussLaser(; a0 = …) instead.

See references/lg_pulse_energy_a0.tex for the derivation.

accumulate_field(trajs, screen, alg; mode = Val(:split), solve_kwargs...)

Compute the radiated electromagnetic field on screen from electron trajs.

For each electron and pixel, solves the retarded-time ODE (as in accumulate_potential), builds the Liénard–Wiechert Faraday tensor at each observer-time sample via lienard_wiechert_F_split, and coherently sums the fields over electrons. The far (1/R) and near (1/R²) pieces are accumulated in separate buffers: this avoids the (c² − X𝔞) cancellation at small a₀ and keeps the (near-field-dominated) total from swamping the far-field sum. The Faraday tensor is antisymmetric and the fields are linear in it, so each contribution is stored in the compact (E, B) basis rather than the redundant 16-component matrix (Σᵢ extract_EB(Fᵢ) = extract_EB(Σᵢ Fᵢ)).

Returns (; E, B, E_far, B_far), each Array{Float64,4} of shape (N_samples, 3, Nx, Ny): E, B are the total time-domain fields (far + near, including their cross term) for screen_observables; E_far, B_far are the far field alone, whose observables are the radiated energy/angular momentum exactly (the near field is recoverable as E − E_far).

mode = Val(:total) returns only (; E, B) (the total), a type-stable trim for callers that don't need the split; Val(:split) (the default) keeps all four.

accumulate_field(trajs, screen, ::GPUKernelNewton, backend; n_iters = 2, mode = Val(:split), sync_per_electron = true)

Field counterpart of the GPUKernelNewton accumulate_potential method: per-slot Newton light-cone solve instead of the RK4 retarded-time march, otherwise identical in buffers, mode, and streaming to the GPUKernelRK4 accumulate_field method. timer = LaunchTimer() records a device-event pair per launch (see LaunchTimer). coef_reuse = Val(true) holds each slot's spline-interval coefficients in registers across the per-slot evaluations (bit-identical; trades ~4D live registers for the coefficient loads, see IntervalCoefs).

accumulate_field(trajs, screen, ::GPUKernelRK4, backend; n_substeps = 1, mode = Val(:split), sync_per_electron = true)

GPU counterpart of the CPU accumulate_field: per-electron kernel launch performs the retarded-time RK4 integration and Liénard–Wiechert (E, B) accumulation in one pass, coherently summing over electrons on the device.

Mirrors accumulate_potential(trajs, screen, ::GPUKernelRK4, backend), but uploads the acceleration spline (to_gpu(traj; with_acceleration = true)) since the far field needs 𝔞μ, and accumulates the split (far/near) Liénard–Wiechert field instead of the 4-potential. Returns (; E, B, E_far, B_far), each (N_samples, 3, Nx, Ny) — identical shape to the CPU accumulate_field: E, B total, E_far, B_far the far field alone (see lienard_wiechert_F_split). mode = Val(:total) returns only (; E, B) (a type-stable trim); Val(:split) (the default) keeps all four.

n_substeps and sync_per_electron behave exactly as in the potential kernel; see accumulate_potential and recommended_n_substeps.

buffers / finish drive electron BATCHING: finish = false returns the live FieldAccumulator instead of the downloaded cube, and passing it back as buffers accumulates the next batch of electrons into the same device buffers, so the host only ever holds one batch of trajectory splines. The download happens once, on the call with finish = true (or through finish_field); on one device the batched sum is bit-identical to the single call.

timer = LaunchTimer() records a device-event pair per launch (see LaunchTimer). coef_reuse = Val(true) holds each slot's spline-interval coefficients in registers across the per-slot evaluations (bit-identical; trades ~4D live registers for the coefficient loads, see IntervalCoefs).

accumulate_field_sharded(trajs, screen, alg, backend;
                         devices = 1:gpu_device_count(backend),
                         reduce = :device, reduce_workers = min(4, Threads.nthreads()),
                         kwargs...)
    -> (; E, B[, E_far, B_far])

Shard trajs across devices (vendor-native 1-based ids) and run the single-device accumulate_field on each shard CONCURRENTLY — one Threads.@spawn task per device, each pinned with gpu_device!(backend, d) so its buffers + kernels land on that GPU — then sum the per-device partials into ONE host cube set. The sum is exact by linearity; only its summation order differs from the single-device call (last-bit).

reduce = :device (default) sums on the GPUs: the first device to finish keeps its buffers as the accumulator, the others fold theirs into it over the device-to-device path as they finish, and one download + permute produces the host cube (the permute runs on reduce_workers threads, each holding a 1/16-cube staging slab). reduce = :host is the streamed host reduce: every partial is downloaded, permuted and added on the host under a lock (one cube resident, one single-threaded full-cube permute per device). reduce_stats = ReduceStats() receives the wall-clock time of the folds and of the final download (see ReduceStats).

buffers / finish extend the electron BATCHING of accumulate_field to the sharded path: finish = false shards this batch, accumulates it into each device's own buffers and returns the ShardedFieldAccumulator holding them, which the next batch takes as buffers; the reduce and the download run once, on the call with finish = true. Each batch is split over the same devices, so a device that draws nothing from the last batch still folds in the electrons it accumulated earlier. Batching changes which electrons meet on which device, so the summation order — and with it the last bits of the cube — differs from an unbatched run.

Needs ≥length(devices) Julia threads (julia -t): each per-device task is GPU-bound and blocks its thread on the final device→host copy, so they only overlap on separate OS threads. Each device holds a full prod-size buffer set (see the VRAM budget), so this trades device count for memory, not memory for device count. The same device id may appear more than once (e.g. devices = [1, 1]): the shards then time-share that GPU — pointless for throughput but the exactness check the CPU-backend test relies on.

accumulate_potential(trajs, screen, alg, backend::Backend; solve_kwargs...)

Compute the Liénard-Wiechert 4-potential using CPU retarded-time solve and GPU-accelerated accumulation via AcceleratedKernels.

backend is a KernelAbstractions backend (e.g., CUDA.CUDABackend()). Uses the original trajs (CubicSpline-based) for the CPU retarded-time solve, and converts to GPUCubicSpline internally for the GPU accumulation phase.

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
• 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

accumulate_potential(trajs, screen, ::GPUKernelNewton, backend; n_iters = 2, sync_per_electron = true)

Newton light-cone GPU path: like the GPUKernelRK4 unified kernel, but the retarded proper time at each saveat slot is found by solving the light-cone condition directly (n_iters warm-started Newton corrections per slot) instead of integrating the retarded-time ODE between slots. The condition is evaluated in the screen-relative light-front spelling f = tₖ − ψ(τ) − ρ²/(R + d³) (algebraically identical to x⁰_target − x⁰(τ) − |r_obs − x⃗(τ)| = 0; see the file header), with the sample grid carried as offsets tₖ = x⁰_k − z_screen.

Spline-eval budget per slot is n_iters + 1 (the final residual eval doubles as the accumulation eval), vs 4·n_substeps + 1 for the RK4 march; and the per-slot error does not accumulate along the march, so accuracy is set by the convergence of the last Newton step alone.

sync_per_electron as in the GPUKernelRK4 method.

buffers / finish drive electron BATCHING: finish = false returns the live FieldAccumulator instead of the downloaded cube, and passing it back as buffers accumulates the next batch of electrons into the same device buffers, so the host only ever holds one batch of trajectory splines. The download happens once, on the call with finish = true (or through finish_field); on one device the batched sum is bit-identical to the single call.

timer = LaunchTimer() records a device-event pair per launch (see LaunchTimer). coef_reuse = Val(true) holds each slot's spline-interval coefficients in registers across the per-slot evaluations (bit-identical; trades ~4D live registers for the coefficient loads, see IntervalCoefs).

accumulate_potential(trajs, screen, ::GPUKernelRK4, backend; n_substeps = 1)

Unified GPU path: per-electron kernel launch performs retarded-time RK4 integration and Liénard-Wiechert accumulation in one pass. Compared to the two-phase AcceleratedKernels path:

• no τ_all intermediate buffer (saves ~N_samples × Nx × Ny × 8 bytes)
• no CPU retarded-time solve, so Phase-1 wall time scales with GPU rather than CPU thread count
• accumulation order eliminates the round(Int, …) saveat-slot mapping that the CPU+AK path needs (slot index k is the loop variable).

n_substeps controls how many fixed-step RK4 sub-steps are taken between successive saveat slots (and within the bridge step). Required when ω · δx⁰ lies outside RK4's accurate range. The per-step amplitude error of fixed-step RK4 on a mode of frequency ω is |R(iθ)| − 1 ≈ θ⁶/144 with θ = ω · dt (R the RK4 stability polynomial): at θ = π/2 that is ~7.5%/step, at θ ≈ 0.23 it is ~1e-6. So to hold per-step error ≤ ε, take n_substeps = ⌈ ω·δx⁰ / (144 ε)^(1/6) ⌉. For the relativistic Thomson script (a₀ = 10, δt = T/4 ⇒ ω·δx⁰ ≈ π/2), the relevant ω is the highest harmonic present (≈ a₀³ × ωfundamental for nonlinear Thomson), so `nsubsteps ≈ 8` recovers parity with the adaptive-Tsit5 reference.

sync_per_electron (default true) inserts a KernelAbstractions.synchronize before freeing each trajectory's device buffers — safe but serializes the electron loop. Setting it false drops the per-electron sync and relies on stream-ordered async free (the kernel that still reads the buffers is queued ahead of the free on the same stream), letting electron N+1's upload overlap kernel N. Verified correct on CUDA and ROCm backends.

timer = LaunchTimer() records a device-event pair per launch (see LaunchTimer). coef_reuse = Val(true) holds each slot's spline-interval coefficients in registers across the per-slot evaluations (bit-identical; trades ~4D live registers for the coefficient loads, see IntervalCoefs).

angular_momentum_flux_z(T, r) -> Real

z-component of the radiated angular-momentum flux density crossing the screen, at screen point r = (x, y, z), derived from the stress-energy tensor T^{μν}. Summed over the screen (× dA) and observer-time samples (× dt) it gives the total radiated L_z; divided by the radiated energy it gives the OAM per photon.

Uses the exact Maxwell-stress form x Tᶻʸ − y Tᶻˣ (the z-row of T), the covariant angular-momentum flux M^{z x y} = xᵘ Tᶻᵛ − xᵛ Tᶻᵘ; no far-field approximation. With slots 1=time, 2,3,4 = x,y,z: Tᶻʸ = T[4,3], Tᶻˣ = T[4,2].

canonical_state_order(traj::TrajectoryInterpolant) -> Bool

Whether the state spline stores x⁰…x³ at components 1:4 and u⁰…u³ at 5:8 — the layout the GPU kernels assume (they index the state with literal constants). The sol constructor always produces it; a hand-built spline must be arranged that way before to_gpu.

extract_EB(F, c) -> (E, B)

Recover the electric and magnetic 3-vectors from an upper-index Faraday tensor F^{μν} — the inverse of faraday. With (+,−,−,−) indexing (slot 1 is time): Eⁱ = c·Fⁱ⁰ and B = (F⁴³, F²⁴, F³²). Linear in F, hence commutes with the electron sum.

faraday(E, B, c) -> SMatrix{4,4}

Upper-index Faraday tensor F^{μν} from the field 3-vectors, in the package's (+,−,−,−) convention: F^{0i} = −Eⁱ/c, F^{ij} = −ε^{ijk}Bᵏ. Returns a StaticArrays matrix so the same assembly serves both the symbolic models (with Num entries) and the numeric screen reduction. Its inverse is extract_EB.

finish_field(acc::FieldAccumulator; sink = nothing, workers = 1) -> (; E, B[, E_far, B_far])

Download the accumulated cube from a FieldAccumulator without adding more electrons — the explicit form of accumulate_field(no_more_electrons, …; buffers = acc). workers > 1 spreads the download's permute over that many tasks (see _download_permuted). sink receives the device buffers instead, exactly as in accumulate_field.

The accumulator keeps its buffers (and its sum) afterwards; call it again, or keep accumulating into it, if that is what the driver wants.

flop_profile(alg::GPUKernelNewton; mode = Val(:split), n_iters = 2) -> NamedTuple
flop_profile(alg::GPUKernelRK4;    mode = Val(:split), n_substeps = 1) -> NamedTuple

Algorithmic FLOP cost of the accumulate_field kernel for alg with the given accuracy knob and mode (Val(:split) / Val(:total), or the symbol), measured on the host by running the kernel body on CountedFloats.Counted{Float64} (see the file header). Returns

• flop_per_slot — FLOPs per (electron, pixel, observer sample), and per_slot, its per-category breakdown (add, mul, div, sqrt, fma, pow, trans exact; cmp, other are not FLOPs and are rounded means);
• flop_per_pixel_launch / per_pixel_launch — the per-(electron, pixel) setup cost (window edges, RK4 bridge; for RK4 minus the one inter-slot advance the last slot skips);
• bytes_per_slot — device-buffer read-modify-write traffic per slot (12 or 6 doubles), arithmetic_intensity = flop_per_slot / bytes_per_slot;
• convention, alg, mode, n (the accuracy knob), n_name.

A run's total is flop_per_slot · slots_executed + flop_per_pixel_launch · N·Nx·Ny (slots_executed from window_coverage; N·N_samples·Nx·Ny when the window is fully covered). Costs milliseconds; the GPU is never touched.

harmonic_bins(N_samples, δt, ω, harmonics) -> Vector{Int}

The rfft bin indices closest to the harmonics n·ω of the fundamental, for a time series of N_samples samples spaced δt apart. harmonics is any iterable of integers (e.g. (1, 2)). Deduplicates the locator copied across the solver/plot scripts.

Errors when a requested harmonic lies above the sampling Nyquist 1/(2δt): nearest-match would otherwise silently clamp it onto the last rfft bin, and the caller would publish maps labeled with the requested n while holding that bin's content (the inverse-Thomson ≈4γ²ω aliasing trap).

harmonic_colorrange(data) -> (lo, hi)

Default per-panel color range for harmonic maps: the data extrema, guarded against a degenerate/underflowing panel (falls back to (-1, 1)). See also [symmetric_colorrange].

harmonic_maps(field::NamedTuple, bins) -> Array{ComplexF64,4}   # primary: (;E,B) → 6 components
harmonic_maps(cube::AbstractArray{<:Number,4}, bins) -> Array{ComplexF64,4}   # generic cube

Spatial maps of a field/potential at the harmonic bins (see harmonic_bins): rfft along the time (first) axis and slice out the rows in bins.

The field method takes (; E, B) (each (N_samples, 3, Nx, Ny)) and returns (length(bins), 6, Nx, Ny) — components Eˣ Eʸ Eᶻ Bˣ Bʸ Bᶻ (E in 1:3, B in 4:6). This is the primary entry point for the field runs.

The cube method is the generic core: any (N_samples, n_components, Nx, Ny) array → (length(bins), n_components, Nx, Ny); it also serves the 4-component 4-potential A. Both deduplicate the per-component rfft reduction copied across 5 scripts — the transform is done one component at a time (these cubes are tens of GB at full resolution).

iso_mesh(vol, Xs, Ys, Zs, level) -> GeometryBasics.Mesh | nothing

Marching-cubes isosurface of vol on the grids (Xs, Ys, Zs) (dim 1 ↔ Xs), as a GeometryBasics.Mesh ready for mesh!, or nothing when the level cuts no surface. This is the RPR-compatible replacement for contour!/volume! (unsupported / GPU-segfaulting on RPR backends), and renders identically on GLMakie.

lg_peak_factor(p, m) -> F

Peak transverse field of LaguerreGaussLaser(; a0, radial_index = p, azimuthal_index = m) on its focal plane, relative to the Gaussian mode (p = m = 0) with the same a0 and polarization:

F = √((p+1)_|m|) · max_ρ |(√2ρ)^|m| ₁F₁(−p, |m|+1, 2ρ²) e^{−ρ²}|,   ρ = r/w₀.

EDM's a0 scales the underlying Gaussian amplitude E₀ (the LaserTypes convention, see references/lg_pulse_energy_a0.tex), so the ring of a high-|m| mode is far hotter than a0 suggests: F(0, 7) ≈ 1945.5, while F(2, −2) ≈ 0.98. For p = 0 the ring sits at ρ = √(|m|/2) and F = √(|m|!)·|m|^{|m|/2}·e^{−|m|/2}. The polarization vector multiplies every mode alike, so F does not depend on it; E_z (order 1/(k w₀)) is not included.

lienard_wiechert_F(X, u, 𝔞, K, c) -> SMatrix{4,4}

Total Liénard–Wiechert field-strength (Faraday) tensor F^{μν} = F_near + F_far, the sum of the two pieces from lienard_wiechert_F_split:

F^{μν} = K [ (Xᵘ𝔞ᵛ − Xᵛ𝔞ᵘ) / (Xu)²  +  (c² − X𝔞)(Xᵘuᵛ − Xᵛuᵘ) / (Xu)³ ].

Indices are upper, matching faraday and the Lorentz force in ChargedParticle. Use lienard_wiechert_F_split when accumulating, to keep the far and near fields in separate coherent sums (see accumulate_field).

lienard_wiechert_F_split(X, u, 𝔞, K, c) -> (F_near, F_far)

Liénard–Wiechert field-strength (Faraday) tensor split by power of R into its near-field (1/R²) and far-field (1/R) pieces — the form used for numerically clean, separately-coherent accumulation in accumulate_field. With Xu = m_dot(X,u) ∝ R and X𝔞 = m_dot(X,𝔞) ∝ R:

F_near = K c² (Xᵘuᵛ − Xᵛuᵘ) / (Xu)³                            ∝ 1/R²
F_far = K [ (Xᵘ𝔞ᵛ − Xᵛ𝔞ᵘ) / (Xu)²
          − X𝔞 (Xᵘuᵛ − Xᵛuᵘ) / (Xu)³ ]                        ∝ 1/R

X is the retarded displacement 4-vector Xᵘ = xᵘ_obs − xᵘ_src(τᵣ) (= (R, R·n̂)), u the 4-velocity and 𝔞 the 4-acceleration at τᵣ; K = q/(4πε₀c).

The subtlety: the −X𝔞 term rides the velocity bivector (Xᵘuᵛ − Xᵛuᵘ) yet scales as 1/R (because X𝔞 ∝ R), so it belongs to the far field. lienard_wiechert_F groups it with the c² term as (c² − X𝔞), which loses ≈log₁₀(c²/X𝔞) significant digits when X𝔞 ≪ c² (small a₀); computing F_far directly never forms that difference. F_near + F_far is identically lienard_wiechert_F in exact arithmetic.

observer_window_start(τi, Z, half_width, Rmax; c)

Observer time x⁰ at which the LAST retarded image of proper time τi reaches a square screen of half-width half_width at distance Z from a source disc of radius Rmax centred on the axis, for electrons at rest at τi: the corner pixel (√2·halfwidth from the axis) and the electron on the far rim, `c·τi + hypot(Z, √2·halfwidth + Rmax). Opening the observer window here makes its first sample already see every electron at every pixel, so no sample is partially populated. (The edge-midpoint distancehypot(Z, half_width + Rmax)` is NOT a bound: the pixels outside the screen's inscribed circle then miss the far-side electrons' first samples.)

phase_winding_fit(az, phase; weights = nothing) -> (; slope, intercept, unwrapped)

Least-squares fit of the azimuthal phase winding: model phase ≈ slope·az + intercept, with az a vector of azimuths sorted ascending (as from ring_pixels) and phase the wrapped phase angle(F) at those points.

Because angle wraps to (-π, π], a winding field reads as a sawtooth, so a naive fit on the raw series returns slope ≈ 0 regardless of the true winding. The phase is therefore first unwrapped along az (adding ±2π at jumps larger than π) before fitting — then the slope recovers the winding ℓ = d∠F/dφ and intercept the absolute phase offset b. weights (e.g. abs.(F) on the ring) down-weight phase nodes where the angle is ill-defined; nothing ⇒ uniform. Returns the fitted slope, intercept, and the unwrapped series (for plotting).

plot_harmonic_grid(maps, x_grid, y_grid; w₀=1, labels, title="", colormap=:jet,
                   colorrange=harmonic_colorrange, transform=real, colorbar_offset=true,
                   ncols=3, panelsize=300, outfile=nothing)

Draw a grid of per-component heat-maps for one harmonic — the real part of each component of maps, which is (n_components, Nx, Ny) (e.g. a slice fields_h[k, :, :, :]). Axes are shown in units of w₀ (set it to the beam waist; left as 1 ⇒ raw coordinates, label "x"/"y"). colormap and colorrange are configurable: colorrange is any data -> (lo, hi) applied per panel (default [harmonic_colorrange] = guarded extrema; pass [symmetric_colorrange] for a diverging map). With colorbar_offset=true (default) each colorbar gets explicit lo/mid/hi ticks — this silences the PlotUtils "No strict ticks found" warning on the radiated field's tiny ~1e-17 ranges, and lets Makie's native formatter render the labels as m×10ⁿ (RichText superscripts, since the range spans >4 orders); pass false for O(1) data such as the phase grid (keeps Makie's automatic ticks). ncols sets panels-per-row (3 ⇒ 2×3 for E/B, 2 ⇒ 2×2 for the 4-potential). Saves to outfile if given; returns the Figure.

Requires CairoMakie — using CairoMakie activates the implementation (package extension); calling this without it raises a MethodError.

plot_phase_grid(maps, x_grid, y_grid; w₀=1, labels, title="", ncols=3, panelsize=300, outfile=nothing)

Like [plot_harmonic_grid] but plots the phase angle.(maps[c]) of each component on a cyclic :phase colormap over (-π, π) — the (x/w₀, y/w₀, ∠F) view. The panel title still reports each component's peak complex amplitude. Requires CairoMakie (package extension).

plot_phase_polar(maps, x_grid, y_grid; w₀=1, labels, radii, tol, title="", panelsize=300, outfile=nothing)

Polar companion to [plot_phase_with_rings]: for each component of maps ((3, Nx, Ny)), a PolarAxis with angular = azimuth φ and radial = ∠F shifted to [0, 2π) (angle(F)+π), one colour-matched series per test ring R ∈ radii (same pixels as the cartesian view, via ring_pixels). A winding of charge ℓ traces ℓ radial oscillations as φ runs once around. A separate diagnostic from the main phase figure. Requires CairoMakie (package extension).

plot_phase_with_rings(maps, x_grid, y_grid; w₀=1, labels, radii, tol, title="",
                      panelsize=300, outfile=nothing)

Combined per-field phase view for one field type (maps is (3, Nx, Ny) — e.g. the E or B slice fields_h[k, 1:3, :, :] / [k, 4:6, :, :]). One row per component, two columns:

• left — the phase heatmap angle.(maps[c]) on the cyclic :phase colormap over (-π, π), with the test annuli drawn as dashed circles at R ± tol (in w₀ units, to match the axes), one colour per radius;
• right — the azimuthal phase winding: angle(F) of the pixels on each test circle of radius R ∈ radii (a thin annulus, half-width tol, about the grid centre) scattered against azimuth atan(y, x), coloured to match the circles on the left. A vortex of topological charge ℓ shows ℓ phase windings as φ runs once around. Every component gets the phase_winding_fit line overlaid (re-wrapped to ±π), with slope ≈ ℓ and intercept b the phase offset (the longitudinal Eᶻ/Bᶻ wind too — e.g. the ℓ≈3 vortex).

Returns (; fig, fits) where fits[c] (one per component) is (; slope, b) — vectors over radii (NaN where a ring has no pixels), so callers can record the winding/offset in [plot_params]. Requires CairoMakie (package extension).

plot_power_spectrum(freqs, power_spec; ω, labels, marks=nothing, n0=1, colors=…, linestyles=nothing, title="", outfile=nothing)

Log-y plot of per-component power spectra power_spec ((Nf, n), from [power_spectrum]) vs freqs, x-axis in units of the run's radiated fundamental: ω₁ (= ω/2π) when n0 = 1, else the backscattered line ωbs = n0·ω₁ (boosted runs; n0 ≈ 4γ² from the manifest's `backscattern0). Reference gridlines come from Makie's automatic ticks;marks(extracted harmonic bins, in ω₁ units) adds explicit guides. Saves tooutfileif given; returns theFigure`.

Requires CairoMakie — using CairoMakie activates the implementation (package extension).

power_spectrum(cube) -> Matrix{Float64}

Per-component frequency power spectrum of a (N_samples, n_components, Nx, Ny) cube, summed over the screen: out[ν, c] = Σ_pixels |rfft(cube[:,c,:,:], 1)[ν]|². Returns (N÷2+1, n_components); pair with rfftfreq(N_samples, 1/δt) for the frequency axis. One component at a time (memory).

recommended_n_substeps(screen, c; ω_max = π * c / step(screen.x⁰_samples), rtol = 1e-6)

Suggest an n_substeps value for the fixed-step GPU kernels from the screen's own sampling, with an optional override for the highest relevant frequency.

c is the speed of light in the working units (e.g. getdefault(world.c)). The screen stores the time-like coordinate x⁰ = c·t, so δx⁰ = step(x⁰_samples) is a length and the integration time step is δt = δx⁰ / c.

Fixed-step RK4 stepping the retarded-time ODE dτ_r/dx⁰ = f multiplies a mode of angular frequency ω by the stability polynomial |R(iθ)|, θ = ω·δt. Expanding, |R(iθ)| − 1 ≈ θ⁶/144, so at θ = π/2 the per-step amplitude loss is ~7.5%, while at θ ≈ 0.23 it is ~1e-6. Sub-stepping each saveat interval into n pieces drives θ = ω_max·δt/n down into the accurate range.

ω_max is a temporal angular frequency (rad/time) and defaults to the screen's Nyquist angular frequency π·c/δx⁰ = π/δt: a uniform x⁰_samples grid cannot represent anything faster, so resolving the RK4 step to Nyquist guarantees the integration is as accurate as the output grid can carry. With this default c and δx⁰ cancel and the result is a screen-independent constant (≈14 at rtol = 1e-6). A caller who knows their top harmonic sits below Nyquist (e.g. ω_max = n_harmonics · 2π·c/λ) can pass it explicitly to take fewer sub-steps.

rtol is the target per-step amplitude error.

ring_pixels(x_grid, y_grid, R; tol) -> (idxs, az)

Grid pixels whose distance from the centre lies within tol of radius R (a thin annulus about the grid origin), as a Vector{CartesianIndex{2}} sorted by azimuth, together with those azimuths az = atan(y, x) ∈ (-π, π]. x_grid/y_grid must be centred on the beam axis. Reads a field/phase on a test circle for phase_winding_fit; shared by the phase-ring plots and the LPWA-vs-numeric comparison so both sample identical pixels.

rpr_capture(screen; hybrid, exposure = 0.2, sat = 1.0, dump_raw = "") -> Matrix{RGBf}

Capture a rendered RPRMakie Screen to an image. On Northstar (hybrid = false) this is plain colorbuffer. On the Hybrid plugins (hybrid = true), where colorbuffer's resolve returns black, it reads the raw HDR framebuffer, normalizes by per-pixel sample count (alpha), applies an exposure-scaled LUMINANCE Reinhard (per-channel compression drags saturated highlights to white) + sRGB gamma, then an optional post-tonemap saturation sat in gamma space (≈1.25 recovers Northstar's per-channel photolinear vibrancy without blowing highlights). dump_raw serializes the raw HDR buffer for offline regrading.

rpr_enable_multiframe!()

Work around the HybridPro rprSceneClear segfault when rendering many frames in one process: RPRMakie's singleton Context releases the previous context on every new Screen, and HybridPro crashes in the release. After this call, new contexts are created non-singleton — they LEAK (measured ~1.7 GB VRAM/frame on heavy scenes) but never tear down. Bound the leak by chunking long frame ranges across processes (~16 frames/process on a 48 GB card).

rpr_tune!(screen; quality = nothing, denoiser = nothing, ray_depth = nothing)

Set HybridPro context parameters on an RPRMakie Screen (call right after constructing it; Hybrid plugins only — Northstar uses different knobs):

• quality: "low" | "medium" | "high" | "ultra" — RPRCONTEXTRENDER_QUALITY (0x1001), gates HybridPro's algorithmic shortcuts.
• denoiser: "none" | "svgf" | "asvgf" | "ml" — RPRCONTEXTPT_DENOISER (0x102D; ml needs RadeonImageFilters libs the jll doesn't ship).
• ray_depth: Int — maxrecursion + refraction/glossyrefraction depths (glossy capped at 8). Nested translucent shells want 12-16: rays that exhaust the preset budget terminate BLACK. Probed: HybridPro ACCEPTS the recursion/diffuse/glossy/refraction depth family and REJECTS shadow depth, Russian roulette, sampler type, and adaptive sampling.

screen_observables(field, screen; ε₀) -> NamedTuple

Derive radiation diagnostics from the accumulated field = (; E, B) (the output of [accumulate_field]) as functions of the observer-time sample index — all following from the Faraday tensor through the electromagnetic stress-energy tensor Tᵘᵛ:

• S — Poynting vector, (N, 3, Nx, Ny)
• energy_density — u = ½ε₀(E² + c²B²), (N, Nx, Ny)
• Lz_density — z angular-momentum flux dens, (N, Nx, Ny)
• energy_total — ∫∫ Sᶻ dA dt (energy through the screen)
• Lz_total — ∫∫ Lz_density dA dt

Each pixel's field is reassembled into the Faraday tensor faraday, then T^{μν} is formed covariantly via stress_energy; S, energy density, and the angular-momentum flux angular_momentum_flux_z are read off T.

screen_spectrum(field, screen; ε₀, bins = nothing) -> NamedTuple

Frequency-domain radiation diagnostics from the accumulated field = (; E, B).

The FFT is linear, so it acts on the fields E, B (not on the quadratic observables); every component of the stress-energy tensor then has a frequency-domain image T̃(ω) that is a bilinear cross-spectral density Re[F̃*(ω) F̃(ω)] of the transformed field — assembled here per kept bin:

• freqs — frequencies (1/time units) at the kept bins
• E_ω, B_ω — complex field spectra, (n_bins, 3, Nx, Ny)
• energy_ω — ½ε₀(|Ẽ|² + c²|B̃|²), (n_bins, Nx, Ny)
• S_ω — Re[Ẽ*×B̃] / μ₀, (n_bins, 3, Nx, Ny)
• Lz_ω — x T̃ᶻʸ − y T̃ᶻˣ, (n_bins, Nx, Ny)

bins selects which rfft frequency indices to keep (e.g. harmonic bins located via rfftfreq); nothing keeps all. Pass a small bins when only a few frequencies are needed — the complex E_ω, B_ω are full grid-sized per bin, so keeping every bin materialises arrays as large as E itself.

The transform is blocked over screen rows, so the transient complex array is (N÷2+1, 3, Nx), never the full grid. Spectra use the δt-scaled rfft (Ẽ ≈ ∫E e^{iωt} dt); the phase is independent of this scaling, so a phase map of a component at bin m is simply angle.(E_ω[m, i, :, :]).

stress_energy(F, g, μ₀) -> SMatrix{4,4}

Electromagnetic stress-energy tensor T^{μν} from the upper-index Faraday tensor F, metric g, and vacuum permeability μ₀, in the (+,−,−,−) convention (F^{0i} = −Eⁱ/c):

M = F g Fᵀ                       # M^{μν} = F^{μα} g_{αβ} F^{νβ}
F_{αβ}F^{αβ} = tr(g M)           # = 2(B² − E²/c²)   (uses Fᵀ = −F)
T^{μν} = (1/μ₀) [ −M^{μν} + ¼ g^{μν} (F_{αβ}F^{αβ}) ]

The sign is fixed so T⁰⁰ = ½ε₀(E²+c²B²) > 0 (energy density) and Sⁱ = c·T⁰ⁱ (Poynting). Written with StaticArrays matrix ops so one function serves both the symbolic EMFieldDynamics (via @SMatrix of Num) and the numeric reduction in screen_observables.

References: Landau & Lifshitz "Classical Theory of Fields" §33; Jackson §12.10.

sunflower(n, α) -> Vector{Vector{Float64}}

n points on the unit disc in Vogel's sunflower spiral: point k at angle k·2π/ϕ² (golden-angle stride, ϕ = (1 + √5)/2) and radius √(k − ½)/√(n − (b + 1)/2), a uniform areal density. The outermost b = round(Int, α√n) points sit exactly ON the boundary ρ = 1 (α sets how many; the production layout uses α = 2, i.e. 2√n edge points), so tests like "ρ > R" flag them from the start unless excluded.

Scale by the disc radius: Rmax * sunflower(N, 2). This is the electron layout of the production solvers (thomsonscattering.jl, inversethomson_scattering.jl, lpwa.jl); analysis scripts that rebuild initial conditions must use it to reproduce them bit for bit.

symmetric_colorrange(data) -> (-m, m)

A color range symmetric about zero (m = maximum(abs, data)), for diverging colormaps like :seismic; guarded against a degenerate panel. Pass as plot_harmonic_grid's colorrange.

trajectory_interpolants(ensemble_sol::EnsembleSolution)

Return a vector of TrajectoryInterpolant objects, which interpolate the trajectories in a fast and thread safe manner.

trajectory_span_for_window(τi, τf, x⁰_samples, Z, half_width, Rmax; c, margin_samples = 2, stretch = 1)

Proper-time span (τ_lo, τ_hi) ⊇ (τi, τf) a trajectory solve must cover so that every pixel of the screen sees every electron at EVERY sample of x⁰_samples, i.e. the field kernels' strict- interior slot range k_start:k_end is 1:N_samples for every (electron, pixel) pair. The history must begin margin_samples before the window opens at the latest-arriving pair (the corner pixel / far-rim electron, see observer_window_start) and must last until the EARLIEST retarded image of its end — the axis pixel, where the arrival is the light-front coordinate x⁰ − x³ plus Z — reaches the last sample. stretch is the proper time per unit of light-front advance: 1 for electrons at rest, γ(1+β) for force-free flight toward the screen (x⁰ − x³ = c·τ/(γ(1+β))). The physics span (τi, τf) is never shrunk; the extension only adds history where the electron moves freely.

window_coverage(trajs, screen; recount = true) -> NamedTuple

Check on the host that the observer window screen.x⁰_samples lies inside every pixel's arrival window for every electron — the condition under which the GPU field kernels (accumulate_field with GPUKernelNewton / GPUKernelRK4) execute all N_samples slots per (electron, pixel), the cube receives every electron's full history, and the nominal work N·N_samples·Nx·Ny is exact. Uses the kernels' own window arithmetic (_window_edge, same floor/ceil slot formulas) on the host spline, evaluated at the four corner pixels and the pixel nearest to the electron, which bound the arrival offsets exactly (see the file header); clipped electrons are recounted over every pixel when recount = true.

Fields of the result:

• ok::Bool — every electron fully covered.
• slots_nominal, slots_executed, slot_fill = slots_executed / slots_nominal (slots_executed is missing for a clipped run with recount = false).
• electrons_clipped, slots_dropped = slots_nominal - slots_executed.
• lead_margin_samples — samples the window could start earlier before the first pixel is clipped at the start (negative = already clipped by that many samples at the worst pixel); tail_margin_samples — same for the end of the window (a negative value means some trajectory ends before the window does, at some pixel).
• worst_electron — index with the smallest tail margin.
• per_electron::Vector of (; k_start_max, k_end_min, lead_margin, tail_margin, slots_executed).

Cost: ten spline evaluations per electron (milliseconds for 10⁴ electrons), plus 2·Nx·Ny per clipped electron. Rounding note: floor/ceil at an exact sample boundary can differ by one slot between this host evaluation and a device launch; that cannot flip ok except at a measure-zero boundary.