rerun_nonthermals(fname, nparticles)

Rerun the nparticles most energetic particles from the experiment defined by fname.

rerun_maxiters(fname, nparticles)

Rerun the nparticles particles with most timesteps from the experiment defined by fname.

rerun_idxs(fname, idxs)

Rerun the the test particles of index idxs from the experiment defined by fname.

guidingcentreapproximation!(du, u, p, _)

The equations of motion of a charged particle in collisionless plasma under the guiding centre approximation (GCA). Compared to the full motion (described by the Lorentz force), the statevector u is reduced from 6 DoF to 4, namely the guiding centre position and the guiding centre velocity parallel to the magnetic field, [Rx, Ry, Rz, vparal], respectively.

The parameters p are the particle charge, mass and magnetic moment (which is assumed to be constant in the GCA), along with the interpolation functor from Interpolations.jl giving the magnetic and electric field at an arbitrary position. The functor should give a 6-component vector, where the first 3 components represents the magnetic field and the last 3 the electric field.

The function uses ForwardDiff and Interpolations to calculate gradients.

gca_2Dxz!(du, u, p, _)

The equations of motion of a charged particle in collisionless plasma under the guiding centre approximation (GCA) in 2.5D xz-plane. Compared to the full motion (described by the Lorentz force), the statevector u is reduced from 6 DoF to 4, namely the guiding centre position and the guiding centre velocity parallel to the magnetic field, [Rx, Ry, Rz, vparal], respectively. In practice, Ry is necessary.

The parameters p are the particle charge, mass and magnetic moment (which is assumed to be constant in the GCA), along with the interpolation functor from Interpolations.jl giving the magnetic and electric field at an arbitrary position. The functor should give a 6-component vector, where the first 3 components represents the magnetic field and the last 3 the electric field.

The function uses ForwardDiff and Interpolations to calculate gradients.

lorentzforce!(du, u, p, _)

Equations of motion described by the Lorentz Force in 3 dimensions. The parameters p contains the charge and mass of the particle, and an interpolation functor giving the magnetic and electric field by passing the position coordinates.

maxwellianvelocitysample(
    rng::AbstractRNG,
    temperature_itp::Any,
    mass::Real,
    args...
    ;
    precision::DataType=Float64,
)

Draw a velocity component sample from a Maxwellian distribution with temperature given by a function call with aribtrary arguments, and by a mass.

OutOfDomainCondition_2Dxz
    xbounds::Tuple{<:Real, <:Real}
    zbounds::Tuple{<:Real, <:Real}

Returns true if the particle is outside the 2D bounds in x and z.

RelativisticConditionGCA_2Dxz

Callback condition for checking if the GCA particle is relativistic. Returns true if the particle's kinetic energy is a user defined fraction of the rest energy. Perpendicular drifts other than E cross B drift are neglected.

Default fraction is 0.002, which is 1022 eV for an electron.

GCABreakDownCondition_2Dxz

Callback condition for checking GCA assumption. Returns true if the ratio between the particle Larmor radius and the characteristic length of the magnetic field is higher than a tolerance switchtol.

outofdomainaffect!(integrator)

Callback affect which terminates the integration and sets the return message to "OutofDomain".

relativisticaffect!(integrator)

Callback affect which terminates the integration and sets the return message to "Relativistic".

gcabreakdownaffect!(integrator)

Callback affect which terminates the integration and sets the return message to "GCABreakDown".

Draw uniform position, Maxwell-Boltzmann distributed velocity.

Draws x, y, and z positions from a given distribution using rejection sampling.

Draw initial conditions from an array.

output_func_max_lightweight(sol, i)

Output function for ensemble simulations using DifferentialEquations.

Finds the maximum larmor radius, minimum field length and maximum scales ratio of the particle's trajectory using the inherent interpolation function in the particle solution. Also returns the initial and final state of the particle, in addition to its charge, mass and magnetic moment.

We set the required "rerun" return argument to false.

SaveBatchAsCSV(expdir, expname, batchsize, nbatches)

Reduction functor that saves the batch as a DataFrame. The struct instance is callable and its fields are

• expdir::String: The directory where the data is saved.
• expname::String: The name of the experiment.
• batchnr::Int: The current batch number.
• batchsize::Int: The size of the batch.
• nbatches::Int: The total number of batches.

SaveBatchAsHDF5(expdir, expname, batchsize, nbatches)

Reduction functor that saves the batch as an HDF5 file. The struct instance is callable and its fields are

• expdir::String: The directory where the data is saved.
• expname::String: The name of the experiment.
• batchnr::Int: The current batch number.
• batchsize::Int: The size of the batch.
• nbatches::Int: The total number of batches.

get_filename(reduction::SaveBatchAsHDF5)

Returns the filename for the current batch from the reduction functor.

h5_getall(filename; batches=nothing)

Returns all datasats from a set of batches in a HDF5-file as a DataFrame. If batches is not specified, the data from all batches are returned.

h5_getbatch(filename, batchnr)

Returns batch number batchnr as a DataFrame.

h5_getdataset(filename, dataset, batchnr)

Returns dataset dataset from batch number batchnr in a HDF5-file.

h5_getdataset(filename, dataset; batches=nothing)

Returns the dataset dataset multiple batches in a HDF5-file. If batches is not specified, all batches are returned.

