KineticForces Module

METHOD_REGISTRY

Single source of truth for the NTV calculation methods. Each entry is a NamedTuple (name, flag, kind, doc):

• name — short method identifier used as the HDF5 group key and in intr.method
• flag — the KineticForcesControl field symbol that enables the method
• kind — dispatch routing tag consumed by method_kind / Torque.jl (:gar for the GAR/matrix family, :fcgl/:rlar/:clar for the three special-cased methods)
• doc — one-line description printed in verbose output

The method names/docs and the Compute.jl enable list are all derived from this tuple, and Torque.jl routes on kind, so the methods are enumerated in one place. To add a method: append an entry here and add the matching *_flag field to KineticForcesControl.

BounceData

Bounce-averaged quantities as functions of pitch angle λ. Produced by compute_bounce_data(), consumed by pitch integration.

EnergyIntegrationResult

Results from a single energy-space integration at one (ψ, λ, ℓ) point. Trajectory fields are only populated when ctrl.save_records=true.

EnergyParams

Parameters for the energy integrand evaluation.

KineticForcesControl

User-facing control parameters from the TOML [KineticForces] section. Configures which NTV methods to run, species parameters, tolerances, and output options.

Constructed via keyword arguments or from a TOML dict:

ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...)

KineticForcesInternal

Internal working state for KineticForces calculations. Holds equilibrium-derived quantities, profile interpolants, and integration results.

Fields replacing former module-level globals:

• ro, bo, chi1: Equilibrium geometry parameters
• mthsurf, mfac: Poloidal grid info
• dbob_m, divx_m: Perturbation mode interpolants

Equilibrium and kinetic profile data are read directly from the PlasmaEquilibrium (equil.profiles, equil.geometry) and the externally-loaded KineticProfileSplines — no shadow copies are kept on this struct.

KineticForcesInternal(equil; verbose=false)

Construct KineticForcesInternal from a PlasmaEquilibrium, extracting the equilibrium geometry parameters needed for NTV calculations.

KineticForcesState

Accumulated results from all KineticForces computations. Written to gpec.h5 under the "kinetic_forces" group.

MethodResult

Results for one NTV computation method across all flux surfaces.

PitchGARParams

Parameters for the GAR pitch-angle integrand. Uses a unified fbnce interpolant that returns [ωb, ωd, f₁, f₂, ...] at each λ, matching Fortran lambdaintgrl_lsode.

_bisect_vpar(tspl, lmda, bo, θa, θb; tol=1e-12, maxiter=50) → θ

Bisect on θ ∈ [θa, θb] to find a bounce point — the angle where the parallel velocity vanishes, v_par(θ) = 1 - (λ/bo) · B(θ) = 0. The caller is responsible for supplying a bracket where v_par(θa) and v_par(θb) have opposite signs; this routine does not check, it just halves toward the sign change.

Arguments

• tspl: 1D θ-interpolant returning at least [B(θ), …] at first index; tspl(mod(θ, 1.0))[1] extracts the local field magnitude.
• lmda: pitch-angle parameter λ = μ·bo / E.
• bo: on-axis toroidal field used to normalise λ.
• θa, θb: bracket endpoints (normalised θ ∈ [0,1)); typically straddling a B-peak.

Keyword arguments

• tol: termination tolerance — exits when either |v_par(θ_mid)| < tol or the bracket width θb − θa < tol. Default 1e-12 is tight enough that the residual v_par is dominated by the spline-evaluation roundoff of tspl.
• maxiter: hard iteration cap. Default 50 halves the bracket by ~10⁻¹⁵, well past tol on a unit-scale bracket; the cap is a runaway guard rather than the expected exit. On exhaustion the midpoint of the final bracket is returned.

Returns

• θ::Float64: the converged (or capped) bounce-point angle.

_bounce_integrate(...)

Perform bounce integrals over θ sub-grid. Computes ωbbar, ωdbar, |δJ|², and optionally W matrix outer products. Ports Fortran torque.F90 lines 674-793.

Build λ grid based on method character (f/t/p). Returns (lambda, dlambda) with endpoints excluded.

