diff --git a/CLAUDE.md b/CLAUDE.md
index 7d7c8c9d3..c001177fb 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -475,8 +475,10 @@ This format is used for compiling release notes, so tags should be human-readabl
 
 ### General
 - **Julia version**: 1.11 is the target version
+- **Never remove packages from Project.toml** - If a package fails to load or resolve, run `Pkg.add(...)` or `Pkg.instantiate()` to fix the local environment. Do NOT remove the package from `Project.toml`. The developer works across multiple branches and machines, so environment drift is expected — the right fix is always to update the environment to satisfy the toml, not to trim the toml to match the current environment state.
 - **Indexing**: The codebase uses 0-based indexing in many places to match Fortran conventions, then converts to 1-based Julia indexing
 - **No step numbering in code comments** - Avoid annotations like "Step 1: do this" followed by "Step 2: do that". These get out of sync as code changes. Just describe the action without numbering.
+- **Documentation coverage** - When adding a new module or submodule with public docstrings, add a corresponding `@autodocs` block in `docs/src/`. Documenter CI will fail with a `missing_docs` error if any exported docstring is not covered. The analysis submodule docs live in `docs/src/analysis.md`.
 - **Keep code comments concise** - A comment should be one line where possible. Do not write multi-line block comments explaining the current session's investigation, what was tried, what was wrong before, or why a specific file/path behaves differently. State what the code does and why at a general level. Example of too much detail: a 6-line block explaining that efit_by_inversion uses psilow>0 while CHEASE starts at 0, that the old code was removed, and that spline spikes result. Preferred: `# Replicate Fortran inverse.f: overwrite deta at axis (r²=0) by extrapolating from innermost surfaces.`
 
 ### Output Files
diff --git a/Project.toml b/Project.toml
index 65c64c3ae..9e480fa04 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,9 +15,11 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
@@ -39,9 +41,11 @@ FastGaussQuadrature = "1.1.0"
 FastInterpolations = "0.2.13"
 HDF5 = "0.17.2"
 JLD2 = "0.6.3"
+LaTeXStrings = "1.4.0"
 LinearAlgebra = "1"
 OrdinaryDiffEq = "6.102.0"
 Pkg = "1"
+PlotlyJS = "0.18.17"
 Plots = "1.40.15"
 Printf = "1"
 Roots = "2.2.13"
diff --git a/docs/src/analysis.md b/docs/src/analysis.md
index 8940a3e4b..0e7f3c100 100644
--- a/docs/src/analysis.md
+++ b/docs/src/analysis.md
@@ -6,6 +6,7 @@ The Analysis module provides post-processing and visualization utilities for GPE
 
 - `ForceFreeStates`: Plotting functions for DCON-style ideal MHD stability results
 - `Equilibrium`: Plotting functions for equilibrium objects
+- `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium results
 
 ```@docs
 GeneralizedPerturbedEquilibrium.Analysis
@@ -22,3 +23,9 @@ Modules = [GeneralizedPerturbedEquilibrium.Analysis.ForceFreeStates]
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.Analysis.Equilibrium]
 ```
+
+## PerturbedEquilibrium
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.Analysis.PerturbedEquilibrium]
+```
diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl
index 192537588..4c5a8eb7c 100644
--- a/examples/DIIID-like_ideal_example/analyze_example.jl
+++ b/examples/DIIID-like_ideal_example/analyze_example.jl
@@ -1,10 +1,16 @@
 using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
-using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS
+using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-plotlyjs()
+isinteractive() ? plotlyjs() : gr()
 
 h5path = joinpath(@__DIR__, "gpec.h5")
-p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5)
-p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path)
-p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path)
+
+# Summary plots
+p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path)
+p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path)
+p_pe  = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path)
+
+display(p_eq);  Plots.savefig(p_eq,  joinpath(@__DIR__, "equilibrium_summary.png"))
+display(p_ffs); Plots.savefig(p_ffs, joinpath(@__DIR__, "ffs_summary.png"))
+display(p_pe);  Plots.savefig(p_pe,  joinpath(@__DIR__, "pe_summary.png"))
\ No newline at end of file
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 7e8cf0216..12f073263 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -1,68 +1,68 @@
 [Equilibrium]
-eq_type = "efit"                        # Type of the input 2D equilibrium file
 eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
-jac_type = "hamada"                     # Coordinate system (hamada, pest, boozer, equal_arc)
-power_bp = 0                            # Poloidal field power exponent for Jacobian
-power_b = 0                             # Toroidal field power exponent for Jacobian
-power_r = 0                             # Major radius power exponent for Jacobian
-grid_type = "log_asymptotic"            # Radial grid packing type
-psilow = 1e-4                           # Lower limit of normalized flux coordinate
-psihigh = 0.993                         # Upper limit of normalized flux coordinate
-mpsi = 0                                # Number of radial grid points (0 = auto-compute from psi_accuracy)
-psi_accuracy = 0.001                    # Target absolute error in q for auto-mpsi
-mtheta = 256                            # Number of poloidal grid points
-newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
-etol = 1e-7                             # Error tolerance for equilibrium solver
-force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
+eq_type = "efit"               # Type of the input 2D equilibrium file
+jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_arc)
+power_bp = 0                   # Poloidal field power exponent for Jacobian
+power_b = 0                    # Toroidal field power exponent for Jacobian
+power_r = 0                    # Major radius power exponent for Jacobian
+grid_type = "log_asymptotic"   # Radial grid packing type
+psilow = 1e-4                  # Lower limit of normalized flux coordinate
+psihigh = 0.993                # Upper limit of normalized flux coordinate
+mpsi = 0                       # Number of radial grid points (0 = auto-compute from psi_accuracy)
+psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+mtheta = 256                   # Number of poloidal grid points
+newq0 = 0                      # Override for on-axis safety factor (0 = use input value)
+etol = 1e-7                    # Error tolerance for equilibrium solver
+force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
 
 [Wall]
-shape = "nowall"                        # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
-a = 0.2415                              # Distance from plasma (conformal) or shape parameter
-aw = 0.05                               # Half-thickness parameter for Dee-shaped walls
-bw = 1.5                                # Elongation parameter for wall shapes
-cw = 0                                  # Offset of wall center from major radius
-dw = 0.5                                # Triangularity parameter for wall shapes
-tw = 0.05                               # Sharpness of wall corners (try 0.05 as initial value)
-equal_arc_wall = true                   # Equal arc length distribution of nodes on wall
+shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a = 0.2415                     # Distance from plasma (conformal) or shape parameter
+aw = 0.05                      # Half-thickness parameter for Dee-shaped walls
+bw = 1.5                       # Elongation parameter for wall shapes
+cw = 0                         # Offset of wall center from major radius
+dw = 0.5                       # Triangularity parameter for wall shapes
+tw = 0.05                      # Sharpness of wall corners (try 0.05 as initial value)
+equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
-bal_flag = false              # Ideal MHD ballooning criterion for short wavelengths
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true              # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
-vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
-mer_flag = true               # Evaluate the Mercier criterian
+bal_flag = false               # Ideal MHD ballooning criterion for short wavelengths
+mat_flag = true                # Construct coefficient matrices for diagnostic purposes
+ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
+vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
+mer_flag = true                # Evaluate the Mercier criterian
 
