Analysis Module

Analysis

Post-processing and visualization utilities for GPEC simulation outputs.

Submodules

• ForceFreeStates: Plotting functions for ForceFreeStates (DCON-style ideal MHD stability) results
• Equilibrium: Plotting functions for equilibrium profiles and flux surfaces
• CoilForcing: Plotting functions for coil forcing spectra
• PerturbedEquilibrium: Plotting functions for perturbed equilibrium and singular coupling results
• PerturbedEquilibriumModes: Data helpers to convert modal GPEC output to (ψ, θ) and (ψ, θ, φ) grids
• ErrorFields: Plotting functions for the error-field assessment (coil sensitivities, tolerance Monte Carlo, locking risk, phasing maps)

ForceFreeStates

Post-processing and visualization functions for ForceFreeStates (DCON-style ideal MHD stability) results stored in GPEC HDF5 output files.

plot_ballooning_alpha_boundaries(bnd, dpmap; save_path=nothing, psi_min=0.0)

Same diagram drawn over a heatmap of the signed Δ' from LocalStability.ballooning_delta_prime_map: each surface's Δ'(α) is oriented by the sign of its α=0 (stable) value so that positive is stable everywhere, regridded from its native physical α = αref*scale onto a shared uniform α axis, and shown with the Δ'=0 contour, the extracted boundaries, and the scan cap `maxalphascale*αexp`. Color limits are set to the 90th percentile of |Δ'| so the pole regions inside the unstable band do not wash out the marginal structure.

plot_ballooning_alpha_boundaries(bnd; save_path=nothing, psi_min=0.0)

Plot the BALOO-style infinite-n ballooning stability diagram with first and second stability boundaries: experimental α (solid black), 1st stability boundary αcrit1 (dashed red), and 2nd stability boundary αcrit2 (dashed blue) versus ψ_N. NaN entries in the boundary arrays leave natural gaps over always-stable surfaces without requiring explicit masking.

Arguments

• bnd: NamedTuple with fields psi, alpha, alpha_critical1, alpha_critical2 (as returned by LocalStability.ballooning_alpha_boundaries)

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)
• psi_min: Lower ψ_N axis limit (default: 0.0)

Returns

A Plots.jl plot object.

plot_ballooning_alpha_boundary(h5path; save_path=nothing, psi_min=0.0)

Plot the BALOO-style infinite-n ballooning stability diagram: the experimental pressure gradient α (solid) and the first stability boundary αcrit (dashed) versus normalized poloidal flux ψN. Surfaces where the experimental α lies above the boundary are ballooning-unstable. Reads LocalStability/ballooning_psi, LocalStability/alpha, and LocalStability/alpha_critical (populated when ForceFreeStates runs with local_stability_flag = true).

Arguments

• h5path: Path to a GPEC HDF5 output file (e.g. "gpec.h5")

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)
• psi_min: Lower ψ_N axis limit; set near the edge (e.g. 0.9) to focus on the pedestal (default: 0.0)

Returns

A Plots.jl plot object.

plot_cond_fbar(h5path; save_path=nothing, zoom=false)

Plot cond(F̄) vs ψ from the kinetic-singular-surface scan stored in SingularSurfaces/Kinetic/ (populated when ForceFreeStates runs with kinetic_factor > 0, singfac_min > 0).

F̄ is the kinetic Euler-Lagrange matrix formed by Schur-reducing the six kinetic matrices against the ideal A/B/C/D/E/H blocks (Logan 2015 Appendix C). Peaks in cond(F̄) locate "kinetically-displaced" singular surfaces — roots of det(F̄) that are not at ideal rational surfaces. When a peak exceeds the threshold stored in scan_threshold the ODE integrator stops there and steps across trapezoidally, mirroring Fortran ode_kin_cross.

The plot overlays the ideal rational surfaces (dotted grey, labelled with their q value) and any accepted kinetic singular surfaces (solid crimson). If no peak exceeds the threshold, kmsing = 0 and the kinetic ODE runs as a single chunk. This diagnostic is useful for anyone asking where the kinetic resonances land relative to the ideal ones.

Arguments

• h5path: Path to a GPEC HDF5 output file produced with kinetic mode enabled

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)
• zoom: If true, auto-scale the y-axis to the scan data (threshold shown as annotation only); if false (default), always include the threshold line in the y-range

Returns

A Plots.jl plot object, or nothing if no kinetic scan is stored in the file.

plot_delta_prime(h5path; save_path=nothing)

