Internal

Public APIs

TestParticle.AdaptiveHybrid Type

julia

AdaptiveHybrid{T}

Hybrid algorithm that switches between adaptive Boris (FO) and RK45 (GC).

Mode switching uses hysteresis with two adiabaticity thresholds to avoid rapid GC ↔ FO chattering near a single bound:

threshold_gc_to_fo (α): switch GC → FO when ε ≥ α.
threshold_fo_to_gc (β): switch FO → GC when ε < β.

With α ≥ β the band β ≤ ε < α forms a buffer where the current mode is retained, so a particle hovering near the boundary does not oscillate.

The adiabaticity parameter ε compared against the thresholds is selected by the adiabaticity field, which may be:

:curvature (default) → ε_curv = ρ_L / R_c (curvature drift; reproduces the legacy get_adiabaticity behavior exactly),
:gradB → ε_gradB = ρ_L / L_B, with L_B = |B| / |∇B|,
:both → OR of the two criteria: the solver switches to the full orbit whenever either ε_curv ≥ α or ε_gradB ≥ α (equivalently max(ε_curv, ε_gradB) ≥ α),
:jacobian → ε_jac = ρ_L · ‖JB‖_F / |B|, CHIMP's original all-in-one criterion (Bates & Fuentes, 2020), where JB is the Jacobian of B and ‖JB‖_F its Frobenius norm. It captures curvature, grad-B, torsion and shear at once, so it switches to full orbit in non-adiabatic regions that :curvature/:gradB miss (e.g. a field whose direction rotates in space while keeping constant magnitude).

Each check point records only the adiabaticity value used by the selected adiabaticity mode (a single scalar per check) together with the active mode in sol.stats.adiabaticity. This diagnostic is only saved when save_adiabaticity = true (the default); set it to false to skip the buffers and reduce memory/allocation overhead.

TestParticle.Field Type

julia

Field{itd, F} <: AbstractField{itd}

A representation of a field function f, defined by:

time-independent field

F = F (x),

time-dependent field

F = F (x, t) .

Arguments

field_function::Function: the function of field.
itd::Bool: whether the field function is time dependent.
F: the type of field_function.

TestParticle.GPUGrid1D Type

julia

GPUGrid1D{T, B, V, A<:AbstractArray{V, 1}} <: AbstractFieldInterpolator

Uniform 1D Cartesian grid interpolator compatible with GPU device execution.

TestParticle.GPUGrid2D Type

julia

GPUGrid2D{T, B, V, A<:AbstractArray{V, 2}} <: AbstractFieldInterpolator

Uniform 2D Cartesian grid interpolator compatible with GPU device execution.

TestParticle.GPUGrid3D Type

julia

GPUGrid3D{T, B, V, A<:AbstractArray{V, 3}} <: AbstractFieldInterpolator

Uniform 3D Cartesian grid interpolator compatible with GPU device execution.

TestParticle.GPUNonUniformAxis Type

julia

GPUNonUniformAxis{T, V<:AbstractVector{T}, B}

Non-uniform 1D grid axis for GPU-compatible interpolators.

TestParticle.GPUSphericalGrid Type

julia

GPUSphericalGrid{T, Ar, Aθ, Aϕ, V, A, B} <: AbstractFieldInterpolator

Spherical grid interpolator compatible with GPU device execution. Supports both uniform and non-uniform grids in r, θ, and ϕ.

TestParticle.GPUUniformAxis Type

julia

GPUUniformAxis{T, B}

Uniform 1D grid axis for GPU-compatible interpolators.

TestParticle.LazyTimeInterpolator Type

julia

LazyTimeInterpolator{T, F, L}

A callable struct for handling time-dependent fields with lazy loading and linear time interpolation.

Fields

times::Vector{T}: Sorted vector of time points.
loader::L: Function i -> field that loads the field at index i.
buffer::Dict{Int, F}: Cache for loaded fields.
lock::ReentrantLock: Lock for thread safety.

TestParticle.PotentialField Type

julia

PotentialField{TP, TA, TGP, TGA}

A representation of scalar potential phi and vector potential A, along with optional spatial gradient functions grad_phi and grad_A.

TestParticle.PotentialField Method

julia

PotentialField(phi, A; grad_phi = nothing, grad_A = nothing)
PotentialField(A; phi = ZeroField(), grad_phi = nothing, grad_A = nothing)

Construct a PotentialField for canonical tracing. When grad_phi or grad_A are omitted, they are computed automatically via ForwardDiff.

TestParticle.Species Type

Type for the particles: Proton, Electron.

TestParticle.TraceCanonicalProblem Type

julia

TraceCanonicalProblem{uType, tType, isinplace, P, F <: AbstractODEFunction, PF}

Problem specification for particle tracing in canonical coordinates (x, p).

TestParticle.TraceGCProblem Type

julia

TraceGCProblem{uType, tType, isinplace, P, F <: AbstractODEFunction, PF} <: AbstractODEProblem{uType, tType, isinplace}

Problem type for tracing guiding centers.

TestParticle.TraceHybridProblem Type

julia

TraceHybridProblem{uType, tType, isinplace, P, F <: AbstractODEFunction, PF}

Problem type for hybrid GC-FO tracing. Initial condition u0 should be the full orbit state (6-element).

TestParticle.TraceProblem Method

julia