-set_psilim_via_dmlim = true   # Safety factor (q) limit determined as q_ir+dmlim...
-dmlim = 0.2                   # See set_psilim_via_dmlim
-psiedge = 1.00                # If less then psilim, calculates dW(psi)...
-qlow = 1.02                   # Integration initiated at q determined by min(q0, qlow)...
-qhigh = 1e3                   # Integration terminated at q limit determined by min(qa, qhigh)...
-sing_start = 0                # Start integration at the sing_start'th rational from the axis (psilow)
+set_psilim_via_dmlim = true    # Safety factor (q) limit determined as q_ir+dmlim...
+dmlim = 0.2                    # See set_psilim_via_dmlim
+psiedge = 1.00                 # If less then psilim, calculates dW(psi)...
+qlow = 1.02                    # Integration initiated at q determined by min(q0, qlow)...
+qhigh = 1e3                    # Integration terminated at q limit determined by min(qa, qhigh)...
+sing_start = 0                 # Start integration at the sing_start'th rational from the axis (psilow)
 
-nn_low = 1                    # Smallest toroidal mode number to include
-nn_high = 1                   # Largest toroidal mode number to include
-delta_mlow = 8                # Expands lower bound of Fourier harmonics
-delta_mhigh = 8               # Expands upper bound of Fourier harmonics
-delta_mband = 0               # Integration keeps only this wide a band...
-mthvac = 512                  # Number of points used in splines over poloidal angle at plasma-vacuum interface.
-thmax0 = 1                    # Linear multiplier on the automatic choice of theta integration bounds
+nn_low = 1                     # Smallest toroidal mode number to include
+nn_high = 1                    # Largest toroidal mode number to include
+delta_mlow = 8                 # Expands lower bound of Fourier harmonics
+delta_mhigh = 8                # Expands upper bound of Fourier harmonics
+delta_mband = 0                # Integration keeps only this wide a band...
+mthvac = 512                   # Number of points used in splines over poloidal angle at plasma-vacuum interface.
+thmax0 = 1                     # Linear multiplier on the automatic choice of theta integration bounds
 
-kin_flag = false              # Kinetic EL equation (default: false)
-con_flag = false              # Continue integration through layers (default: false)
-kinfac1 = 1.0                 # Scale factor for energy contribution (default: 1.0)
-kinfac2 = 1.0                 # Scale factor for torque contribution (default: 1.0)
-kingridtype = 0               # Regular grid method (default: 0)
-passing_flag = true           # Includes passing particle effects (default: false)
-ktanh_flag = true             # Ignore kinetic effects in the core smoothly (default: false)
-ktc = 0.1                     # Parameter for ktanh_flag (default: 0.1)
-ktw = 50.0                    # Parameter for ktanh_flag (default: 50.0)
-ion_flag = true               # Include ion dW_k when kin_flag is true
-electron_flag = false         # Include electron dW_k when kin_flag is true
+kin_flag = false               # Kinetic EL equation (default: false)
+con_flag = false               # Continue integration through layers (default: false)
+kinfac1 = 1.0                  # Scale factor for energy contribution (default: 1.0)
+kinfac2 = 1.0                  # Scale factor for torque contribution (default: 1.0)
+kingridtype = 0                # Regular grid method (default: 0)
+passing_flag = true            # Includes passing particle effects (default: false)
+ktanh_flag = true              # Ignore kinetic effects in the core smoothly (default: false)
+ktc = 0.1                      # Parameter for ktanh_flag (default: 0.1)
+ktw = 50.0                     # Parameter for ktanh_flag (default: 50.0)
+ion_flag = true                # Include ion dW_k when kin_flag is true
+electron_flag = false          # Include electron dW_k when kin_flag is true
 
 eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
-save_interval = 3            # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
-singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
-ucrit = 1e4                   # Maximum fraction of solutions allowed before re-normalized
+save_interval = 3              # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
+singfac_min = 1e-4             # Fractional distance from rational q at which ideal jump enforced
+ucrit = 1e4                    # Maximum fraction of solutions allowed before re-normalized
 
 [ForcingTerms]
 forcing_data_file = "forcing.dat"       # Path to forcing data file (n, m, complex amplitude)
diff --git a/examples/Solovev_ideal_example/analyze_example.jl b/examples/Solovev_ideal_example/analyze_example.jl
index 192537588..4c5a8eb7c 100644
--- a/examples/Solovev_ideal_example/analyze_example.jl
+++ b/examples/Solovev_ideal_example/analyze_example.jl
@@ -1,10 +1,16 @@
 using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
-using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS
+using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-plotlyjs()
+isinteractive() ? plotlyjs() : gr()
 
 h5path = joinpath(@__DIR__, "gpec.h5")
-p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5)
-p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path)
-p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path)
+
+# Summary plots
+p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path)
+p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path)
+p_pe  = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path)
+
+display(p_eq);  Plots.savefig(p_eq,  joinpath(@__DIR__, "equilibrium_summary.png"))
+display(p_ffs); Plots.savefig(p_ffs, joinpath(@__DIR__, "ffs_summary.png"))
+display(p_pe);  Plots.savefig(p_pe,  joinpath(@__DIR__, "pe_summary.png"))
\ No newline at end of file
diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
index d9e8b6635..01961b4bc 100644
--- a/examples/Solovev_ideal_example_3D/gpec.toml
+++ b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -52,7 +52,7 @@ eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler
 singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
 ucrit = 1e3                   # Maximum fraction of solutions allowed before re-normalized
 force_wv_symmetry = true      # Forces vacuum energy matrix symmetry
-save_interval = 10            # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
+save_interval = 3            # Save every Nth ODE step (1=all, 10=every 10th). Always saves near rational surfaces.
 
 [Wall]
 shape = "conformal"                     # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
diff --git a/examples/Solovev_ideal_example_multi_n/analyze_example.jl b/examples/Solovev_ideal_example_multi_n/analyze_example.jl
index 724ab4d35..3a1e46134 100644
--- a/examples/Solovev_ideal_example_multi_n/analyze_example.jl
+++ b/examples/Solovev_ideal_example_multi_n/analyze_example.jl
@@ -1,17 +1,25 @@
 using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
-using GeneralizedPerturbedEquilibrium, Plots, PlotlyJS
+using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-plotlyjs()
+isinteractive() ? plotlyjs() : gr()
 
 # Analyze the multi-n run
 h5path = joinpath(@__DIR__, "gpec.h5")
