Add rho^|m| near-axis regularization for VMEC harmonics

[krystophny]

Problem

VMEC full-grid Fourier amplitudes near the axis violate the analytic regularity condition c_m(rho) ~ rho^|m| (with rho = sqrt(s)). The existing axis healing (determine_nheal_for_axis, 30% tolerance) is a smoothness test: it is structurally blind to the smooth parity violation and lets sub-tolerance grid noise through. On the W7-X high-mirror equilibrium the measured rmnc/rho^m quotient of the dominant m=3 mode flips sign at the innermost surface and the m=4 quotient is 6.3x its asymptote, the textbook near-axis VMEC signature.

Downstream, SIMPLE's symplectic integrator converts that near-axis noise into secular radial transport: r_negative ~ 600 in a 0.01 s internal-Boozer trace, with a creeping non-prompt loss tail absent from RK and from external booz_xform chartmaps on the same equilibrium.

Fix

s_to_rho_power_law continues each harmonic to the axis with the analytic power law:

c(rho) = c_anchor * (rho/rho_anchor)^|m|      (rho < rho_axis_heal)

Below the anchor the surfaces are discarded; outside it the spline acts on c/rho^|m|, an even analytic function of rho, built only from reliable surfaces. m is clamped at 50. This is the booz_xform_to_boozer_chartmap.py continuation (which mirrors booz_xform's own rmnc/sqrt(s) extrapolation) applied at the VMEC harmonic level, so both the Boozer and the canonical transforms inherit a clean near-axis field from one shared evaluator chain (vmec_field -> splint_vmec_data). Neither transform needs its own change.

The regularity c_m ~ rho^|m| is standard (Garren-Boozer 1991; Landreman-Sengupta 2018/2019; DESC and booz_xform both enforce it).

Interface

New flags in new_vmec_stuff_mod:

• axis_healing_power_law (default .false.): opt in to the power-law continuation, bypassing determine_nheal_for_axis.
• rho_axis_heal (default 0.1): anchor radius. At ns=501 this replaces the inner 5 of 501 surfaces (inner 1% of flux); the field is unchanged beyond rho ~ 0.14.

The default keeps the legacy healing, so existing chartmap and field consumers are byte-for-byte unchanged. Consumers that need the regularization (SIMPLE near-axis symplectic tracing) set the flag from their own namelist.

Verification

W7-X high mirror, internal Boozer field, symplectic Euler1, 0.01 s, 1000 fixed-start alphas:

field	r_negative
internal Boozer, legacy healing	~606
internal Boozer, this PR (flag on)	2
booz_xform chartmap (reference)	21

The flag is opt-in because the continuation changes the near-axis field for every consumer. On the low-resolution CI fixture (ns=50) the c/rho^|m| division at the second surface amplifies the startup-layer noise instead of removing it. That distorts the near-axis geometry enough to diverge the chartmap Newton inversion (chartmap from_cyl failed ierr=1), which is what broke test_chartmap_matches_vmec and test_chartmap_vmec_mapping when the flag defaulted on. With the default .false. both chartmap tests pass unchanged; the regularization is exercised only by consumers that request it.

1 s benchmark runs (symplectic Boozer and canonical on the regularized field) are in progress in the orbit-benchmark repo.

Note for maintainers

The default .false. leaves the legacy near-axis field bit-reproducible, so existing golden records are unaffected. Enabling the flag shifts the near-axis harmonics by design; golden records generated under it differ from the legacy ones, and regenerating them is the intended follow-up for consumers that adopt it.

[krystophny]

Added: fix for the near-axis continuation + regression coverage

While adding the regression test that pins this regularization (it had no coverage), the test caught a correctness bug in s_to_rho_power_law.

The continuation is documented as c(rho) = c_anchor*(rho/rho_anchor)**|m|, but the code multiplied the reduced coefficient splcoe(0,1) = arr_in(i_anchor)/rho_anchor**|m| by (rho/rho_anchor)**|m|, leaving a spurious 1/rho_anchor**|m| factor. This is a discontinuity at the anchor surface (~10x for m=1, ~99x for m=2 at the default rho_axis_heal=0.1) that corrupts exactly the near-axis field the regularization cleans. Fixed by using splcoe(0,1)*rho**|m| (equal to arr_in(i_anchor)*(rho/rho_anchor)**|m|, continuous).

Verification

# before the fix (test asserts exact amp*rho**m reproduction)
m = 1   max |arr_out - amp*rho**m| = 0.61   (amp = 0.75)   -> ERROR STOP

# after the fix
test_axis_power_law_regularization ... Passed

The new test pins: axis vanishing for m>0, exact reproduction of an amp*rho**m harmonic on the full rho grid (m=1,2,4), anchor continuity, and the i_anchor index formula. The regularization stays opt-in (default off), so default behavior is unchanged.