Keep the default near-axis regularization radius on rho-native charts (#616)

## Cause
Since b4778e4 the defaults of `orbit_timestep_axis` (`axis_floor` 1e-8,
`near_axis_limit` 0.01) are read in the native radial unit of the chart.
For flux-native charts this is a regularization disc of Cartesian radius
1e-4 and a near-axis zone of radius 0.1, but on rho-native REFCOORDS
charts the same numbers shrank to 1e-8 and 0.01. An orbit crossing the
axis at Cartesian radius of about 1e-6 then meets the raw
polar-to-Cartesian pullback outside the regularized disc, and the RK45
integration exhausts its step budget (`ierr = 2`, libneo odeint max
steps).

## Fix
The defaults are now defined in normalized toroidal flux s (1e-8, 0.01)
and converted to the native coordinate (rho = sqrt(s)) when
`axis_radial_coordinate` reports a rho chart, so both chart types keep
the same Cartesian disc. Explicit `axis_floor` / `near_axis_limit`
arguments stay in native units. The radial-coordinate query moves before
the argument validation, nothing else changes.

## Test
`check_default_axis_disc` in `test_magfie_refcoords_splines` (NCSX rho
chart): six orbits starting at rho 1e-8 and 1e-6 with three pitches, one
step through the axis with default arguments, compared with the same
orbit in 40 short steps with explicit native limits (1e-4, 0.1). Without
the fix all six default-argument steps return `ierr = 2`; with it they
agree to about 1e-13.

Verified on macOS (Release, gfortran 16):
`test_magfie_refcoords_splines`, `test_axis_regular_rhs` and `ctest -R
'refcoords|axis|chartmap|orbit'` pass except four Python tests that need
numpy/netCDF4 missing on that host. The golden record has not been run.

Chris&AI

Author: krystophny

src/sub_alpha_lifetime_can.f90

@@ -901,8 +901,18 @@ subroutine orbit_timestep_axis(z, dtau, dtaumin, relerr, ierr, step_limit, &
                 end if
                 !
                 ierr = 0
+                call axis_radial_coordinate(radial_is_rho, ierr)
+                if (ierr /= 0) return
+                ! Defaults are given in normalized toroidal flux s and converted to
+                ! the native chart coordinate (rho = sqrt(s)), so that the
+                ! regularization disc has the same Cartesian radius for flux- and
+                ! rho-native charts. Explicit arguments stay in native units.
                 axis_floor_local = 1.d-8
                 snear_axis = 0.01d0
+                if (radial_is_rho) then
+                    axis_floor_local = sqrt(axis_floor_local)
+                    snear_axis = sqrt(snear_axis)
+                end if
                 if (present(axis_floor)) axis_floor_local = axis_floor
                 if (present(near_axis_limit)) snear_axis = near_axis_limit
                 if (.not. is_finite(axis_floor_local) .or. &
@@ -912,8 +922,6 @@ subroutine orbit_timestep_axis(z, dtau, dtaumin, relerr, ierr, step_limit, &
                     ierr = 2
                     return
                 end if
-                call axis_radial_coordinate(radial_is_rho, ierr)
-                if (ierr /= 0) return
              near_radius_squared = axis_native_radius_squared(snear_axis, radial_is_rho)
                 z_initial = z
                 y(1) = z(1)

test/tests/magfie/test_magfie_refcoords_splines.f90

@@ -76,6 +76,7 @@ program test_magfie_refcoords_splines
         error stop 'inactive reference field changed the RK axis chart'
     call init_magfie(REFCOORDS)
     call check_active_rho_rhs(n_failed)
+    call check_default_axis_disc(n_failed)
 
     phi_period = twopi/real(nper, dp)
 
@@ -231,6 +232,55 @@ subroutine check_active_rho_rhs(n_failed)
         end if
     end subroutine check_active_rho_rhs
 
+    subroutine check_default_axis_disc(n_failed)
+        !> The default regularization disc and near-axis zone are defined in flux
+        !> units (1e-8, 0.01) and must keep their Cartesian size (1e-4, 0.1) on a
+        !> rho-native chart. A step that starts within 1e-6 of the axis crosses it
+        !> at Cartesian radius ~1e-6 and used to exhaust the RK step budget with
+        !> the unconverted defaults (ierr = 2). The reference integrates the same
+        !> orbit in many short steps with the equivalent explicit native limits.
+        integer, intent(inout) :: n_failed
+        type(tracer_t) :: orbit
+        real(dp), parameter :: rho0(2) = [1.0e-8_dp, 1.0e-6_dp]
+        real(dp), parameter :: pitch(3) = [0.3_dp, -0.6_dp, 0.05_dp]
+        real(dp), parameter :: tol = 1.0e-9_dp
+        real(dp) :: z(5), zref(5), dist, dtau
+        integer :: irho, ipitch, istep, status, status_ref
+        integer, parameter :: n_sub = 40
+
+        call init_params(orbit, 2, 4, 3.5e6_dp, 256, 1, 1.0e-10_dp)
+        dtau = 0.2_dp*orbit%dtau
+        do irho = 1, size(rho0)
+            do ipitch = 1, size(pitch)
+                z = [rho0(irho), 0.4_dp, 0.1_dp, 1.0_dp, pitch(ipitch)]
+                zref = z
+                call orbit_timestep_axis(z, dtau, dtau, 1.0e-10_dp, status)
+                status_ref = 0
+                do istep = 1, n_sub
+                    if (status_ref == 0) &
+                        call orbit_timestep_axis(zref, dtau/n_sub, dtau/n_sub, &
+                                                 1.0e-10_dp, status_ref, &
+                                                 axis_floor=1.0e-4_dp, &
+                                                 near_axis_limit=0.1_dp)
+                end do
+                if (status /= 0 .or. status_ref /= 0 .or. &
+                    .not. all(is_finite([z, zref]))) then
+                    print *, 'FAIL: default axis disc orbit failed', irho, ipitch, &
+                        status, status_ref
+                    n_failed = n_failed + 1
+                    cycle
+                end if
+                dist = hypot(z(1)*cos(z(2)) - zref(1)*cos(zref(2)), &
+                             z(1)*sin(z(2)) - zref(1)*sin(zref(2)))
+                if (dist > tol .or. maxval(abs(z(3:5) - zref(3:5))) > tol) then
+                    print *, 'FAIL: default axis disc differs from explicit disc', &
+                        irho, ipitch, dist
+                    n_failed = n_failed + 1
+                end if
+            end do
+        end do
+    end subroutine check_default_axis_disc
+
     subroutine check_point(i, phi_period, f, n_failed)
         integer, intent(in) :: i
         real(dp), intent(in) :: phi_period