q(psi) from a psi map and F(psi)

[HongSik-Yun-Fusion]

The gap

Nothing in VAFT computes q(psi) from a psi map and an F(psi) profile.
Verified by searching claude/vaft-vfit-foundation for safety factor,
q_profile, def q, qpsi, q_95: everything that exists runs the other
way, or consumes a q that came from somewhere else.

• vaft/formula/equilibrium.py:335 q_from_phi(psi, phi) — needs the toroidal
  flux already.
• vaft/omas/update.py:629 update_equilibrium_profiles_1d_toroidal_flux — goes
  q -> Phi -> rho_tor, i.e. it consumes q.
• vaft/omas/update.py:270 update_equilibrium_global_quantities_q_min — reads
  profiles_1d.q.min().
• vaft/data/eqdsk.py:937,964 — q and q_95 are only ever read in from a
  g-file's QPSI.

So any ODS whose equilibrium was computed natively rather than imported from a
g-file has no route to q at all.

Proposed addition

q(psi) = sigma * F(psi) * <1/R^2> * (dV/dpsi) / (2*pi)^2

All three inputs are already outputs of
vaft.process.equilibrium.flux_surface_quantities
(vaft/process/equilibrium.py:1385): gm1 = <1/R^2>, dvolume_dpsi, and
f_profile is already an accepted keyword argument on the same call. So this is
a small, well-posed addition on top of machinery that exists and is tested.

Points to settle in review:

• Units. flux_surface_quantities returns dvolume_dpsi per radian, and
  its psi_grid contract is Wb/rad. The (2*pi)^2 above assumes that; state the
  convention explicitly and test it, since a 2π slip here is invisible in the
  profile shape and only shows up in the magnitude.
• Sign (sigma) follows the COCOS convention already resolved by
  vaft.process.equilibrium.as_equilibrium(...).convention.
• On-axis limit. flux_surface_quantities documents that its innermost one
  or two levels are the least accurate anywhere, and that volume/area are
  exactly 0 at levels_norm == 0. q_axis needs an explicit extrapolation
  policy rather than whatever falls out of a 0/0.

Why now

Three consumers, immediately:

1. VFIT's FluxGSResult.q_profile is a nan placeholder
   (vfit/flux_gs.py:769) with nothing to fill it.
2. update_equilibrium_derived_profiles chains through the toroidal-flux step,
   which consumes q — so a natively-computed equilibrium may not complete the
   chain (see VEST-Tokamak/vfit#43).
3. A q-on-a-surface reconstruction constraint (EFIT-style) needs it.

Per VFIT's CLAUDE.md responsibility split, this is a physics formula useful to any
VEST analysis code, so it belongs here rather than in VFIT.

[HongSik-Yun-Fusion]

Walkthrough: $q(\psi)$ Calculation from 2D Poloidal Flux Map and $F(\psi)$ (Issue #647)

1. Overview

We implemented direct safety factor profile $q(\psi)$ calculation from a 2D poloidal flux map $\psi(R, Z)$ and toroidal flux function $F(\psi)$, benchmarked against an analytic Solov'ev equilibrium with full COCOS compliance ($1 \leftrightarrow 11$ and sign conventions).

• Issue: #647
• PR: #680 (Merged into develop)

---

2. Changes Implemented

Formula Layer (vaft/formula/equilibrium.py)

• q_from_flux_surface_averages(gm1, dvolume_dpsi, f, cocos=None, ...):
  • Computes $q = \frac{\sigma}{(2\pi)^{2 - e_{B_p}}} F \left\langle R^{-2} \right\rangle \frac{dV}{d\psi}$
  • Supports $e_{B_p} = 0$ (per-radian, Wb/rad) and $e_{B_p} = 1$ (full Weber, Wb).
  • Complies with Sauter Eq. 23 sign conventions: $\sigma = \sigma_{Bp} \cdot \sigma_R \cdot \sigma_Z \cdot \sigma_\varphi = \sigma_{Ip} \cdot \sigma_{B0}$.

Process Layer (vaft/process/equilibrium.py)

• flux_surface_quantities(...):
  • Automatically calculates "q" across all reconstructed surfaces when f_profile is supplied.
  • Implements axis extrapolation using a degree-1 polynomial fit over innermost surfaces ($\psi_N \le 0.25$) to filter out discrete polygon vertex jitter on tiny near-axis contours.
  • Linearly interpolates any interior missing points along $\psi_N$.
• calculate_q_profile_from_psi(psi_grid, R, Z, f_profile, ...):
  • User-facing routine accepting flexible $F(\psi)$ inputs: scalar, 1D array on levels_norm, callable f(psi), or (psi_f, f_raw) tuple.
  • Automatic magnetic axis and boundary detection with custom overrides (psi_axis, psi_boundary, axis_rz, boundary).
  • Infers $\sigma_{B0} = -1$ when negative $F(\psi)$ values are supplied with default $\sigma_{B0}=1$.
  • Supports return_details=True returning dict with q, q_axis, q_95, psi_levels, and reconstructed surfaces, robust against unsorted/descending levels.

Tests & Quality Assurance (test/test_calculate_q_profile_from_psi.py, test/test_formula_catalog.py)

• Analytic Solov'ev on-axis limit: Extrapolated $q_0$ matches analytic Hessian $q_0 = \frac{|F_0|}{R_0 \sqrt{\psi_{RR} \psi_{ZZ}}}$ within 0.2% on a $129 \times 129$ grid.
• Adaptive 1D contour quadrature: Matches high-precision adaptive quadrature across surfaces within $5 \times 10^{-4}$.
• Grid convergence: Relative error drops monotonically from $65 \times 65$ to $129 \times 129$ to $257 \times 257$ grid resolutions ($< 5 \times 10^{-5}$ at 257).
• COCOS invariance & sign rules: Verified identical safety factor values under COCOS 1 and COCOS 11, with correct orientation signs under reversed current/field.
• Formula catalog & docstrings: Verified all 107 equilibrium formulas pass catalog discovery and conforming docstring validation.