Analysis: add on-demand helper for pmodb quantities (Eulerian/Lagrangian δ|B|, ∇·ξ⊥, κ·ξ⊥) instead of new HDF5 datasets

[logan-nc]

Summary

Julia GPEC has no equivalent of the Fortran pmodb_flag outputs, and it should not
get one by adding datasets to gpec.h5. Instead, add an Analysis helper that
reconstructs these quantities on demand from the response fields the file already
carries. This keeps gpec.h5 small while still letting users plot δ|B| and
∇·ξ⊥ profiles whenever they want them.

What Fortran provides

gpout_pmodb (gpec/gpout.f:2745) writes GPEC_pmod_n#.out plus these netCDF
variables, on both the (ψ, m) mode grid and the (ψ, θ) function grid:

Fortran variable	Quantity
b_eul / b_eul_fun	Eulerian δ|B| = χ₁(b_θ + q·b_ζ)/(J·B) [T]
b_lag / b_lag_fun	Lagrangian δ|B| = b_eul + ξ^ψ ∂B/∂ψ + (ξ^α/(χ₁q)) ∂B/∂θ [T]
Bdivxi_perp / _fun	B·∇·ξ⊥
Bkappaxi_perp / _fun	B·κ·ξ⊥

The Fortran docs describe this file as "necessary for NTV post processing" — its
original job was to hand δ|B| from GPEC to PENTRC across a file boundary.

Current state in Julia

The physics exists, but only as in-memory interpolants inside KineticForces:

• set_perturbation_data! (src/KineticForces/KineticForcesStructs.jl:293) builds
  kf_intr.dbob_m (δB/B) and kf_intr.divx_m (∇·ξ⊥) as CubicSeriesInterpolants
  over ψ, [npsi, mpert].
• This is a port of PENTRC's set_peq (pentrc/inputs.f90:826-877) — the JBB
  deweighting route (apply S,T,X,Y,Z in m-space → IDFT → divide by J·B² → DFT),
  see _jbb_deweight! at line 412 — not a port of gpout_pmodb, which builds
  δ|B| directly from the covariant perturbed field.
• Neither is written to HDF5. src/KineticForces/Output.jl writes only
  kinetic_forces/<method>/ (torques, energies, dTdψ, matrices, optional
  per-record bounce diagnostics).
• Because KineticForces runs in-process, the file hand-off Fortran needed is
  obsolete. The Eulerian/Lagrangian split is never formed anywhere in Julia.

So a user who wants to look at δ|B|(ψ, m) today has no supported path.

Proposal

Add an Analysis helper — e.g. compute_pmodb(h5path; ...) in
src/Analysis/PerturbedEquilibrium.jl, with a companion plot function following the
existing plot_mode_spectrogram (src/Analysis/PerturbedEquilibrium.jl:353)
conventions — that reconstructs the gpout_pmodb quantities from datasets already
present in gpec.h5.

Everything required is already in the file:

• perturbed_equilibrium/response/b_cova_theta, b_cova_zeta → b_θ, b_ζ
• perturbed_equilibrium/response/clebsch_psi → ξ^ψ
• perturbed_equilibrium/response/clebsch_alpha → ξ^α/χ₁ (already divided by χ₁)
• integration/psi, integration/q → radial axis and q
• splines/rzphi/* and the equilibrium field data → J, B, ∂B/∂ψ, ∂B/∂θ
• info/mlow, info/mpert, info/nlow, info/npert → mode metadata

Analysis.PerturbedEquilibriumModes.modes_to_theta
(src/Analysis/PerturbedEquilibriumModes.jl:41) already handles the m → θ
transform, so the helper should reuse it rather than reimplementing the DFT.

Explicit non-goal: do not add b_eul/b_lag/divxi/kappaxi datasets to
gpec.h5. These are [npsi, mpert] (and potentially [npsi, mtheta]) complex
arrays per quantity per n, which is real disk cost for a diagnostic most runs never
look at. Everything needed to regenerate them is already stored; recomputing on
demand is cheap.

Open question for the implementer

Which route to use for δB/B:

1. gpout_pmodb route — from the covariant field; gives both Eulerian and
   Lagrangian, and is exactly reproducible against the Fortran output.
2. set_peq route — reuse the existing _jbb_deweight! machinery that
   KineticForces already runs, giving the single δB/B the action integral consumes.

Route 1 is the better fit for a user-facing plotting helper (it is what the Fortran
diagnostic actually emits, and it separates Eulerian from Lagrangian). Route 2 has
the advantage of exercising code already on the NTV hot path. Cross-checking the two
against each other would itself be a useful result and would surface any discrepancy
between the two formulations.

Benchmarking requirement

Any PR addressing this issue must include a benchmark against Fortran GPEC.
The helper is only useful if its numbers match, and there is a nontrivial risk of
convention mismatch: the Julia response fields live in the working jac_type
coordinate, whereas Fortran applies gpeq_bcoordsout with jac_out/tmag_out
before writing, and the θ-space *_fun outputs additionally carry the
exp(i·n·Δφ) factor to machine toroidal angle.

Concretely, a PR should:

• Run a matched case in both codes (e.g. examples/DIIID-like_ideal_example against
  the corresponding Fortran gpec.in with pmodb_flag=t), with jac_out set so the
  coordinates genuinely correspond.
• Compare b_eul, b_lag, Bdivxi_perp, and Bkappaxi_perp in both mode space
  and θ space, over the full ψ range — not just at the edge surface.
• Report amplitude and phase agreement per m, and state the achieved tolerance.
• Include the comparison figure(s) in the PR description, following the
  benchmarks/ conventions (inputs referenced from examples/, outputs written into
  benchmarks/, nothing committed).
• Explicitly document any remaining discrepancy and its cause, rather than tuning a
  tolerance until it passes.

Note that the existing cylindrical-component helpers (xi_R, xi_Z, xi_phi,
b_R, b_Z, b_phi) are documented as beta with up to ~20% disagreement against
Fortran; this helper should not ship in that state.

Acceptance criteria

• Analysis helper returns Eulerian and Lagrangian δ|B|, B·∇·ξ⊥, and B·κ·ξ⊥ in
  mode space, with an option for θ space.
• Companion plot function following existing Analysis plotting conventions.
• No new datasets written to gpec.h5.
• Docstrings added and covered by an @autodocs block in docs/src/analysis.md.
• Fortran benchmark performed, tolerance stated, figures in the PR description.
• Regression harness run per project policy (this should be output-neutral —
  confirm it is).