Scatter plot of Re(Δ') per singular surface vs ψN, computed from the stored asymptotic coefficients `caleftandca_right`. Points are colored red (tearing unstable, Re(Δ') > 0) or blue (tearing stable). Integer-valued q rational surfaces are annotated.

Δ' is computed as (ca_right[resnum,resnum,2,s] - ca_left[resnum,resnum,2,s]) / (4π² ψ₀), where resnum is the linear mode index of the (m,n) resonant pair at surface s.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_edge_stability_scan(h5path; save_path=nothing, ylims=(-2, 3), kwargs...)

Plot the edge stability scan energy components (et, ep, ev, evonly) vs ψ_N.

The edge scan evaluates δW_total = δW_plasma + δW_vacuum at each stored integration step in the region [psiedge, psilim], with the plasma boundary swept from psiedge to psilim. A positive et indicates stability; the truncation point is chosen at the peak et.

Four subplots are shown:

• Total energy et = ep + ev: total free-boundary energy eigenvalue
• Plasma energy ep: plasma contribution to δW
• Vacuum energy ev: vacuum (wv) contribution with singfac scaling
• Vacuum-only eigenvalue evonly: smallest eigenvalue of wv alone (no plasma response)

A horizontal dashed line at zero marks the stability boundary. A vertical dashed line marks psilim (the final truncation psi, where the peak et was found).

Arguments

• h5path: Path to a GPEC HDF5 output file produced with psiedge < psilim

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)
• ylims: y-axis limits applied to all panels (default: (-2, 3))
• kwargs...: Additional Plots.jl keyword arguments applied to all line plots (e.g. lw=2)

Returns

A Plots.jl plot object, or nothing if no EdgeScan/ group is present in the file.

plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing)

Scatter plot of energy eigenvalues vs mode index. Points are colored red (unstable, Re > 0) or green (stable, Re < 0), with a dashed reference line at zero.

Arguments

• h5path: Path to a GPEC HDF5 output file with vacuum data (vac_flag = true)

Keyword arguments

• matrix_type: Which eigenvalues to plot: :total (et), :plasma (ep), or :vacuum (ev)
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing)

Heatmap of energy eigenvector magnitudes vs (m, mode index).

Only matrix_type=:total is supported (the total energy eigenvector matrix Wₜ is stored in FreeBoundaryStability/W_freeboundary_eigenmodes). Plasma and vacuum eigenvectors are not stored separately in the HDF5 output.

Eigenvectors are scaled by χ₁ = 2π ψ₀ × 10⁻³ to match GPEC conventions.

Arguments

• h5path: Path to a GPEC HDF5 output file with vacuum data (vac_flag = true)

Keyword arguments

• matrix_type: Energy matrix to plot; only :total is currently supported
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_ffs_summary(h5path; save_path=nothing)

Four-panel summary of ForceFreeStates (DCON-style) stability results, combining:

• Energy eigenvector heatmap (plot_energy_eigenvectors)
• Fixed-boundary stability criterion |Dc| vs ψN (plot_stability_criterion)
• Eigenvalue spectrum (plot_eigenvalues)
• Tearing stability Δ' at each rational surface (plot_delta_prime)

If no vacuum data is present (vac_flag = false), only the stability criterion and Δ' panels are shown.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_fixed_boundary_stability_criterion(h5path; save_path=nothing)

Plot the stability criterion (smallest eigenvalue of W⁻¹, crit) vs ψ_N. A sign change in crit during integration indicates an ideal fixed-boundary instability.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_mode_displacement(h5path; modes=1:5, save_path=nothing)

Plot |ξψ| vs ψN for the least stable eigenmode, showing one curve per requested poloidal mode number m. The title includes the first eigenvalue dW = et[1].

Arguments

• h5path: Path to a GPEC HDF5 output file (e.g. "gpec.h5")

Keyword arguments

• modes: Iterable of m values to plot (default: 1:5)
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

Equilibrium

Post-processing and visualization functions for GPEC equilibrium objects and HDF5 outputs.

plot_equilibrium_summary(h5path; save_path=nothing)

Summary of equilibrium profiles and geometry. The (R, Z) flux surface plot occupies the full left column. Profile plots are stacked in the right column:

• q(ψ) safety factor with rational surface markers (plot_qprofile)
• μ₀p(ψ) pressure profile (plot_pressure_profile)
• 2πF(ψ) toroidal field function (plot_f_profile)

If gse.h5 is present (requires diagnose_src = true), a combined Grad-Shafranov error panel (θ slices + integrated, log scale) is appended at the bottom of the right column.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_f_profile(h5path; save_path=nothing)