-p_modes = Analysis.ForceFreeStates.plot_mode_displacement(h5path; modes=1:5)
-p_eigen = Analysis.ForceFreeStates.plot_eigenmode_summary(h5path)
-p_stab = Analysis.ForceFreeStates.plot_stability_criterion(h5path)
+# Summary plots
+p_eq = Analysis.Equilibrium.plot_equilibrium_summary(h5path)
+p_ffs = Analysis.ForceFreeStates.plot_ffs_summary(h5path)
+p_pe  = Analysis.PerturbedEquilibrium.plot_perturbed_equilibrium_summary(h5path)
+
+display(p_eq);  Plots.savefig(p_eq,  joinpath(@__DIR__, "equilibrium_summary.png"))
+display(p_ffs); Plots.savefig(p_ffs, joinpath(@__DIR__, "ffs_summary.png"))
+display(p_pe);  Plots.savefig(p_pe,  joinpath(@__DIR__, "pe_summary.png"))
 
 # Analyze the single-n runs
 h5path_n1 = joinpath(@__DIR__, "single_n_1", "euler_n1.h5")
 h5path_n2 = joinpath(@__DIR__, "single_n_2", "euler_n2.h5")
-p_modes_n1 = Analysis.ForceFreeStates.plot_mode_displacement(h5path_n1; modes=1:5)
-p_modes_n2 = Analysis.ForceFreeStates.plot_mode_displacement(h5path_n2; modes=1:5)
+p_n1 = Analysis.ForceFreeStates.plot_ffs_summary(h5path_n1)
+p_n2 = Analysis.ForceFreeStates.plot_ffs_summary(h5path_n2)
+
+display(p_n1); Plots.savefig(p_n1, joinpath(@__DIR__, "single_n_1", "ffs_summary.png"))
+display(p_n2); Plots.savefig(p_n2, joinpath(@__DIR__, "single_n_2", "ffs_summary.png"))
diff --git a/src/Analysis/Analysis.jl b/src/Analysis/Analysis.jl
index 1a770894e..0a2a5765f 100644
--- a/src/Analysis/Analysis.jl
+++ b/src/Analysis/Analysis.jl
@@ -5,12 +5,14 @@ Post-processing and visualization utilities for GPEC simulation outputs.
 
 ## Submodules
 
-  - `ForceFreeStates`: Plotting functions for DCON-style ideal MHD stability results
-  - `Equilibrium`: Plotting functions for equilibrium objects
+  - `ForceFreeStates`: Plotting functions for ForceFreeStates (DCON-style ideal MHD stability) results
+  - `Equilibrium`: Plotting functions for equilibrium profiles and flux surfaces
+  - `PerturbedEquilibrium`: Plotting functions for perturbed equilibrium and singular coupling results
 """
 module Analysis
 
 include("ForceFreeStates.jl")
 include("Equilibrium.jl")
+include("PerturbedEquilibrium.jl")
 
 end # module Analysis
diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl
index 19e745c8c..c82c419cb 100644
--- a/src/Analysis/Equilibrium.jl
+++ b/src/Analysis/Equilibrium.jl
@@ -1,51 +1,205 @@
 """
     Equilibrium
 
-Post-processing and visualization functions for GPEC equilibrium objects.
+Post-processing and visualization functions for GPEC equilibrium objects and HDF5 outputs.
 """
 module Equilibrium
 
+using HDF5
+using LaTeXStrings
 using Plots
 
 """
-    plot_flux_surfaces(plasma_eq; n_psi=11, n_theta=13)
+    plot_qprofile(h5path; show_singular=true, save_path=nothing)
 
-Plot flux surface contours (constant ψ, blue) and field-line angle contours (constant θ, red)
-in physical (R, Z) space.
+Plot the safety factor q(ψ) profile, with optional vertical markers at each rational surface
+and horizontal reference lines at q0 and q95.
 
 ### Arguments
 