TraceProblem(u0, tspan, p; prob_func = DEFAULT_PROB_FUNC)
TraceProblem(f, u0, tspan, p; prob_func = DEFAULT_PROB_FUNC)

Problem specification for particle tracing. Supports both Boris solvers and standard ODE solvers (e.g. Tsit5, Vern7). When f is not provided, defaults to trace! (in-place) for mutable arrays or trace (out-of-place) for StaticArrays.

CommonSolve.solve Method

julia

solve(prob::TraceProblem, alg::AbstractBoris, ensemblealg; kwargs...)

Trace trajectories particles with a Boris method and return an EnsembleSolution. A shorthand for

julia

solve(EnsembleProblem(prob; prob_func = prob.prob_func), alg, ensemblealg;
    trajectories, kwargs...)

which is the form to use for anything it does not cover, such as an output_func, a reduction, or a prob_func other than the one carried by prob.

Keywords

trajectories::Int=1: number of particles.
safetycopy::Bool=false: whether to deepcopy prob per trajectory.
seed: master seed. SciML derives one reproducible generator per trajectory from it and hands it to prob_func as ctx.rng.
batch_size, pmap_batch_size: SciML's own ensemble keywords, controlling how trajectories are grouped before the reduction and how many are sent to a worker at a time.
Every keyword of the single trajectory solve above.

CommonSolve.solve Method

julia

solve(prob::TraceProblem, alg::AbstractBoris; kwargs...)

Trace one particle with a Boris method through the SciML loop and return an ODESolution. Trace several by wrapping the problem in a SciML EnsembleProblem, which solve(prob, alg, ensemblealg) does for convenience. trajectories belongs to the ensemble; it is not a keyword here.

Keywords

dt: time step. Optional for adaptive methods, which otherwise start from safety * 2π / |q B / m|.
saveat: times to save at, as a collection or as an interval. The state is interpolated linearly inside a step, the dense output these methods admit.
isoutside: boundary check (u, p, t); the trace terminates when it holds.
save_start::Bool=true, save_end::Bool=true, save_everystep::Bool=true. With saveat, save_start and save_end add the ends of the time span to the requested times.
save_fields::Bool=false: append E and B to every saved state.
save_work::Bool=false: append the work rates to every saved state.
maxiters::Int=1_000_000: maximum number of steps.
rng, seed: accepted and ignored; a SciML ensemble passes them to every trajectory it solves.

CommonSolve.solve Method

julia

solve(prob::TraceGCProblem; trajectories::Int=1, dt::AbstractFloat,
    saveat=nothing, isoutside::Function=ODE_DEFAULT_ISOUTOFDOMAIN,
    alg::Symbol=:rk4, abstol=1e-6, reltol=1e-6, maxiters=10000,
    save_fields::Bool=false, save_work::Bool=false)

Trace guiding centers using the RK4 method with specified prob. If alg is :rk45, uses adaptive time stepping.

saveat names the times at which to report the state, as a collection or as an interval; the state is interpolated linearly between the two steps that bracket each time, leaving the trajectory unchanged.

TestParticle.BiKappa Function

julia

 BiKappa(args...; kw...)

Construct a BiKappa distribution. Forwards to VelocityDistributionFunctions.BiKappa.

julia

 BiKappa(B, u0, ppar, pperp, n, kappa; m=mᵢ)

Construct a BiKappa distribution with magnetic field B, bulk velocity u0, parallel thermal pressure ppar, perpendicular thermal pressure pperp, number density n, and spectral index kappa in SI units. The default particle is proton.

TestParticle.BiMaxwellian Function

julia

 BiMaxwellian(args...; kw...)

Construct a BiMaxwellian distribution. Forwards to VelocityDistributionFunctions.BiMaxwellian.

julia

 BiMaxwellian(B, u0, ppar, pperp, n; m=mᵢ)

Construct a BiMaxwellian distribution with magnetic field B, bulk velocity u0, parallel thermal pressure ppar, perpendicular thermal pressure pperp, and number density n in SI units. The default particle is proton.

TestParticle.Kappa Function

julia

 Kappa(args...; kw...)

Construct a Kappa distribution. Forwards to VelocityDistributionFunctions.Kappa.

julia

 Kappa(u0, p, n, kappa; m=mᵢ)

Construct a Kappa distribution with bulk velocity u0, thermal pressure p, number density n, and spectral index kappa in SI units. The default particle is proton.

TestParticle.Maxwellian Function

julia

 Maxwellian(args...; kw...)

Construct a Maxwellian distribution. Forwards to VelocityDistributionFunctions.Maxwellian.

julia

 Maxwellian(u0, p, n; m=mᵢ)

Construct a Maxwellian distribution with bulk velocity u0, thermal pressure p, and number density n in SI units. The default particle is proton.

TestParticle.TerminateOutside Method

julia

TerminateOutside(isoutofdomain)

Create a DiscreteCallback that terminates the simulation if isoutofdomain(u, p, t) is true. The condition follows the same (u, p, t) convention as the isoutside keyword, and is adapted here to the (u, t, integrator) signature that SciML callbacks use.

TestParticle.TraceFieldlineProblem Method

julia

TraceFieldlineProblem(u0, p, tspan::Tuple{<:Real, <:Real}; mode::Symbol=:both)

Helper function to create ODEProblems for tracing magnetic field lines.