Plot the toroidal field function 2πF(ψ) profile (F = RBφ/(2π)).

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing)

Plot flux surface contours (constant ψ, blue) and field-line angle spokes (constant θ, red) in physical (R, Z) space, reading nodal grid data directly from HDF5.

Psi contours are drawn at n_psi evenly spaced values between psilow and psihigh. Theta spokes are drawn at n_theta evenly spaced values.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• n_psi: Number of constant-ψ contours to draw (default: 11)
• n_theta: Number of constant-θ/θ spokes to draw (default: 18)
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing)

Plot Grad-Shafranov error vs ψN (log scale) for several θ slices, with the flux-surface-integrated error overplotted as a thick black line. Reads gse.h5 (and optionally gsei.h5) from the same directory as h5path. Returns nothing if gse.h5 is not found (requires `diagnosesrc = true` in the equilibrium configuration).

Arguments

• h5path: Path to a GPEC HDF5 output file (used to locate gse.h5)

Keyword arguments

• n_theta_lines: Number of evenly-spaced θ slices to overlay (default: 8)
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object, or nothing if gse.h5 is absent.

plot_gse_integrated(h5path; save_path=nothing)

Plot the flux-surface-integrated Grad-Shafranov error (log scale) vs ψN, reading from gsei.h5 in the same directory as h5path. Returns nothing if gsei.h5 is not found (requires `diagnosesrc = true` in the equilibrium configuration).

Arguments

• h5path: Path to a GPEC HDF5 output file (used to locate gsei.h5)

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object, or nothing if gsei.h5 is absent.

plot_pressure_profile(h5path; save_path=nothing)

Plot the μ₀p(ψ) pressure profile. Vertical dashed lines mark rational surfaces if present.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_qprofile(h5path; show_singular=true, save_path=nothing)

Plot the safety factor q(ψ) profile, with optional vertical markers at each rational surface and horizontal reference lines at q0 and q95.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• show_singular: If true, overlay rational surface locations (default: true)
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

CoilForcing

Visualization functions for coil geometry and the resulting normal field perturbations on the plasma boundary.

Two 3D entry points cover the common cases: plot_coil_geometry_3d draws the coils alone, and plot_surface_3d draws the plasma boundary surface — optionally coloured by any control-surface field (via control_surface_scalar) and/or with coils overlaid — the Julia analog of GPEC's plot_control_3d. The :RdBu colour scale is shared, so a blue (positive-current) coil and a blue (positive) field read identically.

All functions return Plots.jl plot objects and accept an optional save_path keyword argument to write the figure to disk.

_overlay_coils!(p, coil_sets; clim=nothing)

Draw each conductor of every coil set onto plot p as a 3D path3d curve, coloured on a symmetric :RdBu scale by its current (positive = blue, negative = red — the same mapping used for surface fields, so a blue coil and a blue field read identically).

control_surface_scalar(equil, modes, n; mtheta=180, nphi=180, psi=nothing,
                       helicity=nothing) -> Matrix{Float64}

Inverse-Fourier-transform a spectral control-surface quantity into a real-space [mtheta, nphi] field on the plasma boundary, ready to colour plot_surface_3d.

modes are complex Fourier amplitudes indexed by poloidal mode m for a single toroidal mode n — a Vector{ForcingMode} (filtered to n) or a (m_vals, amplitudes) tuple. The quantity is arbitrary (e.g. applied normal field b_n_x, total b_n, displacement xi_n); the caller supplies whichever modes they want to view.

The transform inverts the SFL Fourier decomposition used by ForcingTerms.fourier_decompose_bn: f(θ) = Σ_m amp_m e^{i m θ}, the per-θ correction e^{i n ν(ψ,θ)} from equil.rzphi_nu, then the toroidal extension in the machine angle c(θ,φ) = -Re[f(θ) e^{i n·hel·φ}]. The e^{i n·hel·φ} follows from the SFL→machine map φ = -hel(2π ζ + ν) (ζ the SFL toroidal angle the modes were decomposed against); substituting 2π ζ = -hel·φ - ν is what turns the decomposition's ζ exponential into the ν correction times e^{i n·hel·φ}.

Keyword arguments

• mtheta, nphi: poloidal/toroidal grid resolution
• psi: flux surface (default equil.rzphi_xs[end], the boundary)
• helicity: override sign(Bt)·sign(Ip)

plot_bn_contour(bn, mtheta, nzeta; n=nothing, save_path=nothing, kwargs...)

Plot the normal magnetic field bn(θ, ζ) on the plasma boundary as a 2D heatmap.

