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.

BounceScratch(ntheta, mpert)

Per-surface scratch for the bounce-averaging inner loops, allocated once in compute_bounce_data and reused across every λ. Sizes are fixed for a flux surface (ntheta sub-grid points, mpert Fourier modes). Buffers the loops populate only partially are fill!-reset per λ.

Fields

• g_wb::Vector{Float64}: length ntheta — bounce-action integrand samples
• g_wd::Vector{Float64}: length ntheta — drift integrand samples
• cum_wb_arr::Vector{Float64}: length ntheta — cumulative bounce-action integral
• jvtheta::Vector{ComplexF64}: length ntheta — action integrand
• bj_samples::Vector{ComplexF64}: length ntheta — action bounce-integral samples
• wsamp::Vector{ComplexF64}: length ntheta — per-mode W bounce-integral samples
• wmu_mt::Matrix{ComplexF64}: mpert × ntheta — W_μ per θ
• wen_mt::Matrix{ComplexF64}: mpert × ntheta — W_E per θ
• expm::Vector{ComplexF64}: length mpert — Fourier basis at a θ
• pl::Vector{ComplexF64}: length ntheta — bounce phase factor
• wmu_ba::Vector{ComplexF64}: length mpert — bounce-averaged W_μ
• wen_ba::Vector{ComplexF64}: length mpert — bounce-averaged W_E
• wmats_lmda::Vector{ComplexF64}: length nqty_matrix(mpert) — packed W outer products
• tspl_f::Vector{Float64}: length 5 — in-place tspl(θ) evaluation
• int_w::Vector{Float64}, cumint_W::Matrix{Float64}: precomputed exact-cubic quadrature weights on the fixed unit θ-grid (∫ = int_w·y, cumulative = cumint_W·y); shared read-only across surfaces, see _quadrature_weights

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.

IonSpecies(; z, m, fraction=NaN, density="")

One main-ion species in a multi-ion NTV run. z/m are the charge (e) and mass (proton masses). The density is given by exactly one of fraction or density, which select a fraction of the total n_i profile or an explicit per-species profile in the kinetic file.

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"])...)

Immutable: vary a field by building a new control rather than assigning to one (the multi-species loop does this per species, and check_psi_quadrature_convergence's test builds a second control for its differing tolerance).

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
• sing_psis: Rational-surface ψ locations (sorted, from the stability analysis), used as panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection

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 "KineticForces" 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.

SurfaceBField

The periodic cubic B(θ) of one flux surface, decomposed once and reused for every λ. Holds the per-cell polynomial coefficients, B at every knot, and the stationary points of B. Consecutive stationary points bound intervals on which B is monotone, so each holds at most one bounce point; within such an interval the cached knot values locate the cell by bisection and the cell's cubic is then solved in closed form.

Fields

• knot::Vector{Float64}: cell boundaries, ascending, knot[1] = 0, knot[end] = 1
• bknot::Vector{Float64}: B at each knot
• poly::Vector{NTuple{4,Float64}}: per-cell (d, c, b, a) of S(u) = d + cu + bu² + au³, in the cell-local coordinate u = θ − knot[i]
• theta::Vector{Float64}: stationary points of B, ascending, in [0,1)
• bval::Vector{Float64}: B at each stationary point

Evaluate B at θ from the cached cell polynomials.

Bisection on the span's monotone knot values, for when the hint does not pay off.

_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.

Bounce points of v_par(θ) = 1 − (λ/bo)·B(θ) for a trapped particle. B = bo/λ at a bounce point, so a monotone interval between stationary points of B contains one iff bo/λ lies strictly between its endpoint B values — a scalar test against cached values. The crossing is then solved directly from the cell's cubic coefficients: no iteration, no spline evaluation, and no blind search, since hints carries each interval's cell over from the previous λ.

Returns roots sorted descending, the order the deepest-well and marginally-trapped logic downstream assumes; fewer than two roots signals a degenerate λ and sends the caller to the fallback. hints is sized length(bf.theta) + 1, the extra slot being the far half of the interval that wraps through the seam.

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

Cell index holding θ ∈ [0,1], from the ascending knot vector.

Solve B(θ) = btarget inside cell ic, restricted to θ ∈ [θlo, θhi]. The cubic is solved in closed form and polished with one Newton step, which recovers the digits the closed form loses when two of its roots are nearly coincident. Returns NaN if no root lies in the restricted range.

Where btarget sits relative to cell ic of a monotone span: -1 before it, 0 inside, +1 past it. s carries the span's direction so one comparison serves both. The span's first and last cells are entered part-way, at the stationary points bounding it, so their outer edge value comes from ba/bb rather than from a knot.

_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 the deepest potential well (largest midpoint v_par) among bounce-point pairs, handling pairs that wrap through θ = 0/1.

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

_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.

Cell holding the B = btarget crossing on a monotone span running from cell ia to cell ib.

λ advances monotonically through compute_bounce_data, so btarget = bo/λ falls monotonically and each interval's crossing walks steadily along the cells in one direction. Resuming from the previous λ's cell therefore costs a step or two, the same hint idiom FastInterpolations uses for its own searches. Bisection stays as the fallback for the first λ of a surface and for the sweep near a stationary point, where B is flat and the crossing can cross many cells between consecutive λ.

Solve the single B = btarget crossing on [θa, θb], a span on which B is monotone and which does not cross the θ = 0/1 seam. Returns (θ, cell); the cell feeds back as the next λ's hint. θ is NaN when no root lies in the span.

_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 roots of a·u² + b·u + c = 0 in the (count, r1, r2, r3) shape of _real_cubic_roots, r3 always NaN. A discriminant within rounding of zero is a double root rather than none; a degenerate leading coefficient falls through to the linear case.

_quadrature_weights(ntheta) → (int_w, cumint_W)

Exact integral of the CubicFit-endpoint spline on the fixed grid range(0,1,ntheta) is a constant linear functional of the node samples, so ∫ = int_w·y and the cumulative integral is cumint_W·y. The weights are obtained once per ntheta by evaluating the public integrate / cumulative_integrate! on the unit basis vectors — bit-faithful to fitting and integrating each sample vector directly, but reducing the per-λ hot loop to a dot/mul!. Cached (build guarded by a lock); the returned arrays are read-only.

Real roots of a·u³ + b·u² + c·u + d = 0 as (count, r1, r2, r3), unused slots NaN. The count comes from a conditioned ratio test rather than the sign of the discriminant, which cancels to rounding noise when two roots nearly coincide; three real roots use the trigonometric form. A degenerate leading coefficient falls through to _quadratic_real_roots.

_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}).