-  - `plasma_eq`: A `GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium` object
+  - `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.
+"""
+function plot_qprofile(h5path; show_singular=true, save_path=nothing)
+    xs, q, q0, q95 = h5open(h5path, "r") do fid
+        read(fid["splines/profiles/xs"]), read(fid["splines/profiles/q"]),
+        read(fid["equil/q0"]), read(fid["equil/q95"])
+    end
+
+    p = plot(
+        xs, q;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="q",
+        title="",
+        legend=false,
+        xlims=(0, 1),
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm,
+        right_margin=8Plots.mm  # annotations at x=1.0 need room on the right
+    )
+    hline!(p, [q0, q95]; linestyle=:dot, color=:gray, label=nothing)
+    annotate!(p, 1.0, q0, text(" q0=$(round(q0, digits=2))", 7, :left, :gray))
+    annotate!(p, 1.0, q95, text(" q95=$(round(q95, digits=2))", 7, :left, :gray))
+
+    if show_singular
+        msing, psi_sing, q_sing = h5open(h5path, "r") do fid
+            read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"])
+        end
+        for s in 1:msing
+            vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing)
+            annotate!(p, psi_sing[s], q_sing[s],
+                text("  q=$(round(q_sing[s], digits=2))", 7, :left, :red))
+        end
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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.
+"""
+function plot_pressure_profile(h5path; save_path=nothing)
+    xs, mu0p, msing, psi_sing = h5open(h5path, "r") do fid
+        read(fid["splines/profiles/xs"]), read(fid["splines/profiles/mu0p"]),
+        read(fid["singular/msing"]), read(fid["singular/psi"])
+    end
+
+    p = plot(
+        xs, mu0p;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="μ₀p",
+        title="",
+        legend=false,
+        xlims=(0, 1),
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    for s in 1:msing
+        vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing)
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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.
+"""
+function plot_f_profile(h5path; save_path=nothing)
+    xs, twopif, msing, psi_sing = h5open(h5path, "r") do fid
+        read(fid["splines/profiles/xs"]), read(fid["splines/profiles/2piF"]),
+        read(fid["singular/msing"]), read(fid["singular/psi"])
+    end
+
+    p = plot(
+        xs, twopif;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="2πF",
+        title="",
+        legend=false,
+        xlims=(0, 1),
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    for s in 1:msing
+        vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing)
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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-θ contours to draw (default: 13)
+  - `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.
 """
-function plot_flux_surfaces(plasma_eq; n_psi=11, n_theta=13)
-    n_psi_grid = size(plasma_eq.rzphi_rsquared.nodal_derivs.partials, 2)
-    n_theta_grid = size(plasma_eq.rzphi_rsquared.nodal_derivs.partials, 3)
+function plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing)
+    rcoords, offset_data, xs_rz, ys_rz, ro, zo, msing, psi_sing, q_sing = h5open(h5path, "r") do fid
+        read(fid["splines/rzphi/rcoords"]), read(fid["splines/rzphi/offset"]),
+        read(fid["splines/rzphi/xs"]), read(fid["splines/rzphi/ys"]),
+        read(fid["equil/ro"]), read(fid["equil/zo"]),
+        read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"])
+    end
 
-    # Build R and Z on the full nodal grid
+    n_psi_grid = length(xs_rz)
+    n_theta_grid = length(ys_rz)
+
+    # Build R, Z on the full nodal grid
     R_grid = Matrix{Float64}(undef, n_psi_grid, n_theta_grid)
     Z_grid = Matrix{Float64}(undef, n_psi_grid, n_theta_grid)
     for ipsi in 1:n_psi_grid
-        rfac = @. sqrt(max(0.0, plasma_eq.rzphi_rsquared.nodal_derivs.partials[1, ipsi, :]))
-        angle = @. 2π * (plasma_eq.rzphi_ys + plasma_eq.rzphi_offset.nodal_derivs.partials[1, ipsi, :])
-        R_grid[ipsi, :] = plasma_eq.ro .+ rfac .* cos.(angle)
-        Z_grid[ipsi, :] = plasma_eq.zo .+ rfac .* sin.(angle)
+        for itheta in 1:n_theta_grid
+            rfac = sqrt(max(0.0, rcoords[ipsi, itheta]))
+            eta = 2π * (ys_rz[itheta] + offset_data[ipsi, itheta])
+            R_grid[ipsi, itheta] = ro + rfac * cos(eta)
+            Z_grid[ipsi, itheta] = zo + rfac * sin(eta)
+        end
     end
 
     p = plot(;
-        title="Flux Coordinate System Contours in (R, Z)",
+        title="Flux surfaces in (R, Z)",
         xlabel="R [m]",
         ylabel="Z [m]",
         aspect_ratio=:equal,
-        legend=:outertopright
+        legend=:outertopright,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
     )
 
     psi_indices = round.(Int, range(1, n_psi_grid; length=n_psi))
@@ -54,15 +208,189 @@ function plot_flux_surfaces(plasma_eq; n_psi=11, n_theta=13)
     for (i, ipsi) in enumerate(psi_indices)
         label = i == 1 ? "Constant ψ" : ""
         plot!(p, [R_grid[ipsi, :]; R_grid[ipsi, 1]], [Z_grid[ipsi, :]; Z_grid[ipsi, 1]];
-            color=:blue, linewidth=1.5, label=label)
+            color=:steelblue, linewidth=1.5, label=label)
     end
-
     for (i, itheta) in enumerate(theta_indices)
         label = i == 1 ? "Constant θ" : ""
         plot!(p, R_grid[:, itheta], Z_grid[:, itheta];
-            color=:red, linewidth=1.0, label=label)
+            color=:tomato, linewidth=0.8, label=label)
+    end
+    # Draw rational surface flux contours in red, annotated at the bottom (min Z) of each
+    for s in 1:msing
+        idx = argmin(abs.(xs_rz .- psi_sing[s]))
+        q_label = abs(q_sing[s] - round(q_sing[s])) < 0.05 ?
+            "q=$(round(Int, q_sing[s]))" : "q=$(round(q_sing[s], digits=2))"
+        plot!(p, [R_grid[idx, :]; R_grid[idx, 1]], [Z_grid[idx, :]; Z_grid[idx, 1]];
+            color=:red, linewidth=1.5,
+            label=s == 1 ? "Rational surface" : "")
+        itheta_bot = argmin(Z_grid[idx, :])
+        annotate!(p, R_grid[idx, itheta_bot], Z_grid[idx, itheta_bot],
+            text(" $q_label", 7, :left, :red))
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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 `diagnose_src = 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.
+"""
+function plot_gse_by_theta(h5path; n_theta_lines=8, save_path=nothing)
+    gse_path = joinpath(dirname(h5path), "gse.h5")
+    isfile(gse_path) || return nothing
+
+    gse_data = h5open(gse_path, "r") do fid
+        read(fid["gse_data"])  # (npsi, ntheta, 7): cols = θ, ψ, flux1, flux2, source, total, error
+    end
+
+    _, ntheta, _ = size(gse_data)
+    theta_indices = round.(Int, range(1, ntheta; length=n_theta_lines))
+
+    p = plot(;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="GSE error",
+        title="",
+        yscale=:log10,
+        xlims=(0, 1),
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    for (i, itheta) in enumerate(theta_indices)
+        psi_col = gse_data[:, itheta, 2]
+        err_col = gse_data[:, itheta, 7]
+        plot!(p, psi_col, err_col; label=i == 1 ? "Const. θ" : "", alpha=0.6)
+    end
+
+    # Overplot flux-surface-integrated error as thick reference line
+    gsei_path = joinpath(dirname(h5path), "gsei.h5")
+    if isfile(gsei_path)
+        xs_int, errori = h5open(gsei_path, "r") do fid
+            read(fid["xs"]), read(fid["errori"])
+        end
+        plot!(p, xs_int, vec(errori);
+            color=:black, linewidth=2.5, linestyle=:solid, label="integrated")
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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 `diagnose_src = 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.
+"""
+function plot_gse_integrated(h5path; save_path=nothing)
+    gsei_path = joinpath(dirname(h5path), "gsei.h5")
+    isfile(gsei_path) || return nothing
+
+    xs, errlogi = h5open(gsei_path, "r") do fid
+        read(fid["xs"]), read(fid["errlogi"])
+    end
+
+    p = plot(
+        xs, vec(errlogi);
+        xlabel="Norm. Poloidal Flux",
+        ylabel="log₁₀(integrated GSE)",
+        title="Flux-surface-integrated Grad-Shafranov error",
+        legend=false,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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.
+"""
+function plot_equilibrium_summary(h5path; save_path=nothing)
+    q0, q95, betat, betan, kappa, li1 = h5open(h5path, "r") do fid
+        read(fid["equil/q0"]), read(fid["equil/q95"]),
+        read(fid["equil/betat"]), read(fid["equil/betan"]),
+        read(fid["equil/kappa"]), read(fid["equil/li1"])
+    end
+
+    title_str = "q0=$(round(q0,digits=2))  q95=$(round(q95,digits=2))  βₜ=$(round(betat,digits=3))  βₙ=$(round(betan,digits=3))  κ=$(round(kappa,digits=2))  li1=$(round(li1,digits=3))"
+
+    p_rz   = plot_flux_surfaces(h5path)
+    p_q    = plot_qprofile(h5path; show_singular=true)
+    p_pres = plot_pressure_profile(h5path)
+    p_f    = plot_f_profile(h5path)
+    p_gse  = plot_gse_by_theta(h5path)  # includes integrated overplot; nothing if absent
+
+    # Suppress x-axis labels/ticks on all but the bottom profile panel — they share the
+    # same ψ_N axis and labeling every panel wastes vertical space.
+    hide_xaxis!(p) = plot!(p; xlabel="", xformatter=_->"", bottom_margin=1Plots.mm)
+
+    if isnothing(p_gse)
+        hide_xaxis!(p_q); hide_xaxis!(p_pres)
+        l = @layout [a{0.38w} [b; c; d]]
+        p = plot(p_rz, p_q, p_pres, p_f; layout=l, size=(1300, 750),
+            plot_title=title_str, plot_titlefontsize=10, top_margin=8Plots.mm)
+    else
+        hide_xaxis!(p_q); hide_xaxis!(p_pres); hide_xaxis!(p_f)
+        l = @layout [a{0.38w} [b; c; d; e]]
+        p = plot(p_rz, p_q, p_pres, p_f, p_gse; layout=l, size=(1300, 1000),
+            plot_title=title_str, plot_titlefontsize=10, top_margin=8Plots.mm)
     end
 
+    isnothing(save_path) || savefig(p, save_path)
     return p
 end
 
diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl
index 1fde16391..8b806a682 100644
--- a/src/Analysis/ForceFreeStates.jl
+++ b/src/Analysis/ForceFreeStates.jl
@@ -7,13 +7,14 @@ results stored in GPEC HDF5 output files.
 module ForceFreeStates
 
 using HDF5
+using LaTeXStrings
 using Plots
 
 """
     plot_mode_displacement(h5path; modes=1:5, save_path=nothing)
 
-Plot Im(ξ_ψ) vs ψ_N for the least stable eigenmode, showing one curve per requested
-poloidal mode number m.
+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
 
@@ -29,38 +30,81 @@ poloidal mode number m.
 A `Plots.jl` plot object.
 """
 function plot_mode_displacement(h5path; modes=1:5, save_path=nothing)
-    mlow, xi_psi, psi = h5open(h5path, "r") do fid
-        read(fid["info/mlow"]), read(fid["integration/xi_psi"]), read(fid["integration/psi"])
+    mlow, xi_psi, psi, et = h5open(h5path, "r") do fid
+        read(fid["info/mlow"]), read(fid["integration/xi_psi"]),
+        read(fid["integration/psi"]), read(fid["vacuum/et"])
     end
 
     mpert = size(xi_psi, 1)
     mhigh = mlow + mpert - 1
+    dW = isempty(et) ? nothing : et[1]
+    title_str = isnothing(dW) ? "Least stable mode" :
+        "Least stable mode, δW = $(round(real(dW), sigdigits=4))"
 
     p = plot(;
-        xlabel="ψ_N",
-        ylabel="Im(ξ_ψ)",
-        title="Least Stable Eigenmode ξ_ψ"
+        xlims=(0, 1),
+        xlabel="Norm. Poloidal Flux",
+        ylabel="|ξ^ψ|",
+        title=title_str,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
     )
     for m in modes
         mlow <= m <= mhigh || continue
-        plot!(p, psi, imag.(xi_psi[m-mlow+1, 1, :]); label="m=$m")
+        plot!(p, psi, abs.(xi_psi[m-mlow+1, 1, :]); label="m=$m")
     end
 
     isnothing(save_path) || savefig(p, save_path)
     return p
 end
 
+
 """
-    plot_eigenmode_summary(h5path; save_path=nothing)
+    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.
+"""
+function plot_fixed_boundary_stability_criterion(h5path; save_path=nothing)
+    psi, crit = h5open(h5path, "r") do fid
+        read(fid["integration/psi"]), read(fid["integration/crit"])
+    end
+
+    p = plot(
+        psi, crit;
+        xlims=(0, 1),
+        xlabel="Norm. Poloidal Flux",
+        ylabel="|Dᶜ|",
+        title="Fixed-Boundary Stability",
+        legend=false,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing)
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
 
-Three-panel summary of the free-boundary energy matrix eigenmodes, analogous to the
-DCON summary plot produced by OMFIT GPEC.
+"""
+    plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing)
 
-Panels:
+Heatmap of energy eigenvector magnitudes vs (m, mode index).
 
-  - Top: |eigenvector| vs index for the least stable mode
-  - Bottom-left: heatmap of |W_t| (eigenvectors) vs mode index
-  - Bottom-right: |eigenvalue| on log scale vs mode index
+Only `matrix_type=:total` is supported (the total energy eigenvector matrix `Wₜ` is stored in
+`vacuum/wt`). Plasma and vacuum eigenvectors are not stored separately in the HDF5 output.
 
 Eigenvectors are scaled by χ₁ = 2π ψ₀ × 10⁻³ to match GPEC conventions.
 
@@ -70,62 +114,104 @@ Eigenvectors are scaled by χ₁ = 2π ψ₀ × 10⁻³ to match GPEC convention
 
 ### 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.
 """
-function plot_eigenmode_summary(h5path; save_path=nothing)
-    wt, et, psio, mlow = h5open(h5path, "r") do fid
-        read(fid["vacuum/wt"]), read(fid["vacuum/et"]),
-        read(fid["equil/psio"]), read(fid["info/mlow"])
+function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing)
+    matrix_type == :total ||
+        error("matrix_type=$matrix_type not supported; only :total has eigenvector matrix stored in HDF5 (ep/ev are eigenvalue vectors, not matrices)")
+
+    wt, psio, mlow = h5open(h5path, "r") do fid
+        read(fid["vacuum/wt"]), read(fid["equil/psio"]), read(fid["info/mlow"])
     end
 
     isempty(wt) && error("No vacuum data in $h5path; rerun with vac_flag = true")
 
     chi1 = 2π * psio
-    wt = wt * (chi1 * 1e-3)
+    wt_scaled = wt * (chi1 * 1e-3)
 
-    nmn = size(wt, 1)
-    nmodes = size(wt, 2)
+    nmn = size(wt_scaled, 1)
+    nmodes = size(wt_scaled, 2)
     m_vals = (0:(nmn-1)) .+ mlow
 
-    p1 = plot(
-        m_vals, abs.(wt[:, 1]);
-        xlabel="m",
-        ylabel="|Eigenvector|",
-        title="Mode 1, |λ₁| = $(round(abs(et[1]), digits=3))",
-        legend=false
+    p = heatmap(
+        m_vals, 1:nmodes, abs.(wt_scaled');
+        xlabel="Poloidal Harmonic",
+        ylabel="Eigenmode Index",
+        title="Total Energy Eigenvectors",
+        colorbar_title="Harmonic Amplitude",
+        left_margin=10Plots.mm,
+        right_margin=20Plots.mm,
+        bottom_margin=5Plots.mm
     )
 
-    p2 = heatmap(
-        m_vals, 1:nmodes, abs.(wt');
-        xlabel="m",
-        ylabel="mode index",
-        colorbar_title="|Wₜ|"
-    )
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
 
-    p3 = scatter(
-        abs.(et), 1:nmodes;
-        xlabel="|Eigenvalue|",
-        ylabel="mode index",
-        xscale=:log10,
-        legend=false
-    )
+"""
+    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`)
 
