Velocity Distribution Functions

This package provides implementations of common velocity distribution functions used in plasma physics.

Maxwellian(args...; u0=nothing, kw...)
Maxwellian(n, args...; u0=nothing, kw...)

BiMaxwellian(args...; u0=nothing, kw...)
BiMaxwellian(n, args...; u0=nothing, kw...)

For Kappa distribution, see Kappa.

Maxwellian Distribution

The isotropic Maxwellian distribution is the most basic equilibrium distribution:

using VelocityDistributionFunctions

v_th = 1.0 # Thermal velocity
u0 = [0.5, 0.0, 0.0] # Drift velocity
vdf = Maxwellian(v_th; u0)

# Sample from the distribution
samples = rand(vdf, 10000)

# Evaluate Theoretical PDF
𝐯 = [1.0, 0.0, 0.0]
vdf(𝐯)

0.13986259134061976

Unit is also supported

using Unitful

T = 30000u"K" # Temperature
vdf = Maxwellian(T)

𝐯 = ones(3) .* 1u"m/s"
@assert vdf(𝐯) ≈ 2.070889302986066e-19 * 1u"s^3/m^3"

# If the density is given, it will return a VelocityDistribution instead of a PDF
n = 1u"m^-3"
p = n * Unitful.k * T
vdf = Maxwellian(n, p)
vdf(𝐯)

2.070889302986066e-19 s^3 m^-6

Strip units from a distribution using ustrip.

vdf_unitless = ustrip(vdf)
vdf_unitless(ustrip(𝐯))

2.070889302986066e-19

Visualization

The following examples demonstrate sampling from distributions and comparing with theoretical PDFs using Makie.

using CairoMakie
using LinearAlgebra: norm
using VelocityDistributionFunctions: V

# Create Maxwellian distribution
vdf = Maxwellian(1.0)

# Generate samples
n_samples = 100000
vs = rand(vdf, n_samples)
speeds = norm.(vs)

# Theoretical speed distribution: f(v) = 4π v² f₃D(v)
# where f₃D is the 3D Maxwellian
v_range = range(0, 4, length=200)
speed_pdf = vdf.(V.(v_range)) # using V constructor to denote speed

# Plot
let fig = Figure()
    ax = Axis(fig[1, 1],
        xlabel="Speed",
        ylabel="Probability Density",
        title="Maxwellian Speed Distribution")

    hist!(ax, speeds, bins=100, normalization=:pdf, label="Samples")
    lines!(ax, v_range, speed_pdf, color=:red, linewidth=2, label="Theory")
    axislegend(ax, position=:rt)
    fig
end

BiMaxwellian Distribution

The BiMaxwellian distribution has different thermal velocities parallel and perpendicular to a magnetic field:

vdf = BiMaxwellian(0.5, 2.0)

vs = rand(vdf, 10000)
v_pars = getindex.(vs, 3)
v_perps = [sqrt(v[1]^2 + v[2]^2) for v in vs]

v_range = range(-4, 4, length=200)
v_par_pdf = vdf.(VPar.(v_range))
v_perp_range = range(0, 4, length=200)
v_perp_pdf = vdf.(VPerp.(v_perp_range))

let fig = Figure()
    ax = Axis(fig[1, 1], xlabel=L"v_\parallel", ylabel="Probability Density", title="Bi-Maxwellian Parallel")

    hist!(ax, v_pars, bins=100, normalization=:pdf, label="Samples")
    lines!(ax, v_range, v_par_pdf, color=:red, linewidth=2, label="Theory")

    ax2 = Axis(fig[1, 2], xlabel=L"|\mathbf{v}_\perp|", title="Bi-Maxwellian Perpendicular")
    hist!(ax2, v_perps, bins=100, normalization=:pdf, label="Samples")
    lines!(ax2, v_perp_range, v_perp_pdf, color=:red, linewidth=2, label="Theory")
    axislegend(ax2, position=:rt)
    fig
end

ShiftedPDF

Any velocity PDF can be drifted by wrapping it in ShiftedPDF.

base = MaxwellianPDF(v_th)
u0 = [0.5, 0.0, 0.0]
vdf = ShiftedPDF(base, u0)

ShiftedPDF{Float64, MaxwellianPDF{Float64}, Vector{Float64}}(
base: MaxwellianPDF{Float64}(vth=1.0)
u0: [0.5, 0.0, 0.0]
)