_resonance_nodes_from_frequencies(wbhat_f, welec_f, wdhat_f, grid; n, nl, xeval=2.5) → Vector{Float64}

Scan grid for the zeros of the resonance operator Ω_ℓ(x; ψ) = ℓ·ω_b(ψ)·√x + n·(ω_E(ψ) + ω_d(ψ)·x) at energy x = xeval, for every bounce harmonic ℓ ∈ −nl:nl, given per-ψ frequency callables. xeval defaults to 2.5 — the peak of the Maxwellian-weighted drive x^2.5·e^−x, where the resonance overlaps the most particles, so the located ψ surfaces sit on the NTV torque-density peaks (on the DIII-D case the x=2.5 nodes land ~3× closer to the measured dT/dψ spikes than the thermal-energy x=1 estimate). Roots from all harmonics are concatenated (deduplication against coincident surfaces happens in psi_panel_points).

_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.

Decompose the cubic B_vpar by walking its cells once, recording each cell's polynomial and endpoint value and solving the quadratic dS/dθ = 0 on each to get the exact stationary points. Uses only the public coeffs/CellPoly interface, so it holds for whatever θ grid the surface interpolant was built on.

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

Parallel-velocity factor v_par = 1 − (λ/bo)·B(θ) from the periodic cubic of B (B_vpar, built where the surface interpolants are constructed), keeping v_par consistent with the bounce-point roots as in Fortran's vspl.

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]

check_psi_quadrature_convergence(total, quad_err, ctrl, method)

Warn when the ψ torque quadrature terminated without satisfying the requested tolerances (hit maxevals_psi), or when a nonzero user-set atol_psi dominated termination — the silent-garbage scenario for weak applied fields, since NTV scales as δB².

combine_species_states(states) -> KineticForcesState

Sum per-species KineticForcesState results into a single total (τ = Σs τs). Per method: the scalar total_torque/total_energy are summed exactly; the dT/dψ profile is summed by linearly interpolating each species' own (psi_grid, dtdpsi) arrays onto the sorted union of the species ψ grids (zero outside a species' range), and the cumulative T(ψ) is re-integrated (trapezoid) from that summed profile. All species share the same ψ-integration range (ctrl.psilims), so the grids differ only in adaptive nodes. The interpolated/trapezoid t_cumulative is a diagnostic profile — its endpoint need not equal the exactly-summed Gauss-Kronrod total_torque, especially near sharp resonances. The combined MethodResult carries only the summed scalars and profile; per-species diagnostics (torque_profile, records, panel_psis, resonance_psis) are not aggregated and are left at their defaults.