Arguments

u0: Initial position.
p: Parameter tuple containing the magnetic field, or a function (analytical/interpolator).
tspan: Time span (arc length) for tracing. typically (0.0, L).
mode: Tracing mode, one of :forward, :backward, :both.

Returns

If mode is :forward or :backward, returns a single ODEProblem.
If mode is :both, returns a Tuple of two ODEProblems (forward and backward).

TestParticle.adiabaticity_components Function

julia

adiabaticity_components(x, Bfunc, q, m, μ, t = 0.0)
adiabaticity_components(x, Bfunc, μ, t = 0.0; species = Proton)

Return the tuple (ε_curv, ε_gradB, ε_jac) of first-order adiabaticity parameters at position x and time t:

ε_curv = ρ / R_c (curvature drift; the definition used by get_adiabaticity),
ε_gradB = ρ / L_B (grad-B drift; CHIMP's eGC), with L_B = |B| / |∇B| the magnetic-field length scale,
ε_jac = ρ · ‖JB‖_F / |B| (CHIMP's all-in-one criterion; JB is the Jacobian of B and ‖JB‖_F its Frobenius norm). Because ‖JB‖_F² = |∇B|² + |B|²|κ|² + |B|²τ² + shear, this folds curvature, grad-B, torsion and shear into one number, so it switches to full orbit in non-adiabatic regions that ε_curv/ε_gradB miss (e.g. a field whose direction rotates in space with constant magnitude). It always satisfies ε_jac ≥ max(ε_curv, ε_gradB).

ρ = √(2 μ m / |B|) / |q| is the gyroradius. When |B| = 0 all three components are Inf.

TestParticle.analytic_projection Method

julia

analytic_projection(f, k, vi, vj, vk, dvk) -> Matrix

Analytic 2-D projection f(vᵢ, vⱼ) = ∫ f dvₖ in [s²/km⁵] of a 3-D phase-space density, evaluated on the velocity grids vi and vj in [km/s] while the k-th axis is integrated over vk in [km/s] with the midpoint rule at spacing dvk. The two remaining axes follow from k as i < j: k = 3 projects onto Vx–Vy, k = 2 onto Vx–Vz and k = 1 onto Vy–Vz.

f takes a velocity in [m/s] and returns f in [s³/m⁶].

TestParticle.bin_centers Method

julia

bin_centers(edges)

Centers of the 1-D bins delimited by edges, in the same unit as edges.

TestParticle.bin_velocity_space Method

julia

bin_velocity_space(vs, weights, edges; vx_source = nothing, average = false)

Accumulate velocity-space samples into a 3-D grid.

vs are sample velocities in [m/s] as returned by get_particle_crossings and weights the weight of each sample. edges are the bin edges in [km/s], shared by all three velocity axes. Samples outside the binned range are discarded, and the returned array carries the unit of weights. Empty bins are zero.

Two estimators are available, and they differ in how the samples that land in a bin are combined:

average = false (default) sums the weights, which is the histogram estimator of the bin-averaged f: feeding f·ΔV/dv³ with f in [s³/km⁶] returns f in [s³/km⁶], and feeding the constant n/(N·dv³) turns binned counts into a Monte-Carlo estimate of f. Its error carries the counting noise ∝ 1/√n of the samples per bin.
average = true takes the mean of the weights, a ratio estimator of the local f that divides out the counting noise. It is the better choice when the sampling spreads only a handful of samples over each bin, e.g. forward Liouville tracking of a uniformly sampled velocity sphere onto a fine detector grid.

Either way project_vdf turns the result into 2-D projections in [s²/km⁵].

If vx_source is given, each sample is rescaled by |vx_source| / |vx|. A trajectory maps a source velocity-space volume onto a detector volume that differs by exactly that ratio, so the rescale turns a source-referred weight into the detector-referred one. This is what makes a flux-weighted launch ensemble consistent with the crossing flux recorded at the detector.

TestParticle.build_interpolator Function

julia

build_interpolator(gridtype, A, grids..., order::Int=1, bc=FillExtrap(NaN);
    coeffs=OnTheFly(), store=StorePolicy())
build_interpolator(A, grids..., order::Int=1, bc=FillExtrap(NaN);
    coeffs=OnTheFly(), store=StorePolicy())

Return a function for interpolating field array A on the given grids.

Arguments

gridtype: CartesianGrid, RectilinearGrid or StructuredGrid. Usually determined by the number of grids.
A: field array. For vector fields, the first dimension should be 3 unless the array contains SVectors.
order::Int=1: order of interpolation in [0,1,3].
bc=FillExtrap(NaN): boundary condition type from FastInterpolations.jl.
- FillExtrap(NaN): Fill with NaN (default).
- ClampExtrap(): Clamp (flat extrapolation).
- WrapExtrap(): Exclusive periodic wrapping ( $L = N Δ x$ ).
coeffs=OnTheFly(): coefficient strategy for cubic interpolation (order=3).
store=StorePolicy(): storage policy forwarded to FastInterpolations.jl. The default copies the grid and field data. Use StorePolicy(copy=false) to alias compatible inputs and avoid the construction-time data copy.

Notes

With StorePolicy(copy=false), the interpolator may alias A and the grid arrays. Do not mutate or resize aliased inputs while the interpolator is in use. When input conversion is required, FastInterpolations.jl may still copy the data.

TestParticle.cart2sph Method

julia

cart2sph(x, y, z)

Convert from Cartesian to spherical coordinates vector.

TestParticle.cart2sphvec Method

julia

cart2sphvec(vx, vy, vz, x, y, z)

Convert a vector from Cartesian to spherical coordinates at the Cartesian position (x, y, z).

TestParticle.embed_vdf Method

julia

embed_vdf(g1, g2, full_centers, M) -> Matrix

Embed a 2-D sub-grid projection M onto a canvas defined by full_centers. g1 and g2 are bin centers along the two axes in [km/s], and full_centers is the target full grid of bin centers.

TestParticle.energy2velocity Method

Return velocity magnitude from energy in [eV].

TestParticle.full_to_gc Function

julia

full_to_gc(xu, param)

Convert full particle state xu to guiding center state state_gc and magnetic moment μ. Returns (state_gc, μ), where state_gc = [R..., vpar].

TestParticle.gc_to_full Function

julia

gc_to_full(state_gc, param, μ, phase=0)

Convert guiding center state state_gc to full particle state xu. Returns xu = [x, y, z, vx, vy, vz].

TestParticle.generate_sphere Function

julia

generate_sphere(nθ=64, nϕ=64, radius=1.0)

Generate a sphere with radius and matching grid points nθ and nϕ. Returns a tuple of (x, y, z) matrices.

TestParticle.get_B_parameters Method

julia

get_B_parameters(x, t, Bfunc)

Evaluate the magnetic field Bfunc at (x, t) and return the field vector B, its magnitude Bmag, the unit vector b̂ = B / Bmag, the gradient ∇B, and the Jacobian JB of the field.

TestParticle.get_adiabaticity Function

julia

get_adiabaticity(r, Bfunc, q, m, μ, t=0.0)
get_adiabaticity(r, Bfunc, μ, t=0.0; species = Proton)

Calculate the adiabaticity parameter ϵ = ρ / Rc at position r and time t. ρ is the gyroradius and Rc is the radius of curvature of the magnetic field.

TestParticle.get_adiabaticity Method

julia

get_adiabaticity(sol::AbstractODESolution, t)

Calculate the adiabaticity parameter ϵ from a solution sol at time t. t must be a scalar.

TestParticle.get_adiabaticity Method

julia

get_adiabaticity(sol::AbstractODESolution)

Calculate the adiabaticity parameter ϵ from a solution sol at all time steps. The parameters q, m, Bfunc are extracted from sol.prob.p.

TestParticle.get_cell_centers Function

julia

get_cell_centers(grid::Union{CartesianGrid, RectilinearGrid})

Return the cell center coordinates of a Meshes.jl grid, one range per dimension. Requires Meshes.jl to be loaded.

TestParticle.get_curvature_radius Method

julia

get_curvature_radius(x, t, Bfunc)

Calculate the radius of curvature of the magnetic field at position x and time t. Returns Inf if the field is zero or the field lines are straight.

TestParticle.get_dv! Method

julia

get_dv!(dv, v, x, p, t)

In-place solver components for DynamicalODEProblem (velocity).

TestParticle.get_dx! Method

julia

get_dx!(dx, v, x, p, t)

In-place solver components for DynamicalODEProblem (location).

TestParticle.get_energy Method

Calculate the energy [eV] of a relativistic particle from γv in [m/s].

TestParticle.get_energy Method

Return the energy [eV] from relativistic sol.

TestParticle.get_fields Method

julia

get_fields(sol::AbstractODESolution)

Return the electric and magnetic fields from the solution sol.

TestParticle.get_first_crossing Function

julia

get_first_crossing(sol, surface) -> SVector{6}
get_first_crossing(sols, surface) -> Vector{SVector{6}}

State [x, y, z, vx, vy, vz] at the first crossing of surface, linearly interpolated between the two bracketing saved points. Returns a state filled with NaN for trajectories that never reach surface, which is how unreachable velocities are flagged when tracing a detector grid backward to a source. Requires Meshes.jl to be loaded.

TestParticle.get_gc Method

julia

get_gc(xu, param; use_vE = false)
get_gc(x, y, z, vx, vy, vz, bx, by, bz, q2m)

Calculate the coordinates of the guiding center according to the phase space coordinates of a particle. Reference: wiki

Nonrelativistic definition:

X = x - m \frac{b \times v}{q B}

If use_vE=true, it calculates the guiding center relative to the E × B drift frame:

X = x - m \frac{b \times (v - v_{E})}{q B}

which is essential for correctly identifying the polarization drift in time-varying fields.

TestParticle.get_gc_func Method

julia

get_gc_func(param)

Return the function for plotting the orbit of guiding center.

Example

julia

param = prepare(E, B; species = Proton)
# The definitions of stateinit, tspan, E and B are ignored.
prob = ODEProblem(trace!, stateinit, tspan, param)
sol = solve(prob, Vern7(); dt = 2e-11)

f = Figure(fontsize = 18)
ax = Axis3(f[1, 1], aspect = :data)
gc = param |> get_gc_func
gc_plot(x, y, z, vx, vy, vz) = (gc(SA[x, y, z, vx, vy, vz])...,)
lines!(ax, sol, idxs = (gc_plot, 1, 2, 3, 4, 5, 6))

TestParticle.get_gc_velocity Method

julia

get_gc_velocity(y, p, t)

Get the guiding center velocity.

TestParticle.get_gyrofrequency Function

julia

get_gyrofrequency(B=5e-9; q=qᵢ, m=mᵢ)

Return the gyrofrequency [rad/s].

Arguments

B: Magnetic field magnitude [T]. Default is 5 nT.
q: Charge [C]. Default is proton charge.
m: Mass [kg]. Default is proton mass.

TestParticle.get_gyroperiod Function

julia

get_gyroperiod(B=5e-9; q=qᵢ, m=mᵢ)

Return the gyroperiod [s].

Arguments

B: Magnetic field magnitude [T]. Default is 5 nT.
q: Charge [C]. Default is proton charge.
m: Mass [kg]. Default is proton mass.

TestParticle.get_gyroradius Method

julia

get_gyroradius(V, B; q=qᵢ, m=mᵢ)

Return the gyroradius [m].

Arguments

V: Velocity magnitude [m/s] (usually perpendicular to the magnetic field).
B: Magnetic field magnitude [T].
q: Charge [C]. Default is proton charge.
m: Mass [kg]. Default is proton mass.

TestParticle.get_gyroradius Method

julia

get_gyroradius(sol::AbstractODESolution, t)

Return the gyroradius [m] from the solution sol at time t. The interpolated magnetic field function is obtained from sol.prob.p.

TestParticle.get_magnetic_properties Method

julia

get_magnetic_properties(x, t, Bfunc)

Calculate magnetic field properties at position x and time t. Returns tuple (B, ∇B, κ, b̂, Bmag):

B: Magnetic field vector
∇B: Gradient of magnetic field magnitude
κ: Curvature vector
b̂: Unit magnetic field vector
Bmag: Magnitude of B

TestParticle.get_mean_magnitude Method

julia

get_mean_magnitude(B)

Calculate the Root Mean Square (RMS) magnitude of a vector field B. It is assumed that the first dimension of B represents the vector components. This function is compatible with any spatial dimension.

TestParticle.get_particle_crossings Function

julia

get_particle_crossings(sol, surface, weight = 1.0) -> (velocities, weights)
get_particle_crossings(sols, surface, weights = 1.0) -> (velocities, weights)
get_particle_crossings(sols, surfaces, weights = 1.0) -> (velocities, weights)

Detect every crossing of a virtual detector surface and return the velocity [m/s] at each crossing together with the corresponding weight.

A trajectory can cross the same detector more than once, so a crossing carries the weight of its trajectory rather than of the trajectory set: passing a scalar weight applies it to every trajectory of the ensemble. With a vector of surfaces the two returned vectors hold one vector per surface. Requires Meshes.jl to be loaded.

TestParticle.get_particle_flux Function

julia

get_particle_flux(sol, surface, weight = 1.0) -> (number_flux, velocity_flux)

Number flux density and velocity flux density through surface, obtained by summing the crossing weights of a single trajectory and dividing by the surface area. Requires Meshes.jl to be loaded.

TestParticle.get_particle_fluxes Function

julia

get_particle_fluxes(sols, surface, weights = 1.0) -> (number_flux, velocity_flux)
get_particle_fluxes(sols, surfaces, weights = 1.0) -> (number_fluxes, velocity_fluxes)

Number flux densities and velocity flux densities through one surface or through each of surfaces, summed over an ensemble of trajectories and divided by the surface area. Requires Meshes.jl to be loaded.

TestParticle.get_velocity Method

Return velocity from relativistic γv in sol.

TestParticle.get_work Method

julia

get_work(sol::AbstractODESolution)

Return the work done by the electric field from the solution sol.

TestParticle.get_work_rates_gc Method

julia

get_work_rates_gc(xv, p, t)

Calculate the work rates done by the electric field and the betatron acceleration for guiding center.

TestParticle.makegrid Function

julia

makegrid(grid::Union{CartesianGrid, RectilinearGrid})

Return the coordinate ranges of a Meshes.jl grid, one range per dimension, e.g. (x, y, z) for a 3-D grid. Requires Meshes.jl to be loaded.

TestParticle.prepare Function

julia

prepare(args...; kwargs...) -> (q2m, m, E, B, F)
prepare(E, B, F = ZeroField(); kwargs...)
prepare(grid::CartesianGrid, E, B, F = ZeroField(); kwargs...)
prepare(x, E, B, F = ZeroField(); dir = 1, kwargs...)
prepare(x, y, E, B, F = ZeroField(); kwargs...)
prepare(x, y, z, E, B, F = ZeroField(); kwargs...)
prepare(B; E = ZeroField(), F = ZeroField(), kwargs...)

Return a tuple consists of particle charge-mass ratio for a prescribed species of charge q and mass m, mass m for a prescribed species, analytic/interpolated EM field functions, and external force F.

Prescribed species are Electron and Proton; other species can be manually specified with m and q keywords or species = Ion(m̄, q̄), where m̄ and q̄ are the mass and charge numbers respectively.

Direct range input for uniform grid in 1/2/3D is supported. The grid vectors must be sorted. For 1D grid, an additional keyword dir is used for specifying the spatial direction, 1 -> x, 2 -> y, 3 -> z. For 3D grid, the default grid type is CartesianGrid. To use StructuredGrid (spherical) grid, an additional keyword gridtype is needed. For StructuredGrid (spherical) grid, dimensions of field arrays should be (Br, Bθ, Bϕ).

Keywords

order::Int=1: order of interpolation in [0,1,3].
bc=FillExtrap(NaN): boundary condition type from FastInterpolations.jl.
store=StorePolicy(): interpolation storage policy. Use StorePolicy(copy=false) for zero-copy construction when field and grid inputs will remain unchanged for the interpolator's lifetime.
species=Proton: particle species.
q=nothing: particle charge.
m=nothing: particle mass.
gridtype: CartesianGrid, RectilinearGrid, StructuredGrid.

TestParticle.prepare Method

julia

prepare(pf::PotentialField; species = Proton, q = nothing, m = nothing, c = c, type = nothing)

Return parameters (q, m, c, pf) for TraceCanonicalProblem.

TestParticle.prepare_canonical Method

julia

prepare_canonical(phi, A; kwargs...)
prepare_canonical(A; phi = ZeroField(), kwargs...)

Convenience wrapper to construct a PotentialField and prepare parameters for TraceCanonicalProblem.

TestParticle.prepare_gc Method

julia

prepare_gc(xv, xrange, yrange, zrange, E, B;
    species = Proton, q = nothing, m = nothing, order::Int = 1,
    bc = FillExtrap(NaN), store = StorePolicy())
prepare_gc(xv, E, B; species = Proton, q = nothing, m = nothing)

Prepare the guiding center parameters for a particle.

TODO: support external forces

TestParticle.project_vdf Method

julia

project_vdf(f3d, dv) -> (f_xy, f_xz, f_yz)

Integrate a 3-D phase-space density f3d in [s³/km⁶] over each velocity axis in turn, returning the projections f(vᵢ, vⱼ) = ∫ f dvₖ in [s²/km⁵]. dv is the bin width in [km/s]. The array axes are ordered (v₁, v₂, v₃), so the projections come out as (v₁, v₂), (v₁, v₃) and (v₂, v₃).

TestParticle.refine_vdf_window Method

julia

refine_vdf_window(f, vx, vy, vz, v0, dv; margin = 0.0, relthresh = 1.0e-5,
    bounds = nothing)

Refine a coarse velocity grid down to the spacing dv, keeping only the populated part of f. This is the second pass of the adaptive backward Liouville method: a cheap coarse pass locates the region of velocity space that the detector can actually see, and only that region is retraced on the fine grid.

Returns (vx_new, vy_new, vz_new) in the same unit as the inputs: the smallest box enclosing every grid point where f exceeds relthresh times its maximum, expanded by margin and clamped to bounds. The box bounds are snapped onto v0 + k·dv so that the refined grid shares bin centers with the full grid starting at v0. v0, dv and margin may be scalars or 3-tuples for per-axis values.

bounds is a (lo, hi) pair, or one such pair per axis. It defaults to the extent of the coarse grid, which is the range the coarse pass actually probed; pass the extent of the full display grid instead when the refined result has to live on that grid.

TestParticle.relative_l2 Method

julia

relative_l2(rec, ref; thresh = 1.0e-4)

Relative L2 difference ‖rec − ref‖ / ‖ref‖, evaluated over the entries where ref is above thresh times its own maximum and rec is finite. Restricting to the populated cells keeps the metric from being dominated by the empty tails of a reconstruction.

TestParticle.sample_maxwellian Function

julia

sample_maxwellian([rng], Tn, m; offset=0.0, u0=SA[0.0, 0.0, 0.0])
sample_maxwellian([rng], Tn, species::String; offset=0.0, u0=SA[0.0, 0.0, 0.0])

Sample a velocity [m/s] from a Maxwellian distribution with temperature Tn [K], mass m [kg], energy offset offset [J], and bulk velocity u0 [m/s]. This function requires VelocityDistributionFunctions.jl to be loaded.

TestParticle.sample_unit_sphere Function

julia

sample_unit_sphere(rng::AbstractRNG=default_rng())

Sample a unit vector on a sphere uniformly.

TestParticle.sample_velocity_ball Method

julia

sample_velocity_ball([rng], radius; center = SA[0.0, 0.0, 0.0]) -> SVector{3}

Sample a velocity [m/s] uniformly from the ball of radius [m/s] centered on center [m/s].

Uniformity is what makes the sample usable as a quadrature rule in velocity space: each draw stands for the same volume 4πR³ / (3N), which forward Liouville tracking needs in order to convert a sampled f into a contribution to a detector bin.

TestParticle.sph2cart Method

julia

sph2cart(r, θ, ϕ)

Convert from spherical to Cartesian coordinates vector.

TestParticle.sph2cartvec Method

julia

sph2cartvec(vr, vθ, vϕ, θ, ϕ)

Convert a vector from spherical to Cartesian coordinates. θ and ϕ define the position (polar and azimuthal angles in radians) where the spherical components (vr, vθ, vϕ) are defined.

TestParticle.trace! Method

julia

trace!(dy, y, p, t)

ODE equations for charged particle moving in EM field and external force field with in-place form.

TestParticle.trace Method

julia

trace(y, p, t)::SVector{6}

ODE equations for charged particle moving in EM field and external force field with out-of-place form.

TestParticle.trace_fieldline! Method

julia

trace_fieldline!(dx, x, p, s)

Equation for tracing magnetic field lines with in-place form. The parameter p is the magnetic field function. Note that the independent variable s represents the arc length.

TestParticle.trace_fieldline Method

julia

trace_fieldline(x, p, s)

Equation for tracing magnetic field lines with out-of-place form.

TestParticle.trace_gc! Method

julia

trace_gc!(dy, y, p, t)

Guiding center equations for nonrelativistic charged particle moving in EM field with in-place form. Variable y = (x, y, z, u), where u is the velocity along the magnetic field at (x,y,z).

TestParticle.trace_gc_drifts! Method

julia

trace_gc_drifts!(dx, x, p, t)

Equations for tracing the guiding center using analytical drifts, including the grad-B drift, curvature drift, and ExB drift. Parallel velocity is also added. This expression requires the full particle trajectory p.sol.

TestParticle.trace_gc_exb! Method

julia

trace_gc_exb!(dx, x, p, t)

Equations for tracing the guiding center using the ExB drift and parallel velocity from a reference trajectory.

TestParticle.trace_gc_flr! Method

julia

trace_gc_flr!(dx, x, p, t)

Equations for tracing the guiding center using the ExB drift with FLR corrections and parallel velocity.

TestParticle.trace_normalized! Method

julia

trace_normalized!(dy, y, p, t)

Normalized ODE equations for charged particle moving in EM field with in-place form. If the field is in 2D X-Y plane, periodic boundary should be applied for the field in z via the extrapolation function.

TestParticle.trace_normalized Method

julia

trace_normalized(y, p, t)

Normalized ODE equations for charged particle moving in EM field with out-of-place form.

TestParticle.trace_relativistic! Method

julia

trace_relativistic!(dy, y, p, t)

ODE equations for relativistic charged particle (x, γv) moving in EM field with in-place form.

TestParticle.trace_relativistic Method

julia

trace_relativistic(y, p, t) -> SVector{6}

ODE equations for relativistic charged particle (x, γv) moving in static EM field with out-of-place form.

TestParticle.trace_relativistic_normalized! Method

julia

trace_relativistic_normalized!(dy, y, p, t)

Normalized ODE equations for relativistic charged particle (x, γv) moving in EM field with in-place form.

TestParticle.trace_relativistic_normalized Method

julia

trace_relativistic_normalized(y, p, t)

Normalized ODE equations for relativistic charged particle (x, γv) moving in EM field with out-of-place form.

TestParticle.vdf_backward Method

julia

vdf_backward(sols, source, f_src, dims) -> Array{T, 3}

Phase-space density [s³/km⁶] on a velocity grid of size dims = (nx, ny, nz), obtained by tracing that grid backward in time and evaluating the source distribution on the traced-back states, as prescribed by Liouville's theorem.

Each trajectory of sols is followed back to the surface source with get_first_crossing and f_src is evaluated on the velocity there; grid points whose trajectory never reaches source are left at zero. f_src takes a velocity in [m/s] and returns f in [s³/m⁶].

TestParticle.vdf_backward_trace Method

julia

vdf_backward_trace(
    param, detector, source, f_src;
    v_range = 1000.0e3, vy_range = v_range, vz_range = v_range,
    dv = 20.0e3, dt, tspan = (0.0, -20.0),
    adaptive = true, dv_coarse = 3 * dv, margin = 3 * dv,
    relthresh = 1.0e-5, bounds = nothing, isoutside = nothing,
    alg = Boris(), ensemblealg = EnsembleThreads()
) -> (f_3d, (vx_grid, vy_grid, vz_grid))

Reconstruct the 3-D phase-space density [s³/km⁶] at detector by tracing a velocity grid backward in time to source and evaluating f_src on the traced-back states.

TestParticle.vdf_forward_trace Method

julia

vdf_forward_trace(
    param, source, detector, f_src;
    nparticles = 50000, vradius, tspan = (0.0, 20.0), dt,
    center = SA[0.0, 0.0, 0.0], dv_km = 20.0,
    v_edges = -1000.0:dv_km:1000.0,
    isoutside = nothing, alg = Boris(), ensemblealg = EnsembleThreads(),
    seed = 42
) -> (f_3d, (f_xy, f_xz, f_yz))

Reconstruct the 3-D phase-space density [s³/km⁶] at detector using forward Liouville tracking from a uniformly sampled velocity ball at source.

TestParticle.vdf_grid_problem Method

julia

vdf_grid_problem(vx, vy, vz, x0, param, tspan) -> TraceProblem

Build a TraceProblem that traces a velocity-space grid: trajectory i starts at position x0 [m] with the grid velocity (vx[ix], vy[iy], vz[iz]) [m/s], the z index varying fastest. This is the setup of the backward Liouville method, where the grid sits at the detector and tspan runs backward in time.

TestParticle.velocity_moments Method

julia

velocity_moments(v, A; dv)
velocity_moments(proj::Tuple; dv)

Density n [m⁻³] and momentum density n·V [m⁻³·km/s] of a 2-D velocity distribution A in [s²/km⁵], sampled on the first-axis bin centers v in [km/s] with bin width dv in [km/s]. Since the momentum is taken along the first axis, pass the Vx–Vy or the Vx–Vz projection to obtain n·Vx.

Private types and methods

TestParticle.AbstractFieldInterpolator Type

julia

AbstractFieldInterpolator

Abstract type for all field interpolators.

TestParticle.AbstractTraceSolution Type

Abstract type for tracing solutions.

TestParticle.FieldInterpolator Type

julia

FieldInterpolator{T}

A callable struct that wraps a 3D interpolation object.

TestParticle.FieldInterpolator1D Type

julia

FieldInterpolator1D{T}

A callable struct that wraps a 1D interpolation object.

TestParticle.FieldInterpolator2D Type

julia

FieldInterpolator2D{T}

A callable struct that wraps a 2D interpolation object.

TestParticle.SphericalFieldInterpolator Type

julia

SphericalFieldInterpolator{T}

A callable struct for spherical grid interpolation (scalar or combined vector).

TestParticle.VDF Type

Abstract type for velocity distribution functions.

TestParticle._eval_f_src Method

julia

vdf_forward_trace(
    param, source, detector, vdf, n0;
    nparticles = 50000, vradius, tspan = (0.0, 20.0), dt,
    center = SA[0.0, 0.0, 0.0], dv_km = 20.0,
    v_edges = -1000.0:dv_km:1000.0,
    isoutside = nothing, alg = Boris(), ensemblealg = EnsembleThreads(),
    seed = 42
) -> (f_3d, (f_xy, f_xz, f_yz))

Reconstruct the 3-D phase-space density [s³/km⁶] at detector using forward Liouville tracking from a uniformly sampled velocity ball at source.

TestParticle._gc_direction Method

julia

_gc_direction(prob) -> Int

+1 for a forward time span and -1 for a backward one.

TestParticle._gc_time_type Method

julia

_gc_time_type(prob, dt) -> Type

The type of the time values a solve will produce.

TestParticle._gc_to_full Function

julia

_gc_to_full(state_gc, E_field, B_field, q, m, μ, t, phase = 0)

Internal GC → full-orbit conversion taking the field functions and charges directly, so it can be shared by the GC and hybrid solvers. E_field and B_field must be callable as E_field(x, t) / B_field(x, t); static fields simply ignore t.

TestParticle._get_gc_parameters Method

julia

_get_gc_parameters(xv, E, B, q, m, t) -> (X, vpar, q, q / m, μ, phase)

Extract the guiding-center position X, parallel velocity vpar, magnetic moment μ, and gyrophase phase from a full-particle state xv at time t. E and B must be callable as E(x, t) / B(x, t); static fields simply ignore t.

#TODO : support external forces

TestParticle._hybrid_save! Method

julia

_hybrid_save!(traj, tsave, plan, isave, t_last, y_last, t_out, y_out)

Record an output sample (t_out, y_out). With saveat in force the sample is kept only as the right-hand side of the interpolation, and every requested time it has passed is reported in between, so the integration itself is untouched. Otherwise the sample is appended as it is.

Returns the updated cursor and previous sample.

TestParticle._rk4! Method

Apply RK4 method for particles with index in irange.

TestParticle._rk45! Method

Apply RK45 method for particles with index in irange.

TestParticle._rng_seed Method

julia

_rng_seed(seed::Integer, i::Integer)

Combine the ensemble seed and the trajectory index i into a well-mixed UInt64 used to seed a per-trajectory RNG.

A naive seed + i couples the two linearly: incrementing seed merely rotates which RNG each trajectory index draws from, so two nearby seeds (e.g. 1 and 2) produce ensembles that share almost all trajectories and are therefore not independent. The SplitMix64-style finalizer below is a bijection on UInt64, so any distinct seed decorrelates every trajectory while remaining deterministic and reproducible (no reliance on hash, which is salted per Julia session).

Note: Julia's multi-argument Xoshiro(seed, i, …) constructor does not mix — it embeds the integers directly into the 4×UInt64 state (Xoshiro(s0,s1,s2,s3) → initstate!). Because xoshiro256++'s first output depends only on s0 and s3, seeding with Xoshiro(seed, i, 0, 0) yields the same first draw for every trajectory index i and does not decorrelate adjacent seeds. Robust mixing in Xoshiro only happens via the single-argument constructor (SHA2-256 hash_seed), so seed and i must be folded into one well-mixed integer first — which is what this function does.

TestParticle.get_curvature Method

julia

get_curvature(x, t, Bfunc)

Calculate the curvature vector κ of the magnetic field at position x and time t.

TestParticle.get_perp_vector Method

julia

get_perp_vector(b::AbstractVector)

Obtain two unit vectors e1 and e2 such that (e1, e2, b) form a right-handed orthonormal system.

TestParticle.get_work_rates Function

julia

get_work_rates(xu, p, t)

Calculate the work rates done by the electric field and the betatron acceleration. Returns a tuple (P_par, P_fermi, P_grad, P_betatron).

TestParticle.is_time_dependent Method

julia

is_time_dependent(f::Function)

Judge whether the field function is time dependent.

TestParticle.trace_gc Method

julia

trace_gc(y, p, t)

Guiding center equations for nonrelativistic charged particle moving in EM field with out-of-place form.

TestParticle.update_dp5 Method

julia

update_dp5(y, param, dt, t)::(dy, E)

Update state using the Dormand-Prince 5(4) method. Returns dy as the update and E as the error estimate, both as SArrays.

TestParticle.update_rk4 Method

julia

update_rk4(y, param, dt, t)

Update state using the RK4 method. Returns dy as the update SArray.