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

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 ForceFreeStates.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 ForceFreeStates.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 locstab/psi, locstab/alpha, and locstab/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 singular/kinetic/ (populated when ForceFreeStates runs with kinetic_factor > 0, ode_flag = true, 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 perturbed_equilibrium/ 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 singular_coupling/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 singular/delta_prime_matrix from the HDF5 directly.

Requires perturbed_equilibrium/singular_coupling/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 singular_coupling/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 singular_coupling/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. "perturbed_equilibrium/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π)