compute_bounce_data(psi, n, l, q, bo, bmax, bmin, theta_bmax,
                    tspl, B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f,
                    divxfac, wdfac, mass, chrg, T_s, method;
                    nlmda=128, 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 ψ
• theta_bmax: θ location of Bmax (nodal knot; the passing-transit start)
• tspl: Periodic poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ]
• B_vpar: Periodic cubic of B(θ) used for v_par and the bounce-point roots (the Fortran vspl equivalent)
• 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 128, matching Fortran pentrc nlmda)
• 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, mats;
                                    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)
• mats: MatrixSplines (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 / Fortran 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).

find_sign_change_roots(f, grid) → Vector{Float64}

Locate the zeros of a callable f by scanning consecutive grid nodes for strict sign changes (f(x[i])·f(x[i+1]) < 0) and refining each bracket with Roots.Brent. Returns the refined roots in grid order (empty if f never changes sign; a node value of exactly zero is not treated as a crossing).

Single source of truth for the scan-then-Brent idiom — used for dB/dθ extrema in the bounce averaging and for kinetic-resonance surfaces in the ψ quadrature paneling.

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 Fortran 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
• psi_quad_error::Float64: Quadrature error estimate for the total torque
• panel_psis::Vector{Float64}: Quadrature panel boundaries actually used (bounds + interior resonant surfaces)
• resonance_psis::Vector{Float64}: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting

The integral is paneled at the rational-surface ψ locations (intr.sing_psis) and capped at ctrl.maxevals_psi evaluations; a warning is emitted if the quadrature fails to reach ctrl.rtol_psi/ctrl.atol_psi or if a nonzero atol_psi dominates termination.

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.

kinetic_resonance_psi_nodes(kinetic_profiles, equil; n, nl, zi=1, mi=2, electron=false, wdfac=1.0, xeval=2.5) → Vector{Float64}

ψ_N locations of kinetic resonance surfaces, for use as ψ-quadrature panel boundaries (and, per the kinetic-aware grid-packing plan, as mandatory equilibrium knots).

Locates the zeros of the trapped-branch (leff = ℓ) resonance denominator Ω(x) = leff·ω_b·√x + n·(ω_E + ω_d·x) at energy x = xeval, for every bounce harmonic ℓ ∈ −nl:nl — this is where the energy-space resonance sweeps through the drive-weighted bulk and the NTV torque density peaks (Logan & Park, Phys. Plasmas 20, 122507 (2013), §IV–V). xeval defaults to 2.5, the peak of the Maxwellian-weighted drive x^2.5·e^−x; on the DIII-D case those nodes sit ~3× closer to the measured dT/dψ spikes than the thermal x=1 estimate. The ℓ = 0 node is the ω_d-shifted ExB (superbanana-plateau) resonance.

The frequencies use the pitch-averaged large-aspect-ratio closed forms of the rlar method (tpsi! in Torque.jl / Fortran pentrc torque.F90): ω_b = (π/4)·√(ε/2)·wtran and ω_d = q·T_s/(2·ε·R₀²·Z·e·B₀)·wdfac, with ε = ⟨r⟩/⟨R⟩ clamped away from the axis where the estimate degenerates. Panel placement only needs ~peak-width accuracy, so these cheap estimates (single spline evaluations) are sufficient and no bounce averaging is performed.

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

psi_panel_points(interior, x0, xout) → Vector{Float64}

Build the ψ-quadrature node list [x0, interior points strictly inside (x0, xout), xout]. interior is the raw union of resonant-surface locations (rational surfaces ∪ kinetic resonances); this function owns the ordering: sort, drop near-duplicates (closer than PANEL_MERGE_ATOL, e.g. a kinetic resonance coinciding with a rational), and drop points within PANEL_MERGE_ATOL of a bound to avoid degenerate panels. Paneling the integral at these surfaces puts the resonant torque-density peaks (reg_spot/collisionally broadened, but narrow in ψ) on Gauss-Kronrod interval endpoints, which the rule handles natively instead of hunting them by adaptive bisection.

set_perturbation_data!(kf_intr, pe_state, ffs, 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; dVdpsi_spline=nothing,
               species_label=nothing)

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

Arguments

• h5file::HDF5.File: Open HDF5 file handle
• state::KineticForcesState: Accumulated computation results
• dVdpsi_spline: Optional dV/dψN profile interpolant; when given, dV/dψN is written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) is directly available
• species_label: nothing writes the run total to KineticForces/<method>/; a label (e.g. "ion_z1_m2", "electron") writes one species' contribution to KineticForces/PerSpecies/<label>/<method>/. A multi-species run calls this once per species and once for the summed total, so the group is opened-or-created each time.