_energy_collision_frequency(x::Float64, p::EnergyParams) → Float64

Collision frequency ν(x) at normalized energy x = E/T.

• "zero": collisionless (ν = 0)
• "small": 1e-5 * we
• "krook": unmodified Krook operator
• "harmonic": (1 + (l/2)²) * krook * x^(-3/2)

_energy_integrand_real(x::Float64, p::EnergyParams) → ComplexF64

Physical energy integrand in x-space, N(x)·exp(-x)/denom(x), evaluated on the real axis. Used both by the production integral (integrate_energy via _integrate_energy_resonant) and for diagnostics (evaluate_energy_integrand).

_energy_numerator(x::Float64, p::EnergyParams) → ComplexF64

Numerator N(x) of the energy integrand, without the resonance denominator and without the Maxwellian weight exp(-x). The physical x-space integrand is N(x)·exp(-x)/denom(x); the residue at a pole uses N(xres)·exp(-xres).

For CGL there is no resonance denominator: N_cgl = x^2.5 / (i·n).

_energy_numerator_deriv(x::Float64, p::EnergyParams) → ComplexF64

x-derivative dN/dx of _energy_numerator, used for the analytic regular-part limit at a real-axis (ν=0) resonance pole. CGL is excluded — the real-pole path never carries a CGL numerator (CGL has no pole).

Find bounce points for trapped/passing particles and build θ sub-grid. Returns (t1, t2, thetapoints, thetaweights).

Find deepest potential well among bounce point pairs. Ports Fortran lines 616-639.

Full-block index (column-major) within a non-Hermitian block.

_integrate_energy_resonant(p, leff, wb, n, we, wd, atol, rtol) → ComplexF64

Energy integral in real x-space over [0, X_ENERGY_MAX], matching Fortran PENTRC's real-space integration (energy.f90). One path for any collisionality: each resonance contributes a pole x_pole = x_res − i·ν/Ω′, removed analytically by principal-value + residue (Sokhotski-Plemelj), leaving a smooth integrand for QuadGK. The collisionless case (ν ≡ 0, real-axis pole) is the exact ν→0⁺ limit of this single formula; its causal branch is carried by the signed zero of −pole_offset (see the add-back below), and its on-axis 0/0 window by the analytic regular-part limit.

_jbb_deweight!(out, jbb_modes, ft, psi, equil, mthsurf, theta_buf)

JBB deweighting step: inverse DFT → divide by J(ψ,θ)·B(ψ,θ)² → forward DFT.

Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back.

_pitch_gar_kernel_quadgk!(out::Vector{ComplexF64}, lambda, p::PitchGARParams)

In-place complex-valued kernel for quadgk!. Writes out[i] = fvals[i+2] * xint_decomposed for i in 1..nqty. QuadGK natively handles ComplexF64.

_pitch_gar_kernel_quadgk_wt!(out::Vector{ComplexF64}, lambda, p::PitchGARParams)

Dual-output pitch kernel. Fills a length-2*nqty buffer: out[1:nqty] — fwmm half: fvals * complex(0, imag(xint)) out[nqty+1:2*nqty] — ftmm half: fvals * complex(real(xint), 0)

One energy integration per λ; both halves share it.

_powspace_antideriv(x, pow)

Analytic antiderivative of |(1-x²)|^pow for pow 1-9. Matches Fortran powspace_sub cases exactly.

_real_pole_regular_part(xr, p, leff, wb, n, wd) → ComplexF64

Laurent regular part (finite limit at x → x_res) of the pole-subtracted real-axis (ν=0) integrand N(x)·exp(-x)/(i·Ω(x)) − R/(x − x_res), with Ω(x) = leff·wb·√x + n·(we+wd·x) and R = N(x_res)·exp(-x_res)/(i·Ω′).

Writing h(x) = N(x)·exp(-x), the limit is [h′(x_res) − h(x_res)·Ω″/(2Ω′)] / (i·Ω′), with h′ = (N′ − N)·exp(-x), Ω′ = leff·wb/(2√x) + n·wd, Ω″ = −leff·wb/(4·x^{3/2}).

