FFS - IMPROVEMENT - Ground the density cap in the Frobenius scale and protect rationals

Reframes the cap as log-uniform sampling of the Frobenius region: near the axis every
component is a power law in psi, and cubic interpolation of psi^p on a log-uniform
grid errs by ~(p*c)^4/384, so c = 0.05 resolves even the steepest spectrum component
(p = mmax/2) to ~2e-4 while the physics responds far below that.

Two generalization guards, motivated by cross-equilibrium testing: the capped region
now ends at the innermost rational surface when that sits inside psi = 0.1, and no
knot inside a rational's RATIONAL_RES_RADIUS window is ever removed -- preserving the
Delta'-stencil structure for decks (e.g. higher n) whose rationals reach the core.
Both guards are no-ops on every current case, verified: DIII-D m512 and Solovev m512
reproduce the previous cap's step counts and knot sets exactly.

Cross-equilibrium check of the rule itself: tj_analytic_direct m1024 (analytic,
traced) 1470 -> 959 steps at et[1] 1.3e-9; LAR m1024 (inversion path, clean geometry)
951 -> 875 at et[1] identical to 8 digits; Solovev ldp m512/m1024 ~unchanged steps at
~4.5e-7 absolute et[1] shift (a +-10.4 cancellation amplifies this to 3e-5 relative).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: logan-nc

src/ForceFreeStates/Fourfit.jl

@@ -499,11 +499,22 @@ function make_matrix(equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStates
     # Decouple the coefficient-spline knots from the equilibrium grid in the packed core:
     # a cubic spline's third-derivative jumps scale as (node error)/dpsi^3, so equilibrium-grade
     # core packing amplifies tolerance-level node error into jumps that slave the EL step size.
-    # The coefficients are near-cylindrical there; cap the density at dpsi >= 0.05*psi below psi=0.1.
+    # Near the axis every component is a Frobenius power law in psi, and power laws are scale-free,
+    # so log-uniform sampling (dpsi >= c*psi) resolves them at constant relative accuracy: cubic
+    # interpolation of psi^p on that grid errs by ~(p*c)^4/384, so c = 0.05 resolves even the
+    # steepest spectrum component (p = mmax/2) to ~2e-4 while the physics responds far below that
+    # (the steep components carry vanishing solution amplitude). The capped region ends at the
+    # innermost rational surface (or psi = 0.1, whichever is smaller) and never removes a knot
+    # inside a rational's resolution window, preserving the Delta'-stencil structure the
+    # equilibrium grid encodes (GridRefinement RATIONAL_RES_RADIUS).
+    cap_edge = 0.1
+    rationals = [s.psifac for s in intr.sing]
+    isempty(rationals) || (cap_edge = min(cap_edge, minimum(rationals) - Equilibrium.RATIONAL_RES_RADIUS))
+    in_rational_window(x) = any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rationals)
     keep = Int[1]
     for i in 2:length(metric.xs)-1
         x = metric.xs[i]
-        if x >= 0.1 || (x - metric.xs[keep[end]]) >= 0.05 * x
+        if x >= cap_edge || in_rational_window(x) || (x - metric.xs[keep[end]]) >= 0.05 * x
             push!(keep, i)
         end
     end