Internal
Public APIs
TestParticle.AdaptiveBoris Method
AdaptiveBoris(; dtmin, dtmax, safety=0.1)Adaptive Boris method with adaptive time stepping based on local gyrofrequency. The time step is determined by dt = safety / |qB/m|, clamped by dtmin and dtmax.
TestParticle.AdaptiveHybrid Type
AdaptiveHybrid{T}Hybrid algorithm that switches between adaptive Boris (FO) and RK45 (GC).
sourceTestParticle.LazyTimeInterpolator Type
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: Functioni -> fieldthat loads the field at indexi.buffer::Dict{Int, F}: Cache for loaded fields.lock::ReentrantLock: Lock for thread safety.
TestParticle.TraceGCProblem Type
TraceGCProblem{uType, tType, isinplace, P, F <: AbstractODEFunction, PF} <: AbstractODEProblem{uType, tType, isinplace}Problem type for tracing guiding centers.
sourceTestParticle.TraceHybridProblem Type
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.BiKappa Function
BiKappa(args...; kw...)Construct a BiKappa distribution. Forwards to VelocityDistributionFunctions.BiKappa.
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
BiMaxwellian(args...; kw...)Construct a BiMaxwellian distribution. Forwards to VelocityDistributionFunctions.BiMaxwellian.
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
Kappa(args...; kw...)Construct a Kappa distribution. Forwards to VelocityDistributionFunctions.Kappa.
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
Maxwellian(args...; kw...)Construct a Maxwellian distribution. Forwards to VelocityDistributionFunctions.Maxwellian.
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.full_to_gc Method
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
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.get_adiabaticity Function
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
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
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_curvature_radius Method
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_energy Method
Calculate the energy [eV] of a relativistic particle from γv in [m/s].
sourceTestParticle.get_fields Method
get_fields(sol::AbstractODESolution)Return the electric and magnetic fields from the solution sol.
TestParticle.get_gc Method
get_gc(xu, param)
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:
TestParticle.get_gc_func Method
get_gc_func(param)Return the function for plotting the orbit of guiding center.
Example
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
get_gc_velocity(y, p, t)Get the guiding center velocity.
sourceTestParticle.get_gyrofrequency Function
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
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
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
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_mean_magnitude Method
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_number_density_flux Method
get_number_density_flux(grid::CartesianGrid, sols, dt)Calculate the steady state particle number density flux on a uniform Cartesian grid. The flux is estimated by accumulating the number of particles in each cell at time steps dt, divided by the surface area of each cell.
Arguments
grid: ACartesianGridfromMeshes.jl.sols: Particle trajectory solutions (e.g.EnsembleSolution).dt: Time step for sampling particle positions.
TestParticle.get_work Method
get_work(sol::AbstractODESolution)Return the work done by the electric field from the solution sol.
TestParticle.prepare Function
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. 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 [1,2,3].bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic.species=Proton: particle species.q=nothing: particle charge.m=nothing: particle mass.gridtype:CartesianGrid,RectilinearGrid,StructuredGrid.
TestParticle.prepare_gc Method
prepare_gc(xv, xrange, yrange, zrange, E, B;
species = Proton, q = nothing, m = nothing, order::Int = 1, bc::Int = 1)
prepare_gc(xv, E, B; species = Proton, q = nothing, m = nothing)Prepare the guiding center parameters for a particle.
sourceTestParticle.sample_unit_sphere Method
sample_unit_sphere()Sample a unit vector on a sphere uniformly.
sourceTestParticle.trace! Method
trace!(dy, y, p, t)ODE equations for charged particle moving in EM field and external force field with in-place form.
sourceTestParticle.trace Method
trace(y, p, t)::SVector{6}ODE equations for charged particle moving in EM field and external force field with out-of-place form.
sourceTestParticle.trace_fieldline! Method
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
trace_fieldline(x, p, s)Equation for tracing magnetic field lines with out-of-place form.
sourceTestParticle.trace_fieldline Method
trace_fieldline(u0, p, tspan::Tuple{<:Real, <:Real}; mode::Symbol=:both, kw...)Helper function to create ODEProblems for tracing magnetic field lines.
Arguments
u0: Initial position.p: Parameter tuple containing the magnetic field. It can also be the magnetic field array (numerical) or function (analytical/interpolator).tspan: Time span (arc length) for tracing. typically(0.0, L).mode: Tracing mode, one of:forward,:backward,:both.kw: Keyword arguments passed topreparewhenpis an array.
Returns
If
modeis:forwardor:backward, returns a singleODEProblem.If
modeis:both, returns a vector of twoODEProblems (forward and backward).
TestParticle.trace_gc! Method
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
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
trace_gc_exb!(dx, x, p, t)Equations for tracing the guiding center using the ExB drift and parallel velocity from a reference trajectory.
sourceTestParticle.trace_gc_flr! Method
trace_gc_flr!(dx, x, p, t)Equations for tracing the guiding center using the ExB drift with FLR corrections and parallel velocity.
sourceTestParticle.trace_normalized! Method
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 provided by Interpolations.jl.
sourceTestParticle.trace_normalized Method
trace_normalized(y, p, t)Normalized ODE equations for charged particle moving in EM field with out-of-place form.
sourceTestParticle.trace_relativistic! Method
trace_relativistic!(dy, y, p, t)ODE equations for relativistic charged particle (x, γv) moving in EM field with in-place form.
sourceTestParticle.trace_relativistic Method
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.
sourceTestParticle.trace_relativistic_normalized! Method
trace_relativistic_normalized!(dy, y, p, t)Normalized ODE equations for relativistic charged particle (x, γv) moving in EM field with in-place form.
sourceTestParticle.trace_relativistic_normalized Method
trace_relativistic_normalized(y, p, t)Normalized ODE equations for relativistic charged particle (x, γv) moving in EM field with out-of-place form.
sourcePrivate types and methods
TestParticle.AbstractFieldInterpolator Type
AbstractFieldInterpolatorAbstract type for all field interpolators.
sourceTestParticle.Field Type
Field{itd, F} <: AbstractField{itd}A representation of a field function f, defined by:
time-independent field
time-dependent field
Arguments
field_function::Function: the function of field.itd::Bool: whether the field function is time dependent.F: the type offield_function.
TestParticle.FieldInterpolator Type
FieldInterpolator{T}A callable struct that wraps a 3D interpolation object.
sourceTestParticle.FieldInterpolator1D Type
FieldInterpolator1D{T}A callable struct that wraps a 1D interpolation object.
sourceTestParticle.FieldInterpolator2D Type
FieldInterpolator2D{T}A callable struct that wraps a 2D interpolation object.
sourceTestParticle.SphericalFieldInterpolator Type
SphericalFieldInterpolator{T}A callable struct for spherical grid interpolation (scalar or combined vector).
sourceTestParticle._solve_single_adaptive_boris Method
See _solve_single_boris for the rationale behind the single_prob construction.
TestParticle._solve_single_boris Method
_solve_single_boris(prob, i, ...)Solve a single trajectory i of prob for use with EnsembleDistributed.
_generic_boris! uses the loop index i for two purposes simultaneously: applying prob_func(prob, i, false) to select per-particle initial conditions, and storing the result at sols[i]. For a distributed worker handling only one trajectory at a time, a 1-element local_sols would be out of bounds if i > 1 were passed directly. Decoupling the two uses would require threading a storage offset through the entire _dispatch_boris! → _generic_boris! call chain.
Instead, we pre-apply prob_func to get the correct IC for trajectory i and wrap the result in a fresh TraceProblem with the default (identity) prob_func, so _generic_boris! can safely iterate 1:1 without applying prob_func again. The TraceProblem construction is a negligible struct copy relative to the simulation cost and the serialization overhead inherent in pmap.
TestParticle.adapt_field_to_gpu Method
adapt_field_to_gpu(field::Field, backend::KA.Backend)Adapt interpolation fields to GPU memory using Adapt.jl. Analytic functions are returned unchanged.
sourceTestParticle.boris_velocity_update Method
boris_velocity_update(v, E, B, qdt_2m)Update velocity using the Boris method, returning the new velocity as an SVector. This is the core logic shared between the standard solver and the kernel solver.
sourceTestParticle.cart2sph Method
cart2sph(x, y, z)Convert from Cartesian to spherical coordinates vector.
sourceTestParticle.get_cell_centers Method
Return cell center coordinates from 2D/3D RectilinearGrid.
sourceTestParticle.get_curvature Method
get_curvature(x, t, Bfunc)Calculate the curvature vector κ of the magnetic field at position x and time t.
TestParticle.get_interpolator Method
get_interpolator(A, gridx, gridy, gridz, order::Int=1, bc::Int=1)
get_interpolator(gridtype, A, grid1, grid2, grid3, order::Int=1, bc::Int=1)Return a function for interpolating field array A on the grid given by gridx, gridy, and gridz.
Arguments
gridtype:CartesianGrid,RectilinearGridorStructuredGrid.A: field array. For vector field, the first dimension should be 3.order::Int=1: order of interpolation in [1,2,3].bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic, 3 -> Flat.
Notes
The input array A may be modified in-place for memory optimization.
TestParticle.get_magnetic_properties Method
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 vectorb̂: Unit magnetic field vectorBmag: Magnitude of B
TestParticle.get_number_density Method
get_number_density(sols, grid, t_start, t_end, dt)Calculate time-averaged particle number density from t_start to t_end with step dt.
TestParticle.get_number_density Method
get_number_density(sols, grid, t)Calculate particle number density at time t in a given grid.
TestParticle.get_perp_vector Method
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_rotation_matrix Method
get_rotation_matrix(axis::AbstractVector, angle) :: SMatrix{3,3}Create a rotation matrix for rotating a 3D vector around a unit axis by an angle in radians. Reference: Rotation matrix from axis and angle
Example
using LinearAlgebra
v = [-0.5, 1.0, 1.0]
v̂ = normalize(v)
θ = deg2rad(-74)
R = get_rotation_matrix(v̂, θ)TestParticle.get_work_rates Function
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.get_work_rates_gc Method
get_work_rates_gc(xv, p, t)Calculate the work rates done by the electric field and the betatron acceleration for guiding center.
sourceTestParticle.getinterp Function
getinterp(::Type{<:CartesianGrid}, A, gridx, gridy, gridz, order::Int=1, bc::Int=1)Return a function for interpolating field array A on the grid given by gridx, gridy, and gridz.
Arguments
order::Int=1: order of interpolation in [1,2,3].bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic, 3 -> Flat.dir::Int: 1/2/3, representing x/y/z direction.
Notes
The input array A may be modified in-place for memory optimization.
TestParticle.getinterp_scalar Function
getinterp_scalar(::Type{<:CartesianGrid}, A, gridx, gridy, gridz, order::Int=1, bc::Int=1)Return a function for interpolating scalar array A on the grid given by gridx, gridy, and gridz. Currently only 3D arrays are supported.
Arguments
order::Int=1: order of interpolation in [1,2,3].bc::Int=1: type of boundary conditions, 1 -> NaN, 2 -> periodic, 3 -> Flat.dir::Int: 1/2/3, representing x/y/z direction.
TestParticle.solve Function
solve(prob::TraceGCProblem; trajectories::Int=1, dt::AbstractFloat,
savestepinterval::Int=1, isoutofdomain::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.
TestParticle.solve Method
solve(prob::TraceProblem; trajectories::Int=1, dt::AbstractFloat,
savestepinterval::Int=1, isoutofdomain::Function=ODE_DEFAULT_ISOUTOFDOMAIN,
n::Int=1, save_start::Bool=true, save_end::Bool=true, save_everystep::Bool=true,
save_fields::Bool=false, save_work::Bool=false)Trace particles using the Boris method with specified prob.
keywords
trajectories::Int: number of trajectories to trace.dt::AbstractFloat: time step.savestepinterval::Int: saving output interval.isoutofdomain::Function: a function with input of position and velocity vectorxvthat determines whether to stop tracing.n::Int: number of substeps for the Multistep Boris method. Default is 1 (standard Boris).save_start::Bool: save the initial condition. Default istrue.save_end::Bool: save the final condition. Default istrue.save_everystep::Bool: save the state at everysavestepinterval. Default istrue.save_fields::Bool: save the electric and magnetic fields. Default isfalse.save_work::Bool: save the work done by the electric field. Default isfalse.
TestParticle.sph2cart Method
sph2cart(r, θ, ϕ)Convert from spherical to Cartesian coordinates vector.
sourceTestParticle.sph_to_cart_vector Method
sph_to_cart_vector(vr, vθ, vϕ, θ, ϕ)Convert a vector from spherical to Cartesian coordinates. θ and ϕ are defined in radians.
sourceTestParticle.trace_gc Method
trace_gc(y, p, t)Guiding center equations for nonrelativistic charged particle moving in EM field with out-of-place form.
sourceTestParticle.update_dp5 Method
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.
sourceTestParticle.update_rk4 Method
update_rk4(y, param, dt, t)Update state using the RK4 method. Returns dy as the update SArray.
sourceTestParticle.update_velocity! Method
update_velocity!(xv, paramBoris, param, dt, t)Update velocity using the Boris method, Birdsall, Plasma Physics via Computer Simulation. Reference: DTIC
sourceTestParticle.update_velocity Method
update_velocity(v, r, param, dt, t)Update velocity using the Boris method, returning the new velocity as an SVector.
sourceTestParticle.update_velocity_multistep Method
update_velocity_multistep(v, r, param, dt, t, n)Update velocity using the Multistep Boris method, returning the new velocity as an SVector. Reference: Zenitani & Kato 2025
source