_setup_surface_state(psi, n, l, zi, mi, wdfac, electron,
                      equil, intr, kinetic_profiles) → NamedTuple

Private helper for the calculated-matrix path. Reproduces the per-surface setup in tpsi! (theta-grid sampling, bounce-extremum finding, flux-function evaluation, diamagnetic/drift frequencies) without any of the perturbation- dependent bookkeeping or method dispatch.

This keeps compute_kinetic_matrices_at_psi! structurally independent of tpsi! so the matrix-only path can be evolved (e.g. QuadGK pitch in Phase C) without perturbing the perturbative torque pipeline.

Upper-triangle index (1 ≤ i ≤ j ≤ mpert) within a triangular block.

calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
               bmax, bmin, n_s, T_s, mass, chrg, tspl, dbob_m_f, divx_m_f,
               divxfac, wdfac)::ComplexF64

Calculate CLAR (Circular Large Aspect Ratio) torque. Uses pitch-angle resolved calculations for trapped and passing particles. Includes bounce-averaged integrals over lambda (pitch angle).

Status: Partially implemented (stub for full calculation)

calculate_fcgl(psi, n, l, tspl, dbob_m_f, divx_m_f, divxfac, n_s, T_s,
               equil, intr)::ComplexF64

Calculate FCGL (Full Circular Gyrokinetic Landau) torque. Implements simplified energy balance equation. Only valid for bounce harmonic l=0.

Based on: [Logan et al., Phys. Plasmas 2013]

calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, bmax,
              bmin, n_s, T_s, mass, chrg, tspl, dbob_m_f, divx_m_f,
              divxfac, wdfac, method, op_wmats; kwargs...)::ComplexF64

Calculate GAR (General Aspect Ratio) torque. Fully general method without aspect ratio expansion. Handles variants: FGAR (full), TGAR (trapped), PGAR (passing). Can compute torque (TMM), energy (WMM), or matrix elements (KMM/RMM).

Ports Fortran torque.F90 GAR branch (lines 529-932).

Steps

1. Compute bounce-averaged quantities via compute_bounce_data()
2. Build fbnce interpolant over λ, normalize for numerical stability
3. Integrate over pitch angle via integrate_pitch_gar_quadgk()
4. Apply torque normalization (Eq. 19, Logan et al. 2013)
5. If matrix path: assemble and normalize kinetic matrices

Keyword Arguments (rex/imx override)

• rex_override::Union{Nothing,Float64}: Override real-part multiplier for resonance operator. When both overrides are provided, bypasses method-string derivation.
• imx_override::Union{Nothing,Float64}: Override imaginary-part multiplier. Use rex_override=1.0, imx_override=1.0 to get full complex result for simultaneous kwmat/ktmat extraction via compute_kinetic_matrices_at_psi!.

Reference: [Logan et al., Phys. Plasmas 20, 122507 (2013)]

calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec, wdhat, wbhat,
               nueff, dVdpsi, n_s, T_s, dbob_m_f, bo, bmin)::ComplexF64

Calculate RLAR (Reduced Large Aspect Ratio) torque. Uses energy space integration with pitch angle averaging. Valid for low aspect ratio tokamaks (ε << 1).

Reference: [Logan et al., Phys. Plasmas, 2013]

compute_bounce_data(psi, n, l, q, bo, bmax, bmin, ibmax, theta_bmax,
                    tspl, mfac, chi1, ro, dbob_m_f, divx_m_f,
                    divxfac, wdfac, mass, chrg, T_s, method;
                    nlmda=64, ntheta=128,
                    smat=nothing, tmat=nothing, xmat=nothing,
                    ymat=nothing, zmat=nothing) → BounceData

Compute bounce-averaged quantities as functions of pitch angle λ. This is the core function that sets up all λ-dependent quantities needed by the pitch-angle quadrature.

Ports Fortran torque.F90 lines 530-816 (GAR branch).

Arguments