-    l = @layout [a{0.25h}; b c{0.25w}]
-    p = plot(p1, p2, p3; layout=l, size=(900, 700))
+### 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.
+"""
+function plot_eigenvalues(h5path; matrix_type=:total, save_path=nothing)
+    dataset = Dict(:total => "vacuum/et", :plasma => "vacuum/ep", :vacuum => "vacuum/ev")
+    haskey(dataset, matrix_type) || error("matrix_type must be :total, :plasma, or :vacuum")
+
+    et = h5open(h5path, "r") do fid
+        read(fid[dataset[matrix_type]])
+    end
+
+    isempty(et) && error("No vacuum data in $h5path; rerun with vac_flag = true")
+
+    nmodes = length(et)
+    ev_real = real.(et)
+    colors = [v < 0 ? :red : :blue for v in ev_real]  # red = negative (unstable), blue = positive (stable)
+
+    p = scatter(
+        1:nmodes, ev_real;
+        xlabel="mode index",
+        ylabel="Re(eigenvalue)",
+        title="Eigenvalue spectrum ($matrix_type)",
+        legend=false,
+        color=colors,
+        markerstrokewidth=0,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    hline!(p, [0]; linestyle=:dash, color=:black, label=nothing)
 
     isnothing(save_path) || savefig(p, save_path)
     return p
 end
 
 """