Expects the raw Biot-Savart B·n̂ field in Tesla (before project_normal_flux! converts to flux). If plotting post-projection data, note units are T·m² not T.

Arguments

• bn: normal field matrix [mtheta, nzeta] — typically B·n̂ in Tesla
• mtheta: number of poloidal grid points
• nzeta: number of toroidal grid points

Keyword arguments

• n: toroidal mode number (used for axis label only)
• save_path: file path to save the figure (default: nothing)
• Any extra kwargs are forwarded to Plots.heatmap

Returns

A Plots.jl plot object.

plot_coil_geometry_3d(coil_sets; save_path=nothing, kwargs...)

Plot the 3D wire geometry of all coil sets.

Each coil set is shown with a distinct color; individual conductors are plotted as lines in Cartesian (X, Y, Z) space in meters.

Arguments

• coil_sets: Vector{CoilSet} from ForcingTerms.load_coil_sets

Keyword arguments

• save_path: file path to save the figure (default: nothing)
• Any extra kwargs are forwarded to Plots.plot

Returns

A Plots.jl plot object.

plot_coil_geometry_rz(coil_sets; equil=nothing, psi=nothing, save_path=nothing, kwargs...)

Plot coil cross-sections in the (R, Z) meridional plane.

Each coil strand is shown as a closed curve in (R, Z) = (√(X²+Y²), Z) space. If equil is provided, the plasma boundary is overplotted for reference.

Arguments

• coil_sets: Vector{CoilSet} from ForcingTerms.load_coil_sets

Keyword arguments

• equil: optional PlasmaEquilibrium; if given, overlays the plasma boundary
• psi: flux surface to trace (default: equil.rzphi_xs[end], the outermost grid point)
• save_path: file path to save the figure (default: nothing)
• Any extra kwargs are forwarded to Plots.plot

Returns

A Plots.jl plot object.

plot_mode_spectrum(forcing_modes; mlow=nothing, mhigh=nothing, save_path=nothing, kwargs...)

Plot the Fourier mode spectrum |bmn| vs poloidal mode number m as a step line, grouped by n.

Arguments

• forcing_modes: Vector{ForcingMode} from compute_coil_forcing_modes!

Keyword arguments

• mlow, mhigh: poloidal mode range for the x-axis (e.g. from ForceFreeStates config); if provided, sets xlims to span the full stability range even if some modes are zero
• save_path: file path to save the figure (default: nothing)
• Any extra kwargs are forwarded to Plots.plot

Returns

A Plots.jl plot object.

plot_surface_3d(equil; color_by=nothing, coil_sets=nothing, mtheta=180, nphi=180,
                psi=nothing, clim=nothing, colorbar_title="", title="Plasma boundary",
                camera=(35, 25), save_path=nothing, kwargs...)

Plot the 3D plasma boundary surface, optionally coloured by a control-surface field and/or with coils overlaid — the Julia analog of GPEC's plot_control_3d.

Keyword arguments

• color_by: optional real-space [mtheta, nphi] scalar (any control-surface output — e.g. from control_surface_scalar). Coloured on a symmetric :RdBu scale centred at zero (positive = blue), matching the coil-current colours. nothing draws an uncoloured surface.
• coil_sets: optional Vector{CoilSet} overlaid as 3D curves (see plot_coil_geometry_3d).
• mtheta, nphi, psi: surface grid resolution and flux surface (must match color_by).
• clim: symmetric colour limit (default maximum(abs, color_by)).
• save_path: file path to save the figure; extra kwargs forward to Plots.surface.

Returns

A Plots.jl plot object.

plot_surface_3d(equil, modes, n; coil_sets=nothing, colorbar_title="", kwargs...)

Convenience method: colour the surface by control_surface_scalar(equil, modes, n; …).

PerturbedEquilibrium

Post-processing and visualization functions for GPEC perturbed equilibrium results stored in the PerturbedEquilibrium/ group of a GPEC HDF5 output file.

plot_chirikov_parameter(h5path; save_path=nothing)

Scatter plot of the Chirikov overlap parameter per rational surface vs ψ_N, with a horizontal reference line at K = 1 (island overlap threshold). Points are colored red when K > 1. Integer-valued q rational surfaces are annotated.

Requires SingularCoupling/chirikov_parameter in the HDF5 file.

Arguments

• h5path: Path to a GPEC HDF5 output file with perturbed equilibrium output

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_driven_delta_prime(h5path; save_path=nothing)