• psi: Normalized poloidal flux
• n: Toroidal mode number
• l: Bounce harmonic number
• q: Safety factor at this ψ
• bo: On-axis toroidal field [T]
• bmax, bmin: Max/min of B(θ) at this ψ
• ibmax: Index of Bmax in poloidal grid
• theta_bmax: θ location of Bmax
• tspl: Poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ]
• mfac: Poloidal mode numbers [mlow:mhigh]
• chi1: 2π·ψ₀ flux normalization
• ro: Major radius [m]
• dbob_m_f: δB/B Fourier modes at this ψ (ComplexF64 vector, length mpert)
• divx_m_f: ∇·ξ⊥ Fourier modes at this ψ (ComplexF64 vector, length mpert)
• divxfac, wdfac: Scaling factors
• mass: Particle mass [kg]
• chrg: Particle charge [C]
• T_s: Species temperature at this ψ [J]
• method: Method string (first char: f/t/p determines λ range)

Keyword Arguments

• nlmda: Number of pitch angle grid points (default 64)
• ntheta: Number of poloidal grid points per bounce (default 128)
• smat, tmat, xmat, ymat, zmat: Geometric matrices (mpert×mpert) for kinetic matrix path

compute_calculated_kinetic_matrices(ffs_ctrl, equil, ffs_intr, metric, ffit;
                                    kf_ctrl=KineticForcesControl(),
                                    kinetic_profiles)
    → (kw_flat, kt_flat)

Drive the KineticForces matrix kernel over the ψ grid stored in metric.xs and return (kw_flat, kt_flat) arrays of shape (mpsi, np^2, 6) matching the contract that ForceFreeStates._compute_fkg_matrices! consumes.

The arrays carry the six bounce-averaged kinetic energy / torque matrices (Logan 2015 Eqs 7.30–7.35) for every ψ on the equilibrium grid, packed as block-diagonal matrices over toroidal mode number n ∈ [nlow, nhigh] and flattened to (np = mpert·npert)².

This routine reads equilibrium-derived profiles (q, dV/dψ, ⟨r⟩, ⟨R⟩) directly from named splines on equil.profiles and equil.geometry, and kinetic profiles (n, T, ωE, ν) from the `kineticprofilesargument, avoiding the former shadow-copy pattern in KineticForcesInternal. The perturbation-mode interpolants (kfintr.dbobm,kfintr.divxm`) remain unwired and are tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope" section.

Arguments

• ffs_ctrl: ForceFreeStatesControl (carries kinetic_factor, kinetic_source)
• equil: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
• ffs_intr: ForceFreeStatesInternal (mode indexing)
• metric: MetricData (provides ψ grid via metric.xs)
• ffit: FourFitVars (used only for numpert_total cross-check)

Keyword arguments

• kf_ctrl: KineticForcesControl, defaults to KineticForcesControl(). Used to carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the KineticForces kernel needs but ForceFreeStatesControl does not expose.
• kinetic_profiles::Equilibrium.KineticProfileSplines: Required. Named kinetic- profile splines loaded via Equilibrium.load_kinetic_profiles.

Returns

• kw_flat::Array{ComplexF64,3}: Energy matrices, shape (mpsi, np^2, 6)
• kt_flat::Array{ComplexF64,3}: Torque matrices, shape (mpsi, np^2, 6)

compute_kinetic_matrices_at_psi!(kwmat, ktmat, psi, n, l, zi, mi,
    wdfac, divxfac, electron, equil, intr, kinetic_profiles)

Compute the six kinetic Euler-Lagrange coefficient matrices at a single flux surface and split them into kwmat and ktmat in the convention used by the DCON matrix-assembly path (Logan 2015 Eqs 7.30–7.35).

Rather than running two integrations like the reference Fortran PENTRC (one with rex=0, imx=1 for kwmat and another with rex=1, imx=0 for ktmat), this path integrates once with rex=imx=1 to get the full complex response, then decomposes by real/imag parts — equivalent math at half the work. After the -i/(2n) normalization inside kinetic_energy_matrices_for_euler_lagrange!, the full complex response splits cleanly:

• kwmat ← fwmm half (Fortran rex=0, imx=1 pass)
• ktmat ← ftmm half (Fortran rex=1, imx=0 pass)