h5_getenergies(filename; batches=nothing, units="eV")

Returns the initial and final energies of a set of batches in a HDF5-file. If batches is not specified, the data from all batches are returned.

h5_getinitialstate(filename)

Returns x0, y0, z0, vparal0, and magneticmoment from a test particle ensemble stored in a HDF5-file.

save_gcastates(filename, fields_itp)

Calculates the GCAStates of the initial and final states of an ensamble of test particles.

GCAState(
    state            ::Vector{<:Real},
    charge           ::Real,
    mass             ::Real,
    magneticmoment   ::Real,
    time             ::Real,
    kineticenergy    ::Real,
    perpvelocity     ::Real,
    larmorradius     ::Real,
    gyrofrequency    ::Real,
    cosofpitchangle  ::Real,
    exbdrift         ::Vector{<:Real},
    ∇Bdrift          ::Vector{<:Real},
    curvaturedrift   ::Vector{<:Real},
    polarisationdrift::Vector{<:Real},
    mirroracc        ::Real,
    paralacc         ::Real,
    magneticfield    ::Vector{<:Real},
    electricfield    ::Vector{<:Real},
    magneticfieldlength::Real
    )

A struct to hold the state (guiding centre position and its parallel velocity) of a particle in the guiding centre approximation, all its drifts, the external electromagnetic field, relevant parameters, and auxiliary quantities.

Constructors

GCAState(
    state ::Vector{<:Real},
    q     ::Real,
    m     ::Real,
    μ     ::Real,
    itpvec::Vector{<:AbstractInterpolation}; 
    time=nothing,
    dimensionality="2Dxz"
    )

Constructs a GCAState object from;

• state - The guiding centre state; position and parallel velocity
• vparal - The parallel velocity of the guiding centre
• q - The charge of the particle
• m - The mass of the particle
• μ - The magnetic moment of the particle
• itpvec - A vector of interpolation objects for the electromagnetic fields

Keyword arguments:

• time - The time of the state

• dimensionality - The dimensionality of the problem (default: "2Dxz")

  GCAState( state ::Vector{<:Real}, q ::Real, m ::Real, μ ::Real, bfield::Vector{<:Real}, efield::Vector{<:Real} ; time=nothing )

Constructs a partially filled GCAState object, neglecting all drifts except the ExB-drift.

GCAState(
    solution::DataFrame,
    fields::Vector{<:AbstractInterpolation}
    ;
    kwargs...
    )

From DataFrame test-particle solution and the field interpolation vector used in the simulation, construct two ensembles of GCAStates. One for the initial states of the particles and one for the final states.

Get the observable sym from an ODESolution or a vector of ODESolutions.

Methods

get_observable(vectorofsolutions, symbol; t=nothing, kwargs...)
get_observable(
    solution, symbol;
    times=nothing, 
    tspan=(first(sol.t), last(sol.t)), 
    nsteps=1000, 
    kwargs...

)

get_stateidx(
    sol::ODESolution,
    idx::Int;
    tspan=(first(sol.t), last(sol.t)),
    nsteps=1000,
    times=range(tspan[1], tspan[2], length=nsteps),
)

Return the idx-th component of the state vector of sol at times times.

get_exbdrift(
    sol::ODESolution,
    t::Real,
)

Calculate the magnitude of the ExB-drift at time t of sol.

get_fermi(
    sol::ODESolution,
    t::Real,
)

Calculate the Fermi acceleration at time t of sol, in eV. That is the energy change due to 1st order Fermi acceleration at time t.

get_betatron(
    sol::ODESolution,
    t::Real,
)

Calculate the betatron acceleration at time t of sol, in eV. That is the energy change due to betatron acceleration at time t.

get_polarisationacc(
    sol::ODESolution,
    t::Real,
)

Calculate the polarisation acceleration at time t of sol, in eV. That is the energy change due to polarisation acceleration at time t.

get_driftenergy(
    sol::ODESolution,
    t::Real,
)

Calculate the kinetic energy of the ODESolution sol at time t from perpendicular drifts.

get_perpenergy(
    sol::ODESolution,
    t::Real;
)

Calculate the perpendicular kinetic energy of the ODESolution sol at time t.

get_parallelenergy(
    sol::ODESolution,
    t::Real;
)

Calculate the parallel kinetic energy of the ODESolution sol at time t.

binmap(
    data,
    args...
    ;
    weights=ones(length(data)),
    mapfunc::Function = (x,w) -> sum(w),
    xvalues=nothing,
    xedges=nothing,
    nbins= 100,
    logx=false,
    logf=false,
    normalise=false,
    )

Maps a binning of the data to the function mapfunc. The default mapping is the sum of the weights of the data in each bin. The defualt weighing is 1.

If no edges or data-corresponding x-values are given, we bin the data in a linear, uniform range between its minimum and maximum value.

Returns the midpoints of the bins and the mapped values of the bins.

binmap(
    xvalues,
    yvalues,
    zvalues,
    mapfunc::Function = x -> length(x),
    args...
    ;
    nbinsx=100,
    nbinsy=100,
    logaxes=false,
    logx=false,
    logy=false,
    logz=false,
    normalise=false,
    )

Maps a 2D binning of the zvalues to the function mapfunc. Does not support weighing of the zvalues.

Returns the 2D midpoints of the bins and their mapped values.