Scatter plot of Re(Δ') per rational surface vs ψN, computed by the perturbed equilibrium module (from `PerturbedEquilibrium/SingularCoupling/Deltaprime`). One marker series per toroidal mode n. Integer-valued q rational surfaces are annotated.

This is the forcing-driven Δ' (response to the applied perturbation amplitudes in intr.forcing_modes); for the equilibrium-intrinsic Δ' from the STRIDE BVP, read SingularSurfaces/Delta_prime_matrix from the HDF5 directly.

Requires PerturbedEquilibrium/SingularCoupling/Delta_prime in the HDF5 file.

Arguments

• h5path: Path to a GPEC HDF5 output file

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_island_widths(h5path; save_path=nothing)

Scatter plot of island half-width w/2 per rational surface vs ψ_N. Integer-valued q rational surfaces are annotated.

Requires SingularCoupling/island_half_width in the HDF5 file.

Arguments

• h5path: Path to a GPEC HDF5 output file with perturbed equilibrium output

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing)

Two-panel spectrogram of a perturbed equilibrium response field component:

• Top: |component| vs ψ_N, one curve per poloidal mode m. Only resonant modes (m ∈ [0, nhigh·q95)) are labeled to keep the legend readable.
• Bottom: Heatmap of |component| in (m, ψ_N) space (psi on vertical axis), with white dashed lines at rational surface locations.

Inspired by plot_spectrograms.py from OMFIT GPEC.

Arguments

• h5path: Path to a GPEC HDF5 output file with perturbed equilibrium response output

Keyword arguments

• component: Response field component to plot; one of :xi_psi, :b_psi, :b_theta, :b_zeta (default: :xi_psi). :b_psi reads the area-normalized b_psi_area_weighted dataset.
• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_perturbed_equilibrium_summary(h5path; save_path=nothing)

Three-panel composite summary of perturbed equilibrium results:

• Top-left: Resonant area-weighted field per surface (plot_resonant_area_weighted_field_amplitude)
• Top-right: Edge |b_ψ| spectrum
• Bottom: ξψ mode spectrogram (`plotmode_spectrogram`)

Arguments

• h5path: Path to a GPEC HDF5 output file with perturbed equilibrium output

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_resonant_area_weighted_field(h5path; save_path=nothing)

Five-panel summary of resonant coupling quantities at each rational surface vs ψ_N:

• |b^r|: resonant area-weighted field (plot_resonant_area_weighted_field_amplitude)
• Re(Δ'): tearing stability parameter (plot_driven_delta_prime)
• |I_res|: resonant current
• w/2: island half-width (plot_island_widths)
• K: Chirikov overlap parameter (plot_chirikov_parameter)

Inspired by plot_resonant_flux.py from OMFIT GPEC.

Arguments

• h5path: Path to a GPEC HDF5 output file with perturbed equilibrium output

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

plot_resonant_area_weighted_field_amplitude(h5path; save_path=nothing)

Scatter plot of the resonant area-weighted field |b^r| per rational surface vs ψ_N. One marker series per toroidal mode n. Integer-valued q rational surfaces are annotated. The resonant area-weighted field b^r = Φ^r/A^r is the resonant flux normalized by the scalar surface area, in tesla [Pharr 2026].

Requires the perturbed equilibrium module to have been run and SingularCoupling/resonant_area_weighted_field to be present in the HDF5 file.

Arguments

• h5path: Path to a GPEC HDF5 output file with perturbed equilibrium output

Keyword arguments

• save_path: If provided, save the figure to this path (default: nothing)

Returns

A Plots.jl plot object.

PerturbedEquilibriumModes

Post-processing utilities for converting modal GPEC outputs to theta-space and (theta, phi) grids.

modes_to_theta(h5_file, variable; mtheta=nothing, keep_sfl_phi=true)

Convert modal output (npsi, numpert_total) from a gpec.h5 file to theta-space (npsi, mtheta, npert).

Reads all required metadata (mlow, nlow, mpert, npert) and spline data (ν) from the HDF5 file.

Arguments

• h5_file::String: Path to gpec.h5 output file
• variable::String: HDF5 dataset path, e.g. "PerturbedEquilibrium/Response/xi_R"

Keyword arguments

• mtheta::Int: theta grid resolution (default: max(2*(|mlow|+mpert), 512))
• keep_sfl_phi::Bool: if true (default), output in SFL toroidal angle; if false, apply exp(i*n*ν(ψ,θ)) to convert to machine toroidal angle and conjugate if helicity > 0 (matches Fortran gpout_xbrzphifun)

Returns

• theta_data::Array{ComplexF64,3}: [npsi × mtheta × npert]
• theta_grid::Vector{Float64}: [mtheta] SFL theta ∈ [0, 1)
• n_vals::Vector{Int}: [npert] toroidal mode numbers

theta_to_thetaphi(theta_data, n_vals; nphi=nothing)

Extend theta-space data to a (θ, φ) grid via toroidal inverse DFT.

f(θ, φ) = Σₙ fₙ(θ) exp(i·n·φ)

Arguments

• theta_data::Array{ComplexF64,3}: [npsi × mtheta × npert] from modes_to_theta
• n_vals::Vector{Int}: toroidal mode numbers

Keyword arguments

• nphi::Int: toroidal grid points (default: max(4*maximum(abs.(n_vals)), 64))

Returns

• full_data::Array{ComplexF64,3}: [npsi × mtheta × nphi]
• phi_grid::Vector{Float64}: [nphi] in radians ∈ [0, 2π)

ErrorFields

Plots for the error-field assessment stored under ErrorFields/ of a GPEC HDF5 output: per-coil sensitivities, the tolerance Monte Carlo distributions, the locking risk and its dependence on tolerance, the threshold scaling, the dominant coupling mode, and coil-array phasing maps. Every plot takes a list of label => h5path pairs so design revisions can be overplotted, with a single-path convenience method for the common case.

A single unlabelled source: a gpec.h5 path, or one in-memory result.

A label paired with something to plot: a gpec.h5 path, or a result already in memory (ErrorFields.SensitivityTable, CoilSensitivities, MonteCarloResult, RiskResult, or a vector of CoilOverlap). Analysing coil geometry that was never part of a run produces the latter, so the plots accept both.

The overlap distribution lives on the Monte Carlo result and the penetration threshold on the risk result, so a plot that draws both against each other needs the pair. Pass them together as a tuple rather than separately, which would leave each half unable to find the other.

plot_applied_spectra(ctx, overlaps; normalize=true, save_path=nothing)
plot_applied_spectra(h5path, coil_sets; normalize=true, save_path=nothing, kwargs...)

Each coil set's applied spectrum against the dominant resonant mode.

Both views answer different questions and the plot shows whichever is asked for: normalize = true scales every curve to its own peak, which says whether a coil drives a different part of the spectrum; false leaves the amplitudes, which says whether it simply drives less. A coil that couples weakly for the first reason needs a geometry change; one that couples weakly for the second may just be further away.

plot_coil_sensitivities(sources; quantity=:shift, coils=nothing, yscale=:identity, save_path=nothing)
plot_coil_sensitivities(source; kwargs...)

Grouped bars of the per-coil-set dominant-mode sensitivity across runs, matched by coil set name. quantity is :shift (per millimetre of rigid in-plane shift), :tilt (per degree), :rim (the same tilt as rim displacement, the unit mechanical tolerances arrive in), :as_designed (|δ| as built), or :fraction (the resonant share of the coil's own applied spectrum, 100·|δ|·B_T0/‖b̃‖ in percent, which separates a coil that couples weakly because it is small from one that couples weakly because its spectrum lies off the dominant mode). coils restricts and orders the coil sets shown; a name a source does not carry is an error rather than an empty bar. yscale = :log10 draws the values as stems with markers instead of bars, since a bar's height on a logarithmic axis is set by the axis floor rather than by the data; non-positive values are omitted.

Each source is a gpec.h5 path, a SensitivityTable already in memory, or a vector of CoilOverlap (:as_designed and :fraction only), so a sweep of coil geometry that was never part of a run plots the same way a stored run does. :fraction needs the applied field as well as the overlap, so it is not available from a SensitivityTable.

plot_correction_requirement(req::EF.CorrectionRequirement; save_path=nothing)
plot_correction_requirement(h5path, source_name, array_names; save_path=nothing, kwargs...)

Two panels of a ErrorFields.correction_requirement: the resonant field on each rational surface before correction, after each array's dominant-mode correction alone and after the joint least-squares correction, as stems on a logarithmic axis per surface; then one polar panel per array with the complex current it needs as a phasor, magnitude as radius in kilo-ampere-turns (or as a factor on the array's current as given when any array's ampere-turns are unknown) and toroidal phase as angle: the array's solo dominant-mode factor in grey, the minimum-current member of the set's dominant-mode family in blue, and the joint least-squares factor in red. A joint factor far from the solo one, or a surface left with most of its field, says the single-mode picture of the correction is not the whole story.

plot_dominant_mode_spectrum(sources; mode=1, save_path=nothing)
plot_dominant_mode_spectrum(h5path; kwargs...)

The right singular vector of the full-window dominant coupling mode of each run, |V[m, mode]| against poloidal mode number (one series per toroidal mode), from PerturbedEquilibrium/SingularCoupling/DominantMode/.

plot_efc_ntv_limits(h5path; torque_budget, delta_threshold=nothing, safety_factor=1.0, delta_max=15, profiles=false, save_path=nothing)
plot_efc_ntv_limits(couplings::Vector{EF.EFCCoupling}; delta_threshold, torque_budget, kwargs...)

Correction current against intrinsic overlap for each correction array of a run (ErrorFields/NTV/): the linear single-mode current and the NTV-limited current whose residual torque lowers the threshold, with the largest correctable overlap marked; when the run tabulated the torques against rotation, a second panel shows those tables with the offsets found and the reference rotation, and the curve uses the torque balance (model, rotation_exponent as in efc_current_curve). The zero-rotation limits are marked with the residual torque (circle) and with the whole field's torque (star); nonzero torque_rtol or budget_rtol shade the band between the pessimistic and optimistic curves and the range of the correctable limit. delta_threshold defaults to the run's fitted penetration threshold (ErrorFields/Risk/threshold_fit); torque_budget is the torque, N·m, that would bring the reference rotation to rest. profiles adds a panel of the torque integrated from the axis, T(ψ) per kAt² at zero rotation shift for the whole field and the residual, which shows where in the plasma the torque that survives correction is deposited.

plot_error_field_summary(h5path; save_path=nothing)

Four panels of a run's error-field assessment: coil sensitivities to shift, the tolerance Monte Carlo distributions, the overlap distribution against the threshold, and the locking risk against tolerance scale. Panels whose data the run did not produce say so.

plot_linearity_check(check::EF.LinearityCheck; save_path=nothing)
plot_linearity_check(h5path; scales=(1.0, 2.0), save_path=nothing, kwargs...)

The relative error of the linear sensitivity model at finite displacements, from ErrorFields.linearity_check: one stem per row on a logarithmic axis, grouped by coil set or coherent group, with the shift and tilt rows told apart by colour, the scale by marker, and the sign by fill. The dashed line is the 1 % level the sensitivity calculation warns at for its own step residuals. A row whose linear change vanishes has no relative error and is not drawn.

plot_linearity_residuals(sources; at=:tolerance, coils=nothing, tolerances=nothing, fd_step_shift_m=nothing, fd_step_tilt_deg=nothing, yscale=:identity, save_path=nothing)
plot_linearity_residuals(source; kwargs...)

How linear each coil set's field is in its rigid motions: for the six taps (shift and tilt about x, y, z) the ratio of the central second difference ‖b̃(+h) + b̃(−h) − 2b̃(0)‖ to the largest first difference among the set's taps, as stored in ErrorFields/CoilSensitivities/*_linearity_residual. The normalization is against the set's largest linear response, so a tap whose own first difference vanishes by symmetry is judged against the terms the model keeps, and a weak tap's own curvature is understated. at = :step shows the ratio at the finite-difference steps the run used; at = :tolerance rescales each tap by tolerance / step (the second difference grows as h², the first as h), which is the curvature the linear model neglects over the coil's own tolerance range. Coherent-group amplitudes are not added. Tolerances come from the file's snapshot or tolerances; a coil without a tolerance block draws at zero. The steps come from the file's echoed deck, or from fd_step_shift_m and fd_step_tilt_deg for in-memory CoilSensitivities. The dashed line is the 1 % level at which the sensitivity calculation itself warns; on a linear axis it sets the scale and hides sub-percent residuals, so yscale = :log10 draws them as stems and markers (never bars) over the decades they span. One panel per source. For the model's error at finite displacements rather than at the step, see plot_linearity_check.

plot_locking_risk(sources; corrected=true, target_percent=nothing, save_path=nothing)
plot_locking_risk(h5path; kwargs...)

Locking probability against tolerance scale from each run's ErrorFields/Risk/ToleranceScan/, intrinsic and (when corrected) with error-field correction, batch spread as error bars, and the allowable scale at target_percent marked where the scan reaches it.

plot_needed_current(mc, array::EF.CoilOverlap; coupling=nothing, delta_threshold=nothing, torque_budget=nothing, save_path=nothing)

The distribution of the factor on array's current that the sampled machine needs to cancel its dominant-mode overlap (ErrorFields.needed_current_distribution), with the as-designed factor marked and, when the array's EFCCoupling, the threshold and the torque budget are given, the NTV-limited correctable overlap as a vertical line and the fraction beyond it (ErrorFields.uncorrectable_probability) in the legend.

plot_overlap_contributions(ctx, overlaps; save_path=nothing)
plot_overlap_contributions(h5path, coil_sets; save_path=nothing, kwargs...)

Which poloidal harmonics produced each coil set's overlap, as conj(V[m])·b̃[m] rotated so that the coil's own total is real. The bars therefore sum to its resonant fraction, printed in the legend.

This is what separates a coil that drives the resonant harmonics weakly from one that drives them strongly and cancels against itself. The single complex overlap cannot tell those apart.

plot_overlap_phasors(sources; coils=nothing, save_path=nothing)
plot_overlap_phasors(source; kwargs...)

Each coil set's as-built dominant-mode overlap δ as an arrow in the complex plane, laid head to tail in coils order so the walk ends at the run's total, drawn as one bold arrow from the origin. The picture is rotated so that the total is real: an arrow pointing along it adds to the error field, one pointing against it cancels part of another coil's, and the single number |Σδ| a Monte Carlo starts from cannot tell those apart. One panel per source, since the mode phase is arbitrary between runs. Sources are a gpec.h5 path, a SensitivityTable, or a vector of CoilOverlap; on a file the walk's end is checked against the stored Monte Carlo delta_as_designed and a mismatch is warned about.

plot_phasing_map(h5path, coil_names; psi_low=0.0, psi_high=CORE_PSI_HIGH, mode=1, nphase=180, quantity=:delta_per_kat, save_path=nothing)
plot_phasing_map(map::EF.PhasingMap; quantity=:delta_per_kat, save_path=nothing)

Contour map of the dominant-mode overlap per kilo-ampere-turn (:delta_per_kat) or of the resonant fraction of the applied field (:resonant_fraction_percent) against the relative phases of two or three coil arrays (EF.phasing_map). Two arrays give a line, three a filled contour whose axes are the phase of the middle array relative to the first and of the third relative to the middle; the extreme is marked.

plot_surface_overlay(ctx, overlaps; ntheta=256, nzeta=180, levels=4, save_path=nothing)
plot_surface_overlay(h5path, coil_sets; save_path=nothing, kwargs...)

Each coil set's applied normal field on the control surface, with contours of the dominant resonant mode over it. One panel per coil, each scaled to its own peak so the shapes are comparable.

This is the view that makes a weak overlap physical rather than numerical: a coil driving a different poloidal region from the one the dominant mode occupies couples poorly however strong it is. The poloidal angle is cut at 180° so the outboard midplane sits in the middle of each panel rather than being split across its edges.

plot_threshold_scaling(sources; save_path=nothing)
plot_threshold_scaling(h5path; kwargs...)

The sampled penetration-threshold density and P(lock|δ) of each run against its overlap distributions, on a logarithmic |δ| axis, with the fitted threshold marked.

plot_tolerance_budget(h5path; psi_low=0.0, psi_high=CORE_PSI_HIGH, mode=1, tolerance_scale=1.0, sort=:total, top=nothing, save_path=nothing)
plot_tolerance_budget(table, tolerances, coil_sets; monte_carlo=nothing, kwargs...)
plot_tolerance_budget(terms::NamedTuple; monte_carlo=nothing, sort=:total, top=nothing, save_path=nothing)

The worst-case tolerance budget of a run as stacked horizontal bars, one per coil set, coherent group and the unattributed budget: the as-designed |δ_as_designed|, the shift term and the tilt term of ErrorFields.worst_case_terms, whose sum over every bar is the abs_delta_worst_case the Monte Carlo sizes its histogram by. It shows which coil's tolerance the budget is spent on, and whether it is spent on the coil's shift, its tilt, or on the field it makes as designed. sort orders the bars by :total (largest at the top), :name, or :none (file order); top keeps only that many largest. With a Monte Carlo result (read from the file, or passed as monte_carlo) the sample mean of |δ| and the coherent as-designed total are marked, which puts the worst-case bars against what the sampling typically realises.

plot_tolerance_pdf(sources; corrected=true, normalize=false, xscale=:identity, show_batches=false, save_path=nothing)
plot_tolerance_pdf(h5path; kwargs...)

The Monte Carlo distributions of the dominant-mode overlap |δ| of each run, intrinsic and (when corrected) with error-field correction, with the as-designed overlap marked. show_batches draws each batch's intrinsic histogram faintly under the mean, which is the sampling noise a risk figure inherits, and notes in the legend the fraction of samples that fell beyond the last bin edge.