Each half is complex (not pure real / pure imag), matching Fortran's two independent integration passes at torque.F90:842-847. Per-surface matrix dumps confirm element-by-element agreement with Fortran fourfit.F:1080-1082 (kwmat_l, ktmat_l). This is the Fortran convention required by the adjoint combinations kwmat ± ktmat in ForceFreeStates/Kinetic.jl / ~/Code/gpec/dcon/sing.f:967-1075 for non-Hermitian Bk, Ck, E_k.

Arguments

• kwmat::Array{ComplexF64,3}: Output (mpert×mpert×6), fwmm half, zeroed on entry
• ktmat::Array{ComplexF64,3}: Output (mpert×mpert×6), ftmm half, zeroed on entry
• psi, n, l, zi, mi, wdfac, divxfac, electron: Same as tpsi! (divxfac unused on the matrix path — retained for call-site compatibility)
• equil: PlasmaEquilibrium
• intr::KineticForcesInternal: Internal state with mode indexing, geometric matrices (smats/tmats/xmats/ymats/zmats), and per-surface θ-grid buffers
• kinetic_profiles::Equilibrium.KineticProfileSplines: Named kinetic-profile splines

Reference: [Logan et al., Phys. Plasmas 20, 122507 (2013)]

compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal,
                            ctrl::KineticForcesControl, equil, kinetic_profiles)

Calculate torque/energy for all enabled methods. For each method, integrates over flux surfaces using adaptive QuadGK quadrature via integrate_psi_quadgk. For multi-n calculations, loops over toroidal mode numbers and assembles block-diagonal kinetic matrices.

Arguments

• state::KineticForcesState: Accumulates results for all methods
• intr::KineticForcesInternal: Internal state with equilibrium data
• ctrl::KineticForcesControl: Control parameters specifying which methods to run
• equil: PlasmaEquilibrium with 2D interpolants
• kinetic_profiles::Equilibrium.KineticProfileSplines: Named kinetic-profile splines

energy_integrand_scalar(x::Float64, p::EnergyParams) → ComplexF64

Evaluate the physical energy integrand N(x)·exp(-x)/denom(x) at normalized energy x = E/T. Implements [Logan, Park, et al., Phys. Plasmas, 2013] Eq. (8).

evaluate_energy_integrand(x_grid; wn, wt, we, wd, wb, nuk, leff, n,
                           nutype="harmonic", f0type="maxwellian",
                           nufac=1.0, ximag=0.0, qt=false) → Vector{ComplexF64}

Diagnostic convenience: evaluate the physical x-space energy integrand N(x)·exp(-x)/denom(x) at specified x = E/T values. Returns the integrand value (not the integral) at each point in x_grid. Useful for plotting the energy integrand shape and verifying kinetic resonance resolution.

Example

x = 10 .^ range(-2, stop=2, length=500)
f = KineticForces.evaluate_energy_integrand(x; wn=1e3, wt=2e3, we=5e4,
        wd=1e2, wb=3e4, nuk=1e3, leff=1.0, n=1)
plot(x, real.(f); xscale=:log10, xlabel="x = E/T", ylabel="Re(integrand)")

find_resonance_energies(leff, wb, n, we, wd) → Vector{Float64}

Real positive energies x_res where the resonance condition vanishes:

Ω(x) = leff·wb·√x + n·(we + wd·x) = 0

With s = √x this is the quadratic n·wd·s² + leff·wb·s + n·we = 0. Returns the x = s² values for the positive real roots (the locations of the resonance poles of the energy integrand).

integrate_energy(wn, wt, we, wd, wb, nuk, ell, leff, n, psi, lambda, method;
                 nutype="harmonic", f0type="maxwellian", nufac=1.0,
                 ximag=0.0, qt=false, atol=1e-7, rtol=1e-5) → ComplexF64

Integrate the kinetic resonance operator over normalized energy x = E/T.

The integral ∫₀^∞ N(x)·exp(-x)/denom(x) dx is evaluated in real x-space over [0, X_ENERGY_MAX] (the integrand and its poles decay as x^p·exp(-x), so the tail there is far below any tolerance) via _integrate_energy_resonant. Each resonance pole (root of Ω(x) = leff·wb·√x + n·(we + wd·x), shifted off the real axis by collisions to xpole = xres - i·ν/Ω′) is removed by subtracting its singular part R/(x - xpole) and adding back the analytic principal-value + residue. A single formula handles all collisionalities: the collisionless case (ν ≡ 0) is the exact ν→0 limit, with its real-axis pole resolved analytically (see `integrateenergyresonant`).