-    plot_stability_criterion(h5path; save_path=nothing)
+    plot_delta_prime(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.
+Scatter plot of `Re(Δ')` per singular surface vs ψ_N, computed from the stored asymptotic
+coefficients `ca_left` and `ca_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
 
@@ -139,18 +225,96 @@ A sign change in `crit` during integration indicates an ideal fixed-boundary ins
 
 A `Plots.jl` plot object.
 """
-function plot_stability_criterion(h5path; save_path=nothing)
-    psi, crit = h5open(h5path, "r") do fid
-        read(fid["integration/psi"]), read(fid["integration/crit"])
+function plot_delta_prime(h5path; save_path=nothing)
+    msing, psi_sing, q_sing, ca_l, ca_r, psio, mn_index = h5open(h5path, "r") do fid
+        read(fid["singular/msing"]), read(fid["singular/psi"]), read(fid["singular/q"]),
+        read(fid["singular/ca_left"]), read(fid["singular/ca_right"]),
+        read(fid["equil/psio"]), read(fid["info/mn_index"])
     end
 
-    p = plot(
-        psi, crit;
-        xlabel="ψ_N",
-        ylabel="crit",
-        title="Stability criterion (smallest eigenvalue of W⁻¹) vs ψ_N",
-        legend=false
+    msing == 0 && return plot(; title="No singular surfaces found", legend=false)
+
+    numpert_total = size(ca_l, 1)
+    chi1 = 2π * psio
+
+    dp_vals = ComplexF64[]
+    for s in 1:msing
+        q_s = q_sing[s]
+        resnum = findfirst(1:numpert_total) do j
+            n_j = mn_index[j, 2]
+            n_j != 0 && abs(mn_index[j, 1] / n_j - q_s) < 1e-6
+        end
+        if isnothing(resnum)
+            push!(dp_vals, 0.0 + 0.0im)
+        else
+            dp = (ca_r[resnum, resnum, 2, s] - ca_l[resnum, resnum, 2, s]) / (2π * chi1)
+            push!(dp_vals, dp)
+        end
+    end
+
+    dp_real = real.(dp_vals)
+    colors = [v > 0 ? :red : :steelblue for v in dp_real]
+
+    p = scatter(
+        psi_sing, dp_real;
+        xlims=(0, 1),
+        xlabel="Norm. Poloidal Flux",
+        ylabel="Re(Δ')",
+        title="Tearing stability Δ'",
+        legend=false,
+        color=colors,
+        markersize=7,
+        markerstrokewidth=0,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
     )
+    hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing)
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    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 |D_c| 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.
+"""
+function plot_ffs_summary(h5path; save_path=nothing)
+    has_vac = h5open(h5path, "r") do fid
+        haskey(fid, "vacuum/wt") && !isempty(read(fid["vacuum/wt"]))
+    end
+
+    p_crit = plot_fixed_boundary_stability_criterion(h5path)
+    p_dp = plot_delta_prime(h5path)
+
+    if has_vac
+        p_evec = plot_energy_eigenvectors(h5path; matrix_type=:total)
+        p_modes = plot_mode_displacement(h5path)
+        p = plot(p_evec, p_crit, p_modes, p_dp; layout=(2, 2), size=(1100, 900), right_margin=10Plots.mm)
+    else
+        title!(p_crit, "Stability criterion (no vacuum data — rerun with vac_flag = true)")
+        p = plot(p_crit, p_dp; layout=(1, 2), size=(1100, 500), right_margin=10Plots.mm)
+    end
 
     isnothing(save_path) || savefig(p, save_path)
     return p
diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl
new file mode 100644
index 000000000..6ba24ddca
--- /dev/null
+++ b/src/Analysis/PerturbedEquilibrium.jl
@@ -0,0 +1,482 @@
+"""
+    PerturbedEquilibrium
+
+Post-processing and visualization functions for GPEC perturbed equilibrium results stored
+in the `perturbed_equilibrium/` group of a GPEC HDF5 output file.
+"""
+module PerturbedEquilibrium
+
+using HDF5
+using LaTeXStrings
+using Plots
+
+# Check that a PE dataset exists and is non-empty.
+function _has_pe_data(h5path, key)
+    h5open(h5path, "r") do fid
+        haskey(fid, key) && !isempty(read(fid[key]))
+    end
+end
+
+"""
+    plot_resonant_flux(h5path; save_path=nothing)
+
+Scatter plot of `|Φ_res|` (normalized resonant flux) per singular surface vs ψ_N. One
+marker series per toroidal mode n. Integer-valued q rational surfaces are annotated.
+
+Requires the perturbed equilibrium module to have been run and `singular_coupling/resonant_flux`
+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.
+"""
+function plot_resonant_flux(h5path; save_path=nothing)
+    key = "perturbed_equilibrium/singular_coupling/resonant_flux"
+    _has_pe_data(h5path, key) ||
+        return plot(; title="No resonant flux data — run with perturbed equilibrium enabled", legend=false)
+
+    resonant_flux, psi_sing, q_sing, msing, pe_n = h5open(h5path, "r") do fid
+        read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]),
+        read(fid["singular/msing"]),
+        read(fid["perturbed_equilibrium/forcing_modes/n"])
+    end
+
+    p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|Φ_res|",
+        title="Resonant flux |Φ_res| per surface", legend=:outertopright,
+        left_margin=10Plots.mm, bottom_margin=5Plots.mm)
+
+    n_vals = unique(pe_n)
+    for nn in n_vals
+        n_rows = findall(==(nn), pe_n)
+        rf_n = [sum(abs.(resonant_flux[n_rows, s])) for s in 1:msing]
+        scatter!(p, psi_sing, rf_n; label="n=$nn", markersize=7, markerstrokewidth=0)
+        # Annotate integer-q surfaces only (avoids overcrowding for n>1)
+        for s in 1:msing
+            abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue
+            annotate!(p, psi_sing[s], rf_n[s],
+                text("  q=$(round(Int, q_sing[s]))", 8, :left, :black))
+        end
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    plot_island_widths(h5path; save_path=nothing)
+
+Scatter plot of island half-width `w/2` per singular 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.
+"""
+function plot_island_widths(h5path; save_path=nothing)
+    key = "perturbed_equilibrium/singular_coupling/island_half_width"
+    _has_pe_data(h5path, key) ||
+        return plot(; title="No island width data — run with perturbed equilibrium enabled", legend=false)
+
+    island_hw, psi_sing, q_sing, msing = h5open(h5path, "r") do fid
+        read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]),
+        read(fid["singular/msing"])
+    end
+
+    p = scatter(
+        psi_sing, island_hw;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="w/2",
+        title="Island half-widths",
+        legend=false,
+        color=:steelblue,
+        markersize=7,
+        markerstrokewidth=0,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    for s in 1:msing
+        abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue
+        annotate!(p, psi_sing[s], island_hw[s],
+            text("  q=$(round(Int, q_sing[s]))", 8, :left, :black))
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    plot_chirikov_parameter(h5path; save_path=nothing)
+
+Scatter plot of the Chirikov overlap parameter per singular 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.
+"""
+function plot_chirikov_parameter(h5path; save_path=nothing)
+    key = "perturbed_equilibrium/singular_coupling/chirikov_parameter"
+    _has_pe_data(h5path, key) ||
+        return plot(; title="No Chirikov data — run with perturbed equilibrium enabled", legend=false)
+
+    chirikov, psi_sing, q_sing, msing = h5open(h5path, "r") do fid
+        read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]),
+        read(fid["singular/msing"])
+    end
+
+    colors = [k > 1.0 ? :red : :steelblue for k in chirikov]
+
+    p = scatter(
+        psi_sing, chirikov;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="K (Chirikov)",
+        title="Chirikov overlap parameter",
+        legend=false,
+        color=colors,
+        markersize=7,
+        markerstrokewidth=0,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    hline!(p, [1.0]; linestyle=:dash, color=:black, label=nothing)
+    for s in 1:msing
+        abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue
+        annotate!(p, psi_sing[s], chirikov[s],
+            text("  q=$(round(Int, q_sing[s]))", 8, :left, :black))
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    plot_driven_delta_prime(h5path; save_path=nothing)
+
+Scatter plot of `Re(Δ')` per singular surface vs ψ_N, computed by the perturbed equilibrium
+module (from `singular_coupling/delta_prime`). One marker series per toroidal mode n.
+Integer-valued q rational surfaces are annotated.
+
+This is complementary to `Analysis.ForceFreeStates.plot_delta_prime`, which uses the FFS
+asymptotic coefficients. The PE result includes the vacuum Green's function contribution.
+
+Requires `singular_coupling/delta_prime` 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.
+"""
+function plot_driven_delta_prime(h5path; save_path=nothing)
+    key = "perturbed_equilibrium/singular_coupling/delta_prime"
+    _has_pe_data(h5path, key) ||
+        return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false)
+
+    delta_prime, psi_sing, q_sing, msing, pe_n = h5open(h5path, "r") do fid
+        read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]),
+        read(fid["singular/msing"]),
+        read(fid["perturbed_equilibrium/forcing_modes/n"])
+    end
+
+    p = plot(; xlabel="Norm. Poloidal Flux", ylabel="Re(Δ')",
+        title="Tearing stability Δ' (PE)", legend=:outertopright,
+        left_margin=10Plots.mm, bottom_margin=5Plots.mm)
+    hline!(p, [0.0]; linestyle=:dash, color=:black, label=nothing)
+
+    n_vals = unique(pe_n)
+    for nn in n_vals
+        n_rows = findall(==(nn), pe_n)
+        dp_n = [real(delta_prime[n_rows[1], s]) for s in 1:msing]
+        colors = [v > 0 ? :red : :steelblue for v in dp_n]
+        scatter!(p, psi_sing, dp_n; label="n=$nn", color=colors,
+            markersize=7, markerstrokewidth=0)
+        for s in 1:msing
+            abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue
+            annotate!(p, psi_sing[s], dp_n[s],
+                text("  q=$(round(Int, q_sing[s]))", 8, :left, :black))
+        end
+    end
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    plot_resonant_field(h5path; save_path=nothing)
+
+Five-panel summary of resonant coupling quantities at each singular surface vs ψ_N:
+
+  - `|Φ_res|`: resonant flux (`plot_resonant_flux`)
+  - `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_field.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.
+"""
+function plot_resonant_field(h5path; save_path=nothing)
+    p1 = plot_resonant_flux(h5path)
+    p2 = plot_driven_delta_prime(h5path)
+    p3 = _plot_resonant_current(h5path)
+    p4 = plot_island_widths(h5path)
+    p5 = plot_chirikov_parameter(h5path)
+
+    # Move outer legends inside so all panels have the same right margin
+    for p in (p1, p2, p3)
+        plot!(p; title="", legend=:topright)
+    end
+    for p in (p4, p5)
+        plot!(p; title="")
+    end
+
+    hide_xaxis!(p) = plot!(p; xlabel="", xformatter=_->"", bottom_margin=1Plots.mm)
+    for p in (p1, p2, p3, p4)
+        hide_xaxis!(p)
+    end
+
+    p = plot(p1, p2, p3, p4, p5; layout=(5, 1), size=(900, 1600))
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+# Internal helper — resonant current scatter plot
+function _plot_resonant_current(h5path)
+    key = "perturbed_equilibrium/singular_coupling/resonant_current"
+    _has_pe_data(h5path, key) ||
+        return plot(; title="No resonant current data", legend=false)
+
+    resonant_current, psi_sing, q_sing, msing, pe_n = h5open(h5path, "r") do fid
+        read(fid[key]), read(fid["singular/psi"]), read(fid["singular/q"]),
+        read(fid["singular/msing"]),
+        read(fid["perturbed_equilibrium/forcing_modes/n"])
+    end
+
+    p = plot(; xlabel="Norm. Poloidal Flux", ylabel="|I_res|",
+        title="Resonant current |I_res| per surface", legend=:outertopright,
+        left_margin=10Plots.mm, bottom_margin=5Plots.mm)
+
+    n_vals = unique(pe_n)
+    for nn in n_vals
+        n_rows = findall(==(nn), pe_n)
+        rc_n = [sum(abs.(resonant_current[n_rows, s])) for s in 1:msing]
+        scatter!(p, psi_sing, rc_n; label="n=$nn", markersize=7, markerstrokewidth=0)
+        for s in 1:msing
+            abs(q_sing[s] - round(q_sing[s])) < 0.05 || continue
+            annotate!(p, psi_sing[s], rc_n[s],
+                text("  q=$(round(Int, q_sing[s]))", 8, :left, :black))
+        end
+    end
+
+    return p
+end
+
+"""
+    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`)
+  - `save_path`: If provided, save the figure to this path (default: `nothing`)
+
+### Returns
+
+A `Plots.jl` plot object.
+"""
+function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing)
+    comp_map = Dict(
+        :xi_psi  => ("xi_psi_real",  "xi_psi_imag"),
+        :b_psi   => ("b_psi_real",   "b_psi_imag"),
+        :b_theta => ("b_theta_real", "b_theta_imag"),
+        :b_zeta  => ("b_zeta_real",  "b_zeta_imag"),
+    )
+    haskey(comp_map, component) ||
+        error("component must be one of :xi_psi, :b_psi, :b_theta, :b_zeta")
+
+    rkey, ikey = comp_map[component]
+    base = "perturbed_equilibrium/response/"
+
+    _has_pe_data(h5path, base * rkey) ||
+        return plot(; title="No response data — run with perturbed equilibrium enabled", legend=false)
+
+    data_r, data_i, psi_response, mlow, mhigh, nhigh, q95, msing, psi_sing = h5open(h5path, "r") do fid
+        read(fid[base * rkey]), read(fid[base * ikey]),
+        read(fid["integration/psi"]),
+        read(fid["info/mlow"]), read(fid["info/mhigh"]), read(fid["info/nhigh"]),
+        read(fid["equil/q95"]),
+        read(fid["singular/msing"]), read(fid["singular/psi"])
+    end
+
+    data = complex.(data_r, data_i)  # shape: (npsi, numpert_total)
+    mpert = mhigh - mlow + 1
+
+    # Use only the first n's modes for a clean spectrogram (single-n assumption for display)
+    m_vals = mlow:mhigh
+    data_mn = data[:, 1:mpert]  # (npsi, mpert): first mpert columns = first n's modes
+
+    # Top panel: line plot per mode — only label resonant range m ∈ [0, nhigh·q95)
+    m_max_legend = nhigh * q95
+    p1 = plot(;
+        xlabel="Norm. Poloidal Flux",
+        ylabel="|$(component)|",
+        title="Mode spectrogram: $(component)",
+        legend=:outertopright,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    cmap = cgrad(:roma, mpert; categorical=true)
+    for (i, m) in enumerate(m_vals)
+        show_label = 0 <= m < m_max_legend
+        plot!(p1, psi_response, abs.(data_mn[:, i]);
+            label=(show_label ? "m=$m" : nothing), color=cmap[i], linewidth=1.5)
+    end
+
+    # Bottom panel: heatmap in (m, ψ_N) space — psi on vertical axis, m on horizontal
+    # z must be (n_psi, n_m) = (npsi, mpert) so that rows = psi, cols = m
+    p2 = heatmap(
+        collect(m_vals), psi_response, abs.(data_mn);
+        xlabel="m",
+        ylabel="Norm. Poloidal Flux",
+        title="",
+        colorbar_title="|$(component)|",
+        left_margin=10Plots.mm,
+        right_margin=20Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+    # Overlay rational surface locations as white dashed lines
+    for s in 1:msing
+        hline!(p2, [psi_sing[s]]; linestyle=:dash, color=:white, linewidth=1.5, label=nothing)
+    end
+
+    p = plot(p1, p2; layout=(2, 1), size=(950, 750))
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+"""
+    plot_perturbed_equilibrium_summary(h5path; save_path=nothing)
+
+Three-panel composite summary of perturbed equilibrium results:
+
+  - Top-left: Island half-widths (`plot_island_widths`)
+  - Top-right: Energy breakdown — plasma, vacuum, and total energies
+  - Bottom: ξ_ψ mode spectrogram (`plot_mode_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.
+"""
+function plot_perturbed_equilibrium_summary(h5path; save_path=nothing)
+    p_islands  = plot_resonant_flux(h5path)
+    p_bpsi     = _plot_bpsi_edge_spectrum(h5path)
+    p_spectro  = plot_mode_spectrogram(h5path; component=:xi_psi)
+
+    l = @layout [grid(1, 2){0.35h}; b]
+    p = plot(p_islands, p_bpsi, p_spectro; layout=l, size=(1100, 1100))
+
+    isnothing(save_path) || savefig(p, save_path)
+    return p
+end
+
+# Internal helper — |b_psi(m)| spectrum at the outermost psi surface
+function _plot_bpsi_edge_spectrum(h5path)
+    base = "perturbed_equilibrium/response/"
+    _has_pe_data(h5path, base * "b_psi_real") ||
+        return plot(; title="No b_psi data — run with perturbed equilibrium enabled", legend=false)
+
+    data_r, data_i, mlow, mhigh = h5open(h5path, "r") do fid
+        read(fid[base * "b_psi_real"]), read(fid[base * "b_psi_imag"]),
+        read(fid["info/mlow"]), read(fid["info/mhigh"])
+    end
+
+    mpert  = mhigh - mlow + 1
+    m_vals = mlow:mhigh
+    data   = complex.(data_r, data_i)   # (npsi, numpert_total)
+    bpsi_edge = abs.(data[end, 1:mpert])
+
+    p = bar(
+        collect(m_vals), bpsi_edge;
+        xlabel="m",
+        ylabel="|b_ψ|",
+        title="b_ψ spectrum at ψ edge",
+        legend=false,
+        color=:steelblue,
+        linewidth=0,
+        left_margin=10Plots.mm,
+        bottom_margin=5Plots.mm
+    )
+
+    return p
+end
+
+end # module PerturbedEquilibrium
diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl
index 03aea908d..652f1023b 100644
--- a/src/PerturbedEquilibrium/Utils.jl
+++ b/src/PerturbedEquilibrium/Utils.jl
@@ -61,12 +61,20 @@ perturbed_equilibrium/
 │   ├── amplitude_real # Real parts of forcing amplitudes
 │   └── amplitude_imag # Imaginary parts of forcing amplitudes
 ├── response/
-│   ├── xi_perturbed   # Displacement field
-│   └── b_perturbed    # Magnetic field perturbation
+│   ├── xi_psi_real/imag    # Radial displacement (real/imag parts)
+│   ├── b_psi_real/imag     # Normal field component
+│   ├── b_theta_real/imag   # Poloidal field component
+│   └── b_zeta_real/imag    # Toroidal field component
 ├── singular_coupling/
-│   ├── coupling_coefficient_real
-│   ├── coupling_coefficient_imag
-│   └── resonant_amplitude
+│   ├── coupling_coefficient_real/imag
+│   ├── resonant_amplitude
+│   ├── resonant_flux         # ComplexF64 [numpert_total × msing]
+│   ├── resonant_current      # ComplexF64 [numpert_total × msing]
+│   ├── island_width_sq       # ComplexF64 [numpert_total × msing]
+│   ├── penetrated_field      # ComplexF64 [numpert_total × msing]
+│   ├── delta_prime           # ComplexF64 [numpert_total × msing]; tearing stability Δ'
+│   ├── island_half_width     # Float64 [msing]; actual w/2
+│   └── chirikov_parameter    # Float64 [msing]; island overlap metric
 └── energies/
     ├── plasma_energy
     ├── vacuum_energy
@@ -113,6 +121,13 @@ function write_outputs_to_HDF5(
         coupling_group["coupling_coefficient_real"] = real(state.coupling_coefficient)
         coupling_group["coupling_coefficient_imag"] = imag(state.coupling_coefficient)
         coupling_group["resonant_amplitude"] = state.resonant_amplitude
+        coupling_group["resonant_flux"]      = state.resonant_flux
+        coupling_group["resonant_current"]   = state.resonant_current
+        coupling_group["island_width_sq"]    = state.island_width_sq
+        coupling_group["penetrated_field"]   = state.penetrated_field
+        coupling_group["delta_prime"]        = state.delta_prime
+        coupling_group["island_half_width"]  = state.island_half_width
+        coupling_group["chirikov_parameter"] = state.chirikov_parameter
 
         # Write additional metrics
         for (key, val) in intr.singular_coupling_metrics