Collision operator types (nutype): "zero", "small", "krook", "harmonic". Distribution function types (f0type): "maxwellian", "jkp", "cgl".

ximag is accepted for backward compatibility but no longer used — resonance poles are now handled analytically rather than by contour deformation.

Returns

• ComplexF64: energy integral value

integrate_pitch_gar_quadgk(wn, wt, we, nuk, bobmax, epsr, q, fbnce, fbnce_norm,
                            nqty, ell, n, rex, imx, psi, method; ...) → Vector{ComplexF64}

Integrate the kinetic resonance operator over pitch angle λ using adaptive Gauss-Kronrod quadrature. Uses QuadGK.quadgk! with an in-place ComplexF64 kernel buffer.

The fbnce interpolant returns [ωb, ωd, f₁, f₂, ...] at each λ, where:

• f₁ = ωb|δJ|²/ro² (scalar torque)
• f₂:end = ωb·Wouterproducts/ro² (kinetic matrix elements, if present)

Splits the domain at the trapped/passing boundary so Gauss-Kronrod resolves the kink in leff = ell + n*q (circulating) → ell (trapped). One quadgk! call writes all nqty complex quantities per λ-evaluation.

Returns

• Vector{ComplexF64} of length nqty: integrated pitch-angle results

integrate_pitch_gar_quadgk_wt(wn, wt, we, nuk, bobmax, epsr, q, fbnce, fbnce_norm,
                               nqty, ell, n, psi, method; ...) → Vector{ComplexF64}

Dual-output variant for the kinetic-matrix path. Emits both the wmm half (rex=0, imx=1 → Fortran kwmat) and the tmm half (rex=1, imx=0 → Fortran ktmat) in a single pitch integration, sharing one energy integration per (λ, E).

Returns a length-2*nqty packed buffer: [wmm | tmm]. The two halves each reproduce Fortran's independent-pass result at Fortran's element-by-element convention (verified via matrix-dump comparison vs ~/Code/gpec/dcon/fourfit.F kwmat_l/ktmat_l). Downstream kwmat ± ktmat combinations in ForceFreeStates/Kinetic.jl then reproduce sing.f:967-1075 exactly for the non-Hermitian Bk, Ck, E_k diagonals.

integrate_psi_quadgk(n, nl, zi, mi, wdfac, divxfac, electron, method,
                      equil, intr, ctrl, kinetic_profiles; psi_min, psi_max) → NamedTuple

Integrate torque over ψ using adaptive Gauss-Kronrod quadrature with QuadGK.BatchIntegrand. Every integrand evaluation is logged, giving a diagnostic T(ψ) profile at no extra cost (the values are computed anyway — we just keep them).

Returns

NamedTuple with:

• total::ComplexF64: Total integrated torque
• torque_profile: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points
• matrix_integrated: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method)
• psi_nsteps::Int: Number of integrand evaluations

kinetic_energy_matrices_for_euler_lagrange!(kwmat, ktmat, state, psi, n, l, wdfac, intr;
                                            kwargs...) → nothing

Compute the six kinetic Euler-Lagrange coefficient matrices of Logan 2015 Eqs 7.30–7.35 (Ak, Bk, Ck, Dk, Ek, Hk) at a single (ψ, n, ℓ) and write them into pre-allocated kwmat[mpert, mpert, 6] (fwmm half, Fortran rex=0, imx=1) and ktmat[mpert, mpert, 6] (ftmm half, rex=1, imx=0).

Matrix-only path (no scalar torque slot in the pitch-angle buffer), so nqty = mpert²·6 instead of 1 + mpert²·6.

Uses integrate_pitch_gar_quadgk_wt to emit both halves from a single energy integration per (λ, E), reproducing Fortran's two-pass semantics (torque.F90:842-847). Per-surface matrix dumps confirm element-by-element match against Fortran fourfit.F:1080-1082 (kwmat_l, ktmat_l).

For the Hermitian-outer-product blocks A/D/H stored as upper-triangles, the mirror rule differs between halves: kwmat[j,i] = conj(kwmat[i,j]) — Hermitian (Sw pure imaginary) ktmat[j,i] = -conj(ktmat[i,j]) — anti-Hermitian (St pure real) Derivation: conj(S_w) = -S_w vs conj(S_t) = S_t, combined with conj(factor) = -factor (factor = -i/(2n)). These mirrors recover Fortran's independent-slot computation at the mirrored (j,i) positions.

method_kind(name) -> Symbol

Return the dispatch routing tag (:gar, :fcgl, :rlar, :clar) for the NTV method name, looked up from [METHOD_REGISTRY]. Errors on an unknown method.

Number of packed complex entries per λ for the 6 kinetic matrices.

powspace(xmin, xmax, pow, num, endpoints) → (points, weights)

Generate a grid with power-law concentration near endpoints. Port of Fortran powspace_sub from equil/grid.f90.

Arguments

• xmin, xmax: Grid bounds
• pow::Int: Power of grid concentration (higher = more refined near edges)
• num::Int: Number of grid points
• endpoints::String: Where to concentrate: "lower", "upper", or "both"

Returns

• points::Vector{Float64}: Grid point locations
• weights::Vector{Float64}: Derivatives dx/dnorm (integration weights)

print_summary(state::KineticForcesState; verbose::Bool=false)

Print a summary of KineticForces results to stdout.

Arguments

• state::KineticForcesState: Accumulated computation results
• verbose::Bool: Print detailed per-surface results

set_perturbation_data!(kf_intr, pe_state, ffs_intr, equil, metric)

Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal.

Builds three interpolant sets from PE Clebsch displacements:

1. xs_m — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ
2. dbob_m — δB/B Fourier modes via JBB deweighting (Fortran set_peq)
3. divx_m — ∇·ξ⊥ Fourier modes via JBB deweighting

The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868):

1. Apply geometric matrices S,T,X,Y,Z in m-space
2. Inverse DFT to θ-space
3. Divide by J·B² at each θ
4. Forward DFT back to m-space

tpsi!(tpsi_var, psi, n, l, zi, mi, wdfac, divxfac, electron, method, equil, intr,
      kinetic_profiles; op_wmats=nothing)

Toroidal torque resulting from nonambipolar transport in perturbed equilibrium. Imaginary component is proportional to the kinetic energy Im(T) = 2ndW_k.

Arguments

• tpsi_var: Output complex torque value
• psi::Float64: Normalized poloidal flux
• n::Int: Toroidal mode number
• l::Int: Bounce harmonic number
• zi::Int: Ion charge in fundamental units (e)
• mi::Int: Ion mass (units of proton mass)
• wdfac::Float64: Drift factor
• divxfac::Float64: Divergence factor
• electron::Bool: Calculate quantities for electrons (zi,mi ignored)
• method::String: Integration method (RLAR, CLAR, *GAR, *TMM, *WMM, *KMM) where * = F,T,P for full,trapped,passing
• equil: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
• intr::KineticForcesInternal: Internal state with mode indexing and perturbation splines
• kinetic_profiles::Equilibrium.KineticProfileSplines: Named kinetic-profile splines (ni, ne, Ti, Te, ωE, νi, ν_e) loaded from kinetic.dat

Optional Arguments

• op_wmats::Array{ComplexF64,3}: Store ForceFreeStates matrix elements

Returns

• ComplexF64: Toroidal torque due to nonambipolar transport

write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegrationResult})

Write variable-length integration trajectory records using offset-indexed concatenated arrays. This is the standard HDF5 ragged array pattern for storing variable-length data.

Arguments

• mg::HDF5.Group: HDF5 group for this method
• records::Vector{EnergyIntegrationResult}: Integration records to write

write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState)

Write all KineticForces results to the "kinetic_forces" group in gpec.h5.

Arguments

• h5file::HDF5.File: Open HDF5 file handle
• state::KineticForcesState: Accumulated computation results