diff --git a/POTATO/CMakeLists.txt b/POTATO/CMakeLists.txt
index 25e3b4f4..a6677b30 100644
--- a/POTATO/CMakeLists.txt
+++ b/POTATO/CMakeLists.txt
@@ -5,7 +5,7 @@ enable_language(Fortran)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/../cmake")
 
 set(CMAKE_Fortran_MODULE_DIRECTORY ${PROJECT_BINARY_DIR})
-add_compile_options(-cpp -g -fbacktrace -Wfatal-errors -std=f2018)
+add_compile_options(-g -fbacktrace)
 add_compile_options(-O3 -march=native -mtune=native)
 
 include(Util)
@@ -21,6 +21,33 @@ FetchContent_Declare(
 )
 FetchContent_MakeAvailable(fortplot)
 
+# fortnum provides the variable-order Adams ODE integrator with event/root
+# detection used by find_bounce to locate the bounce return (replacing the
+# fixed-step Poincare-cut crossing search). Same pin as the NEO-RT migration.
+if(NOT TARGET fortnum)
+    FetchContent_Declare(
+        fortnum
+        GIT_REPOSITORY https://github.com/lazy-fortran/fortnum.git
+        GIT_TAG 92de6e949a772cfffc73bb5295fe5e2b056b9c18
+    )
+    FetchContent_MakeAvailable(fortnum)
+endif()
+# fortnum relies on a GNU implicit-typing extension; pin it to the GNU standard.
+if(TARGET fortnum)
+    target_compile_options(fortnum PRIVATE -std=gnu)
+endif()
+
+# libneo supplies the field, magfie and spline code POTATO used to carry as its
+# own copies. Link the upstream library so fixes (e.g. the stretch_coords
+# convex-wall load) propagate instead of drifting in a fork.
+find_or_fetch(libneo)
+
+# Hold POTATO's own sources to f2018, stop at the first diagnostic, and run the C
+# preprocessor. Applied after the dependency fetches above so the upstream
+# libraries (libneo, fortplot, fortnum) build under their own standards instead
+# of inheriting these; -cpp in particular makes CMake mis-parse libneo's .f90.
+add_compile_options(-cpp -Wfatal-errors -std=f2018)
+
 add_subdirectory(SRC)
 
 add_executable (potato.x
diff --git a/POTATO/SRC/CMakeLists.txt b/POTATO/SRC/CMakeLists.txt
index 36673445..cab99073 100644
--- a/POTATO/SRC/CMakeLists.txt
+++ b/POTATO/SRC/CMakeLists.txt
@@ -1,25 +1,11 @@
+# The field, magfie, spline and kinds/constants sources live in libneo and are
+# linked below (LIBNEO::magfie pulls in libneo's core). Only POTATO-specific
+# sources are compiled here. binsrc/plag_coeff stay local: POTATO calls them as
+# bare external subroutines, whereas libneo wraps them in *_sub modules.
 add_library(potato_base
-  libneo_kinds.f90
   logging_mod.f90
   potato_input_mod.f90
-  math_constants.f90
   alpha_lifetime_mod.f90
-  chamb_divB0.f90
-  magfie_cyl.f90
-  period_mod.f90
-  input_files.f90
-  field_c_mod.f90
-  field_mod.f90
-  inthecore_mod.f90
-  field_eq_mod.f90
-  field_divB0.f90
-  bdivfree_mod.f90
-  amn_mod.f90
-  theta_rz_mod.f90
-  extract_fluxcoord_mod.f90
-  bdivfree.f90
-  spl_three_to_five.f90
-  spline5_RZ.f90
   velo.f90
   odeint_allroutines.f
   plag_coeff.f90
@@ -38,10 +24,13 @@ add_library(potato_base
 )
 find_or_fetch(vode)
 # vode (dvode_f90_m) uses labeled DO termination, deleted in f2018; build it with
-# the permissive legacy standard rather than the project-wide f2018.
+# the permissive legacy standard rather than POTATO's f2018.
 target_compile_options(vode PRIVATE -std=legacy)
 target_link_libraries(potato_base PUBLIC
 vode
+fortnum
+LIBNEO::magfie
+mc_efit
 BLAS::BLAS
 LAPACK::LAPACK
 HDF5::HDF5
diff --git a/POTATO/SRC/amn_mod.f90 b/POTATO/SRC/amn_mod.f90
deleted file mode 120000
index ebf63216..00000000
--- a/POTATO/SRC/amn_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/amn_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/bdivfree.f90 b/POTATO/SRC/bdivfree.f90
deleted file mode 120000
index bdc437df..00000000
--- a/POTATO/SRC/bdivfree.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/bdivfree.f90
\ No newline at end of file
diff --git a/POTATO/SRC/bdivfree_mod.f90 b/POTATO/SRC/bdivfree_mod.f90
deleted file mode 120000
index 207ce50f..00000000
--- a/POTATO/SRC/bdivfree_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/bdivfree_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/box_counting.f90 b/POTATO/SRC/box_counting.f90
index b32bfdde..4d2385fb 100644
--- a/POTATO/SRC/box_counting.f90
+++ b/POTATO/SRC/box_counting.f90
@@ -24,7 +24,7 @@ end subroutine timestep_vode
 
 subroutine time_in_box(z, cnt, sbox, taub, tau)
   ! Returns time spent in boxes
-  use dvode_f90_m, only: vode_opts, set_normal_opts, dvode_f90, get_stats
+  use dvode_f90_m, only: vode_opts, set_opts, dvode_f90, get_stats
   use field_sub, only: psif
   use field_eq_mod, only: psi_axis, psi_sep
   use orbit_dim_mod, only: neqm
@@ -37,23 +37,18 @@ subroutine time_in_box(z, cnt, sbox, taub, tau)
   real(8), intent(in)    :: taub       ! Bounce time
   real(8), intent(out)   :: tau(cnt)   ! Time in each box
 
-  real(8) :: delphi
-  real(8) :: sprev, snext ! Previous and next flux radius box
+  real(8) :: smid
   real(8) :: y(neqm)
 
-  integer(4) :: nmax  ! Maximum loop iterations
-
-  integer(4) :: k
+  integer(4) :: k, nsample
   real(8) :: ti
 
-  real(8) :: atol(neqm), rtol, tout, rstats(22)
-  integer(4) :: itask, istate, istats(31), numevents
+  real(8) :: atol(neqm), rtol, tout
+  integer(4) :: itask, istate, method_flag
   type (vode_opts) :: options
   real(8) :: bmod, phi_elec, s
-  real(8) :: sold, told, yold(neqm)
-  integer(4) :: sind, sind0 ! s index
-
-  integer(4) :: jroots(2)
+  real(8) :: sold, told
+  integer(4) :: sind ! s index
 
   external timestep_vode
 
@@ -63,83 +58,58 @@ subroutine time_in_box(z, cnt, sbox, taub, tau)
   atol = 1d-13
   itask = 1
   istate = 1
-  numevents = 2
-  options = set_normal_opts(abserr_vector=atol, relerr=rtol, nevents=numevents)
+  method_flag = 10
+  options = set_opts(method_flag=method_flag, abserr_vector=atol, &
+    relerr=rtol, mxstep=2000)
 
-  nmax = 3*size(sbox)
+  nsample = min(256, max(64, size(sbox)))
   tau = 0d0
   ti = 0d0
 
   call get_bmod_and_Phi(z(1:3), bmod, phi_elec)
 
   s = abs((psif-psi_axis)/(psi_sep-psi_axis))
-  sind = size(sbox)+1
-  sprev = -1e5
-  snext = 1e5
-  do k = 1,size(sbox)
-    if(sbox(k)>s) then
-        sind = k
-        snext = sbox(k)
-        if(k>1) sprev = sbox(k-1)
-        exit
-    end if
-  enddo
-  sind0 = sind
-
+  sold = s
   told = 0d0
-  do k = 1,nmax
-    yold = y
-    sold = s
-    tout = taub
+  do k = 1,nsample
+    tout = taub*dble(k)/dble(nsample)
     call dvode_f90( &
-        timestep_vode, neqm, y, ti, tout, itask, istate, options, &
-        g_fcn = sroots )
-    if (istate == 2) exit
-    if (istate == 3) then
-        tau(sind) = tau(sind) + ti-told
-        told = ti
-        call get_stats(rstats, istats, numevents, jroots)
-        if (jroots(2).ne.0) then
-          sind = sind + 1 ! moving outwards
-          sprev = snext
-          if (sind == size(sbox)+1) then
-              snext = 1e5
-          else
-              snext = sbox(sind)
-          end if
-        end if
-        if (jroots(1).ne.0) then
-          sind = sind - 1 ! moving inwards
-          snext = sprev
-          if (sind == 1) then
-              sprev = -1e5
-          else
-              sprev = sbox(sind-1)
-          end if
-        end if
-    elseif (istate == -1) then
-      print *, 'Error in VODE: k =', k
-      istate = 1
-      exit
+        timestep_vode, neqm, y, ti, tout, itask, istate, options)
+    if (istate == -1) then
+      print *, 'time_in_box: VODE exceeded mxstep at k =', k, &
+        ' ti =', ti, ' taub =', taub
+      tau = 0d0
+      return
+    endif
+    if (istate /= 2) then
+      print *, 'time_in_box: VODE failed with istate =', istate, &
+        ' k =', k, ' ti =', ti, ' taub =', taub
+      tau = 0d0
+      return
     end if
-  end do
-
-  tau(sind) = tau(sind) + taub-told
 
-contains
-
-  subroutine sroots(neqext, t, yext, ng, gout)
-    ! For finding roots between boxes
+    call get_bmod_and_Phi(y(1:3), bmod, phi_elec)
 
-    integer, intent(in) :: neqext, ng
-    real(8), intent(in) :: t, yext(neqext)
-    real(8), intent(out) :: gout(ng)
+    s = abs((psif-psi_axis)/(psi_sep-psi_axis))
+    smid = 0.5d0*(sold+s)
+    sind = radial_box_index(smid)
+    tau(sind) = tau(sind) + ti-told
+    sold = s
+    told = ti
+  end do
 
-    call get_bmod_and_Phi(yext(1:3), bmod, phi_elec)
+contains
 
-    s = abs((psif-psi_axis)/(psi_sep-psi_axis))
+  integer function radial_box_index(sval)
+    real(8), intent(in) :: sval
+    integer :: ibox
 
-    gout(1) = s - sprev
-    gout(2) = s - snext
-  end subroutine sroots
+    radial_box_index = size(sbox)
+    do ibox = 1,size(sbox)
+      if (sval <= sbox(ibox)) then
+        radial_box_index = ibox
+        exit
+      endif
+    enddo
+  end function radial_box_index
 end subroutine time_in_box
diff --git a/POTATO/SRC/chamb_divB0.f90 b/POTATO/SRC/chamb_divB0.f90
deleted file mode 120000
index 6bd4433a..00000000
--- a/POTATO/SRC/chamb_divB0.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/MC/chamb_divB0.f90
\ No newline at end of file
diff --git a/POTATO/SRC/extract_fluxcoord_mod.f90 b/POTATO/SRC/extract_fluxcoord_mod.f90
deleted file mode 120000
index 460a2739..00000000
--- a/POTATO/SRC/extract_fluxcoord_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/extract_fluxcoord_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/field_c_mod.f90 b/POTATO/SRC/field_c_mod.f90
deleted file mode 120000
index 81a304b0..00000000
--- a/POTATO/SRC/field_c_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/field_c_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/field_divB0.f90 b/POTATO/SRC/field_divB0.f90
deleted file mode 120000
index ba979e1a..00000000
--- a/POTATO/SRC/field_divB0.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/field_divB0.f90
\ No newline at end of file
diff --git a/POTATO/SRC/field_eq_mod.f90 b/POTATO/SRC/field_eq_mod.f90
deleted file mode 120000
index b4842173..00000000
--- a/POTATO/SRC/field_eq_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/field_eq_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/field_mod.f90 b/POTATO/SRC/field_mod.f90
deleted file mode 120000
index 4c8e4f4b..00000000
--- a/POTATO/SRC/field_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/field_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/find_all_roots.f90 b/POTATO/SRC/find_all_roots.f90
index ed6dfa08..09dec90a 100644
--- a/POTATO/SRC/find_all_roots.f90
+++ b/POTATO/SRC/find_all_roots.f90
@@ -5,7 +5,8 @@
 !
   module find_all_roots_mod
     logical :: customgrid=.false.
-    integer :: nroots, nsearch_min=100, ncustom, niter=100
+    logical :: fail_fast=.false.
+    integer :: nroots, nsearch_min=100, ncustom, niter=100, max_roots_abort=0
     double precision :: relerr_allroots=1.d-12
     double precision, dimension(:), allocatable :: xcustom,roots
   end module find_all_roots_mod
@@ -25,7 +26,7 @@ subroutine find_all_roots(fun,x1in,x2in,ierr)
 !
   use find_all_roots_mod, only : nsearch_min,niter,relerr_allroots,  &
                                  customgrid,ncustom,xcustom,         &
-                                 nroots,roots
+                                 nroots,roots,fail_fast,max_roots_abort
 !
   implicit none
 !
@@ -33,9 +34,11 @@ subroutine find_all_roots(fun,x1in,x2in,ierr)
   double precision :: x,x1in,x2in,f,df,hx,dx,errdist,xb,xe,dfb,dfe,xxtr,fxtr
   double precision :: x1,x2
   double precision, dimension(:), allocatable :: xarr,farr,dfarr,dummy1d
+  logical :: abort_search
   external :: fun
 !
   ierr=0
+  abort_search=.false.
   x1=x1in
   x2=x2in
 !
@@ -137,6 +140,7 @@ subroutine find_all_roots(fun,x1in,x2in,ierr)
         x=(farr(i)*xarr(i-1)-farr(i-1)*xarr(i))/(farr(i)-farr(i-1))
 !
         call addroot
+        if(abort_search) exit
 !
       endif
     else
@@ -166,6 +170,7 @@ subroutine find_all_roots(fun,x1in,x2in,ierr)
         x=(fxtr*xarr(i-1)-farr(i-1)*xxtr)/(fxtr-farr(i-1))
 !
         call addroot
+        if(abort_search) exit
 !
       endif
       if(fxtr*farr(i).le.0.d0) then
@@ -173,9 +178,11 @@ subroutine find_all_roots(fun,x1in,x2in,ierr)
         x=(farr(i)*xxtr-fxtr*xarr(i))/(farr(i)-fxtr)
 !
         call addroot
+        if(abort_search) exit
 !
       endif
     endif
+    if(abort_search) exit
   enddo
 !
 !------------
@@ -189,6 +196,12 @@ subroutine addroot
 ! Root adjustment by Newton method
 !
   double precision :: xx
+!
+  if(max_roots_abort.gt.0 .and. nroots.ge.max_roots_abort) then
+    ierr=2
+    abort_search=.true.
+    return
+  endif
 !
 ! extend storage arragy for a new root:
   nroots=nroots+1
@@ -222,9 +235,22 @@ subroutine addroot
   enddo
 !
   if(iter.gt.niter) then
+    ierr=1
+    if(fail_fast) then
+      abort_search=.true.
+      nroots=nroots-1
+      if(nroots.eq.0) then
+        deallocate(roots)
+      else
+        dummy1d=roots(1:nroots)
+        deallocate(roots)
+        allocate(roots(nroots))
+        roots=dummy1d
+      endif
+      return
+    endif
     print *,'No convergence of Newton in find_all_roots: error = ',sngl(abs(dx)), &
             ' tolerance = ',sngl(errdist)
-    ierr=1
   endif
 !
   roots(nroots)=x
diff --git a/POTATO/SRC/input_files.f90 b/POTATO/SRC/input_files.f90
deleted file mode 120000
index 1fd65563..00000000
--- a/POTATO/SRC/input_files.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/input_files.f90
\ No newline at end of file
diff --git a/POTATO/SRC/inthecore_mod.f90 b/POTATO/SRC/inthecore_mod.f90
deleted file mode 120000
index 6a056a5d..00000000
--- a/POTATO/SRC/inthecore_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/inthecore_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/libneo_kinds.f90 b/POTATO/SRC/libneo_kinds.f90
deleted file mode 120000
index b880394d..00000000
--- a/POTATO/SRC/libneo_kinds.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/libneo_kinds.f90
\ No newline at end of file
diff --git a/POTATO/SRC/magfie_cyl.f90 b/POTATO/SRC/magfie_cyl.f90
deleted file mode 120000
index a248f9fe..00000000
--- a/POTATO/SRC/magfie_cyl.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/magfie_cyl.f90
\ No newline at end of file
diff --git a/POTATO/SRC/math_constants.f90 b/POTATO/SRC/math_constants.f90
deleted file mode 120000
index b87bd015..00000000
--- a/POTATO/SRC/math_constants.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/math_constants.f90
\ No newline at end of file
diff --git a/POTATO/SRC/period_mod.f90 b/POTATO/SRC/period_mod.f90
deleted file mode 120000
index 266c77e8..00000000
--- a/POTATO/SRC/period_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/period_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/potato_input_mod.f90 b/POTATO/SRC/potato_input_mod.f90
index a72c8e5c..19718e10 100644
--- a/POTATO/SRC/potato_input_mod.f90
+++ b/POTATO/SRC/potato_input_mod.f90
@@ -43,17 +43,37 @@ module potato_input_mod
     integer :: nlagr_sampling = 3
     double precision :: eps_sampling = 1d-2
     integer :: itermax_sampling = 5
+    logical :: clip_resonance_classes = .true.
 
     ! Orbit plotting
     double precision :: toten_plot = 1d0
     double precision :: perpinv_plot = 4.5d-5
 
+    ! Single-orbit trace (itest_type=4): start point on the poloidal plane and
+    ! pitch cosine of one guiding-center orbit, traced over one bounce period.
+    double precision :: orbit_Rstart = 210d0   ! start major radius [cm]
+    double precision :: orbit_Zstart = 0d0     ! start height [cm]
+    double precision :: orbit_lambda = 0.4d0   ! pitch cosine v_par/v at start
+
+    ! Frequency radial scan (itest_type=5): trace one bounce at each of freq_n
+    ! midplane start radii from freq_Rmin to freq_Rmax (pitch orbit_lambda),
+    ! and write rho_pol, omega_b, omega_phi to freq_scan.dat.
+    double precision :: freq_Rmin = 178d0      ! inner start radius [cm]
+    double precision :: freq_Rmax = 221d0      ! outer start radius [cm], into SOL
+    integer          :: freq_n = 40            ! number of surfaces
+
     ! Canonical frequency plot
     double precision :: enkin_over_temp = 1d0
 
     ! Profile input file
     character(len=256) :: profile_file = 'profile_poly.in'
 
+    ! Edge extension: when .true. the orbit integrator accepts guiding-center
+    ! points outside the separatrix (scrape-off layer), so bananas whose tips
+    ! cross the last closed flux surface can be traced in the extended
+    ! equilibrium field. Sets field_eq_mod::allow_sol. Default .false.
+    logical :: edge_extension = .false.
+
     namelist /potato_nml/ &
         itest_type, E_alpha, A_alpha, Z_alpha, &
         rho_pol, rho_pol_max, scalfac_energy, scalfac_efield, &
@@ -61,9 +81,11 @@ module potato_input_mod
         m_min, m_max, n_tor, &
         nenerg, thermen_max, nbox, &
         adaptive_jperp, npoi_init, nlagr_sampling, eps_sampling, &
-        itermax_sampling, &
+        itermax_sampling, clip_resonance_classes, &
         toten_plot, perpinv_plot, enkin_over_temp, &
-        profile_file
+        profile_file, edge_extension, &
+        orbit_Rstart, orbit_Zstart, orbit_lambda, &
+        freq_Rmin, freq_Rmax, freq_n
 
 contains
 
@@ -122,10 +144,18 @@ subroutine print_potato_input(iunit)
         write(iunit, '(A,I0)') '  nlagr_sampling   = ', nlagr_sampling
         write(iunit, '(A,ES12.5)') '  eps_sampling     = ', eps_sampling
         write(iunit, '(A,I0)') '  itermax_sampling = ', itermax_sampling
+        write(iunit, '(A,L1)') '  clip_resonance_classes = ', clip_resonance_classes
         write(iunit, '(A,ES12.5)') '  toten_plot       = ', toten_plot
         write(iunit, '(A,ES12.5)') '  perpinv_plot     = ', perpinv_plot
         write(iunit, '(A,ES12.5)') '  enkin_over_temp  = ', enkin_over_temp
         write(iunit, '(A,A)') '  profile_file     = ', trim(profile_file)
+        write(iunit, '(A,L1)') '  edge_extension   = ', edge_extension
+        write(iunit, '(A,ES12.5)') '  orbit_Rstart     = ', orbit_Rstart
+        write(iunit, '(A,ES12.5)') '  orbit_Zstart     = ', orbit_Zstart
+        write(iunit, '(A,ES12.5)') '  orbit_lambda     = ', orbit_lambda
+        write(iunit, '(A,ES12.5)') '  freq_Rmin        = ', freq_Rmin
+        write(iunit, '(A,ES12.5)') '  freq_Rmax        = ', freq_Rmax
+        write(iunit, '(A,I0)') '  freq_n           = ', freq_n
         write(iunit, '(A)') '================================'
     end subroutine print_potato_input
 
diff --git a/POTATO/SRC/resonant_int.f90 b/POTATO/SRC/resonant_int.f90
index a9ba3dc1..83bc63db 100644
--- a/POTATO/SRC/resonant_int.f90
+++ b/POTATO/SRC/resonant_int.f90
@@ -134,44 +134,57 @@ subroutine integrate_class_resonances
 !
 ! Computes sum over resonances $x=x^{res}_{(\bm,k)}$ in Eq.(104) for a given class $k$
 !
-  use find_all_roots_mod, only : customgrid,ncustom,xcustom,nroots,roots
-  use get_matrix_mod,     only : relmargin,iclass
+  use find_all_roots_mod, only : customgrid,ncustom,xcustom,nroots,roots, &
+                                 relerr_allroots,fail_fast,max_roots_abort
+  use get_matrix_mod,     only : iclass
   use form_classes_doublecount_mod, only : ifuntype,R_class_beg,R_class_end,sigma_class
   use resint_mod,         only : nmodes,marr,narr,twopim2,rm3,delint_mode,respoints_jp
-  use orbit_dim_mod,      only : neqm
-  use global_invariants,  only : toten,perpinv,cE_ref,Phi_eff
-  use sample_matrix_mod,  only : npoi,xarr
+  use orbit_dim_mod,      only : neqm,fb_max_abs_delphi,fb_max_tau
+  use global_invariants,  only : toten,perpinv,cE_ref,Phi_eff,dtau
+! xbeg,xend are the search bounds set (and clipped to the |Delta_phi_b| cap) by
+! the preceding sample_class_doublecount, so the root search runs over exactly
+! the sampled interval.
+  use sample_matrix_mod,  only : npoi,xarr,xbeg,xend,eps
   use logging_mod,        only : tee_message
 !
   implicit none
 !
   double precision, parameter :: pi=3.14159265358979d0, twopi=2.d0*pi, &
                                  pi32_over4m=-0.25d0*pi**1.5d0
+  integer, parameter :: max_roots_per_class_mode=32
 !
-  integer          :: mode,iroot,ierr
-  double precision :: relmargin_loc,widthclass,xbeg,xend
+  integer          :: mode,iroot,ierr,max_roots_abort_save
   double precision :: rescond,dresconddx,psiast,dpsiastdx,taub,delphi
   double precision :: one_res,sigma,delta_R,Rst,xi,dxi_dx,dpsiast_dRst,absHn2
   double precision :: toten_loc,perpinv_loc,fmaxw,A1ast,A2ast
   double precision :: dens,temp,ddens,dtemp,phi_elec,dPhi_dpsi
+  double precision :: relerr_allroots_save,fb_max_abs_delphi_save,fb_max_tau_save
   double precision, dimension(neqm) :: z
+  character(len=256) :: msg
+  logical :: fail_fast_save
 !
   toten_loc=toten
   perpinv_loc=perpinv
 !
   sigma=sigma_class(iclass)
   delta_R=R_class_end(iclass)-R_class_beg(iclass)
-!old=>  relmargin_loc=1.d-8
-!old=>  widthclass=1.d0
-  relmargin_loc=relmargin                                      !<=new
-  widthclass=abs(R_class_end(iclass)/R_class_beg(iclass)-1.d0) !<=new
-!
-  call classbounds(ifuntype(iclass),relmargin_loc,widthclass,xbeg,xend)
 !
-  customgrid=.true.
-  ncustom=npoi
-  allocate(xcustom(ncustom))
-  xcustom=xarr
+  ncustom=count(xarr.ge.xbeg .and. xarr.le.xend)
+  customgrid=ncustom.ge.2
+  if(customgrid) then
+    allocate(xcustom(ncustom))
+    xcustom=pack(xarr,xarr.ge.xbeg .and. xarr.le.xend)
+  endif
+  relerr_allroots_save=relerr_allroots
+  relerr_allroots=max(1.d-8,1.d-2*eps)
+  fail_fast_save=fail_fast
+  fail_fast=.true.
+  max_roots_abort_save=max_roots_abort
+  max_roots_abort=max_roots_per_class_mode+1
+  fb_max_abs_delphi_save=fb_max_abs_delphi
+  fb_max_tau_save=fb_max_tau
+  fb_max_abs_delphi=1.5d0*twopi*dble(maxval(abs(marr)))/dble(abs(narr(1)))
+  fb_max_tau=200.d0*dtau
 !
 !
   do mode=1,nmodes
@@ -180,13 +193,25 @@ subroutine integrate_class_resonances
     delint_mode(mode)=0.d0
 !
     call find_all_roots(get_rescond,xbeg,xend,ierr)
+!
+    if(ierr.eq.2 .or. nroots.gt.max_roots_per_class_mode) then
+      write(msg,'(A,ES12.4,A,ES12.4,A,I0,A,I0,A,I0,A,I0)') &
+        'skipped oscillatory roots: toten=',toten,' perpinv=',perpinv, &
+        ' iclass=',iclass,' mode=',mode,' m=',marr(mode), &
+        ' nroots=',nroots
+      call tee_message(trim(msg))
+      respoints_jp(mode,iclass)%nrespoi=0
+      respoints_jp(mode,iclass)%toten_res=toten
+      respoints_jp(mode,iclass)%perpinv_res=perpinv
+      cycle
+    endif
 !
     if(ierr.ne.0) then
-      call tee_message( &
-        'integrate_class_resonances: error in find_all_roots')
-      customgrid=.false.
-      deallocate(xcustom)
-      return
+      call tee_message('integrate_class_resonances: skipped unpolished roots')
+      respoints_jp(mode,iclass)%nrespoi=0
+      respoints_jp(mode,iclass)%toten_res=toten
+      respoints_jp(mode,iclass)%perpinv_res=perpinv
+      cycle
     endif
 !
     respoints_jp(mode,iclass)%nrespoi=nroots
@@ -254,7 +279,12 @@ subroutine integrate_class_resonances
   enddo
 !
   customgrid=.false.
-  deallocate(xcustom)
+  relerr_allroots=relerr_allroots_save
+  fail_fast=fail_fast_save
+  max_roots_abort=max_roots_abort_save
+  fb_max_abs_delphi=fb_max_abs_delphi_save
+  fb_max_tau=fb_max_tau_save
+  if(allocated(xcustom)) deallocate(xcustom)
 !
 !------------
   contains
@@ -262,25 +292,33 @@ subroutine integrate_class_resonances
 !
   subroutine get_rescond(x,rescond,dresconddx)
 !
-! Computes resonance condition $F(x)=\Delta\varphi_b+2\pi m_2/m_3$
-! and its derivative $F^\prime(x)$ for $F(x)=0$ root finding.
-! Computes as by-products normalized toroidal momentum $\psi^\ast$,
-! bounce time $\tau_b$ and toroidal displacement $\Delta\varphi_b$.
+! Resonance condition in the canonical FREQUENCIES, $F(x)=m\,\Omega_b+n\,\Omega_\varphi$,
+! and its derivative, for $F(x)=0$ root finding.  This is the bounce-drift
+! resonance $m\,\Omega_b+n\,\Omega_\varphi=0$ -- the same roots as
+! $\Delta\varphi_b+2\pi m/n=0$ (multiply by $\tau_b/2\pi$), but $\Omega_b=2\pi/\tau_b$
+! and $\Omega_\varphi=\Delta\varphi_b/\tau_b$ stay finite at the separatrix where
+! $\tau_b\to\infty$, so the interpolated grid (built by get_matrix_doublecount,
+! which now stores the frequencies in rows 2,3) converges instead of diverging.
+! By-products $\psi^\ast$, $\tau_b$, $\Delta\varphi_b$ are recovered for the weight.
 !
   use sample_matrix_mod, only : n1
 !
   implicit none
 !
-  double precision :: x,rescond,dresconddx
+  double precision, parameter :: twopi=6.28318530717958648d0
+  double precision :: x,rescond,dresconddx,omega_b,omega_phi,rm2
   double precision, dimension(n1) :: vec,dvec
 !
   call interpolate_class_doublecount(x,vec,dvec)
 !
-  psiast=vec(1)               ! $\psi^\ast$
-  taub=vec(2)                 ! $\tau_b$
-  delphi=vec(3)               ! $\Delta\varphi_b$
-  rescond=delphi+twopim2/rm3
-  dresconddx=dvec(3)
+  psiast=vec(1)                          ! $\psi^\ast$
+  omega_b=vec(2)                         ! $\Omega_b=2\pi/\tau_b$
+  omega_phi=vec(3)                       ! $\Omega_\varphi=\Delta\varphi_b/\tau_b$
+  taub=twopi/omega_b                     ! recover $\tau_b$ (finite at the root)
+  delphi=omega_phi*taub                  ! recover $\Delta\varphi_b$
+  rm2=twopim2/twopi                      ! poloidal harmonic $m=$ twopim2$/2\pi$
+  rescond=rm2*omega_b+rm3*omega_phi      ! $m\,\Omega_b+n\,\Omega_\varphi$
+  dresconddx=rm2*dvec(2)+rm3*dvec(3)
 !
   end subroutine get_rescond
 !
@@ -421,6 +459,7 @@ subroutine resonant_torque
   double precision :: omdens,trapez_fac,perpinv_max
   double precision :: torque_int,torque_int_loc
   double precision :: xjperp,xenerg,totxint,step_energ
+  double precision :: w_energ,dw_energ
   double precision :: time_beg,time_end
   double precision :: dens, temp, ddens, dtemp
   character(len=256) :: msg
@@ -483,8 +522,10 @@ subroutine resonant_torque
 !
   torquebox=0.d0
 !
-  step_energ=toten_range/dble(nenerg) !integration step over total energy
-!step_energ=0.22521463755624047d0
+! Energy grid uniform in velocity w = sqrt(toten - toten_min) ~ v, so the scan
+! starts at low energy while still covering the same maximum v/vth range as
+! NEO-RT. The energy integral keeps its value through d(toten) = 2 w dw.
+  dw_energ=sqrt(toten_range)/dble(nenerg)
 !
   if(adaptive_jperp) then
     open(1901,file='subint_ofH0int_104_vsJperp_fromresp.dat')
@@ -493,12 +534,10 @@ subroutine resonant_torque
     open(1902,file='subint_ofH0int_104_vsJperp_equi.dat')
   endif
 !
-!  do ienerg=1,nenerg
   do ienerg=2,nenerg
-!  do ienerg= 20,20 !10,10 !20,20 !<=fix energy for debugging
-    xenerg=(dble(ienerg)-0.5d0)/dble(nenerg)
-    toten=toten_min+toten_range*xenerg
-!toten =   -3.1211921097605737d0
+    w_energ=(dble(ienerg)-0.5d0)*dw_energ
+    toten=toten_min+w_energ**2
+    step_energ=2.d0*w_energ*dw_energ
     write(msg, '(A,ES22.14)') 'toten = ', toten
     call tee_message(trim(msg))
 !
diff --git a/POTATO/SRC/rng.f90 b/POTATO/SRC/rng.f90
deleted file mode 120000
index 38273697..00000000
--- a/POTATO/SRC/rng.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/MC/rng.f90
\ No newline at end of file
diff --git a/POTATO/SRC/sample_matrix.f90 b/POTATO/SRC/sample_matrix.f90
index 40c1a33d..fa5e6a71 100644
--- a/POTATO/SRC/sample_matrix.f90
+++ b/POTATO/SRC/sample_matrix.f90
@@ -38,9 +38,12 @@ SUBROUTINE sample_matrix(get_matrix,ierr)
   h=(xend-xbeg)/npoilag/(1.d0+symm_break)
   hh=symm_break*h/npoilag
 !
-  if(allocated(amat)) then
-    DEALLOCATE(amat,xarr,amat_arr)
-  endif
+! amat is the shared sample_matrix_mod array; callers (e.g. the resonance class
+! preclip) may have allocated it alone, so deallocate each array on its own
+! status rather than gating all three on amat being allocated.
+  if(allocated(amat))     DEALLOCATE(amat)
+  if(allocated(xarr))     DEALLOCATE(xarr)
+  if(allocated(amat_arr)) DEALLOCATE(amat_arr)
 !
   ALLOCATE(amat(n1,n2),xarr(npoi),amat_arr(n1,n2,npoi))
 !
diff --git a/POTATO/SRC/spl_three_to_five.f90 b/POTATO/SRC/spl_three_to_five.f90
deleted file mode 120000
index 59516009..00000000
--- a/POTATO/SRC/spl_three_to_five.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/spl_three_to_five.f90
\ No newline at end of file
diff --git a/POTATO/SRC/spline5_RZ.f90 b/POTATO/SRC/spline5_RZ.f90
deleted file mode 120000
index 4496bd49..00000000
--- a/POTATO/SRC/spline5_RZ.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/spline5_RZ.f90
\ No newline at end of file
diff --git a/POTATO/SRC/sub_potato.f90 b/POTATO/SRC/sub_potato.f90
index b4d77f56..6aa93157 100644
--- a/POTATO/SRC/sub_potato.f90
+++ b/POTATO/SRC/sub_potato.f90
@@ -81,6 +81,18 @@ module orbit_dim_mod
     logical             :: write_orb=.false.
     integer             :: iunit1=100,next,numbasef
     double precision    :: Rorb_max
+! Orbit-integrator tolerances used by find_bounce.  6 significant digits is far
+! tighter than the ~10% code-to-code benchmark and the eps~1d-3 class grid need,
+! while keeping orbit integration (the dominant cost) affordable.  The class
+! sampler relaxes these further for its grid build (see sample_class_doublecount).
+    double precision    :: fb_relerr=1d-6, fb_abserr=1d-8
+    double precision    :: fb_max_abs_delphi=huge(1d0)
+    double precision    :: fb_max_tau=huge(1d0)
+! Keep the outer X-point -> rho_pol-boundary classes (clip their search interval
+! to the |Delta_phi_b| resonance cap instead of dropping them) so resonances are
+! found out to the edge, matching NEO-RT.  Default off: keeping them integrates
+! many long-period orbits and needs the find_bounce low-energy speed-up first.
+    logical             :: clip_resonance_classes=.false.
   end module orbit_dim_mod
 !
 !------------------------------------------------------
@@ -111,194 +123,173 @@ subroutine find_bounce(next,velo_ext,dtau_in,z_eqm,taub,delphi,extraset,ierr)
 ! extraset(next) - extra integrals along the orbit (inout)
 ! ierr           - error flag, 0 = success, 1 = orbit left domain (output)
 !
-  use orbit_dim_mod, only : neqm,write_orb,iunit1,Rorb_max
+  use orbit_dim_mod, only : neqm,write_orb,iunit1,Rorb_max,fb_relerr, &
+                            fb_abserr,fb_max_abs_delphi,fb_max_tau
   use field_eq_mod, only : ierrfield
   use logging_mod, only : log_message
+  use fortnum_ode, only : ODE_EVENT_RISING
+  use fortnum_ode_vode, only : vode_state_t, vode_init, vode_integrate_to
+  use fortnum_status, only : fortnum_status_t, FORTNUM_OK
 !
   implicit none
 !
-! nousecut - way to close the orbit during primary search:
-! .true. - without using Poincare cut, can be used for general Phi distributions
-! .false. - with using the cut, valid for Phi=Phi(psi) only
-!  logical, parameter :: nousecut=.true.
-  logical, parameter :: nousecut=.false.
+! relative/absolute error of orbit integrator (module-settable; the sampler
+! relaxes them for speed, see orbit_dim_mod):
+  double precision :: relerr
+  double precision :: abserr
+! Step budget per orbit. A confined orbit closes in O(1e2..1e3) variable-order
+! Adams steps; even the longest kept orbits (|Delta_phi_b| up to ~1.5x the
+! resonance cap, tens of toroidal transits) close well under this. An orbit on
+! the trapped-passing separatrix has a diverging bounce time and never closes, so
+! cap the budget and skip it (ierr=1) instead of grinding to 50000 steps -- those
+! non-closing grazers are the bulk of the sampler cost and are discarded by the
+! class clip anyway.
+  integer, parameter :: max_orbit_steps=10000
+! Finite non-writing chunks let resonance sampling reject non-closing tails once
+! their toroidal displacement is beyond the largest requested resonance.
+  double precision, parameter :: nonwrite_chunk_factor=20.d0
 !
-! maximum number of Newton iterations for closing the orbit:
-  integer, parameter :: niter=20
-!
-! relative error of orbit integrator:
-  double precision, parameter :: relerr=1d-10 !8
-
   integer, intent(in) :: next
   double precision, intent(in) :: dtau_in
   double precision, dimension(neqm), intent(in) :: z_eqm
   double precision, dimension(next), intent(inout) :: extraset
   double precision, intent(out) :: taub, delphi
   integer, intent(out) :: ierr
-
-  logical :: firstpass
-  integer :: ndim, iter
-  double precision :: dtau
-  double precision :: dL2_pol,dL2_pol_start,dtau_newt,r_prev,z_prev
-  double precision :: tau0,RNorm,ZNorm,vnorm,dnorm,vel_pol,dL2_pol_min
-  double precision :: dZ_dR,sign_delZ,Z_tmp
-  double precision, dimension(neqm+next) :: z,z_start,vz
-  character(len=256) :: msg
+!
+  integer :: ndim, nleg
+  double precision :: t_root, dt_leg, tau_max, vel_pol, RNorm, ZNorm
+  double precision, dimension(neqm+next) :: z_start, vz
+  double precision, dimension(neqm+next) :: atol
+  double precision, dimension(:), allocatable :: yout
+  type(vode_state_t) :: vstate
+  type(fortnum_status_t) :: status
+  logical :: root_found
 !
   external velo_ext
 !
   ndim = neqm+next
   ierr = 0
-!
-  z(1:neqm)=z_eqm
-  if(next.gt.0) then
-    z(neqm+1:ndim)=extraset
-  endif
-  Rorb_max=z(1)
-!
-!  dtau=dtau_in/max(z(4),1d-3)
-  dtau=dtau_in
-!
-! Primary search:
-!
-  z_start=z
-!
-  call velo_ext(dtau,z,vz)
-!
-! unit 2D vector in the direction of the guiding center velocity in RZ-plane:
-!
-  vel_pol=sqrt(vz(1)**2+vz(3)**2)
-  RNorm=vz(1)/vel_pol
-  ZNorm=vz(3)/vel_pol
-!
-  tau0=0.d0
-!
-  if(write_orb) write (iunit1,*) z(1:neqm),vz(5)
-!
-! first step:
-!
-  call odeint_allroutines(z,ndim,tau0,dtau,relerr,velo_ext)
-  if(ierrfield.ne.0) then
-    write(msg,'(A,2ES14.6,A,2ES14.6)') &
-      'find_bounce: orbit left domain at R,Z=', &
-      z(1),z(3),' started from R,Z=',z_start(1),z_start(3)
-    write(*,'(A)') trim(msg)
-    call log_message(trim(msg))
-    ierr = 1
+  relerr = fb_relerr
+  abserr = fb_abserr
+  atol = abserr
+!
+  z_start(1:neqm)=z_eqm
+  if(next.gt.0) z_start(neqm+1:ndim)=extraset
+  Rorb_max=z_start(1)
+!
+! Bounce detection by RETURN TO THE START POINT, not by a Poincare-cut crossing.
+! The cut crossing fails for deeply trapped orbits (a small banana never reaches
+! the cut) -- exactly what stalled the old fixed-step search at low v/vth. The
+! orbit closes when its displacement projected on the start velocity direction
+! returns through zero moving in the start direction (a rising zero). This is
+! the original start-point closing condition and is well defined for trapped,
+! deeply trapped and passing orbits alike. fortnum's variable-order Adams
+! integrator locates the root on its own interpolant, so taub is exact and no
+! separate Newton closing is needed.
+  call velo_ext(0.d0, z_start, vz)
+  vel_pol = sqrt(vz(1)**2 + vz(3)**2)
+  RNorm = vz(1)/vel_pol
+  ZNorm = vz(3)/vel_pol
+!
+  call vode_init(vstate, ndim, 0.d0, z_start)
+  vstate%max_steps = max_orbit_steps
+  dt_leg = merge(dtau_in, nonwrite_chunk_factor*dtau_in, write_orb)
+  tau_max = min(2000.d0*dtau_in,fb_max_tau)
+!
+  if(write_orb) write (iunit1,*) z_start(1:neqm), 0.d0
+!
+! Step off the start plane once without event detection (the start sits exactly
+! on it: a rising zero at t=0 that must not be mistaken for the return).
+  call vode_integrate_to(rhs_w, vstate, dtau_in, relerr, atol, yout, status)
+  if(status%code /= FORTNUM_OK .or. ierrfield /= 0) then
+    call report_left(z_start, yout)
     return
   endif
-!
-  taub=dtau
-  if(nousecut) then
-! initialize sqrt(2) of the poloidal length of the step
-    dL2_pol=2.d0*(z(1)-z_start(1))**2+(z(3)-z_start(3))**2
-  else
-! initialize Poincare cut crossing check
-!
-    call get_poicut(z(1),Z_tmp,dZ_dR)
-!
-    sign_delZ=sign(1.d0,z(3)-Z_tmp)
-    firstpass=.true.
-    dL2_pol=2.d0*(z(1)-z_start(1))**2+(z(3)-z_start(3))**2
-  endif
-!
-  if(write_orb) then
-!
-    call velo_ext(dtau,z,vz)
-!
-    write (iunit1,*) z(1:neqm),vz(5)
-  endif
-!
-  do
-    r_prev=z(1)
-    z_prev=z(3)
-!
-    call odeint_allroutines(z,ndim,tau0,dtau,relerr,velo_ext)
-    if(ierrfield.ne.0) then
-      write(msg,'(A,2ES14.6,A,2ES14.6)') &
-        'find_bounce: orbit left domain at R,Z=', &
-        z(1),z(3),' started from R,Z=',z_start(1),z_start(3)
-      write(*,'(A)') trim(msg)
-      call log_message(trim(msg))
-      ierr = 1
+  Rorb_max=max(Rorb_max, yout(1))
+  if(write_orb) call write_orbit_point(yout)
+!
+! Integrate to the return. write_orb dumps the orbit, advancing in dtau_in-sized
+! legs and writing each; otherwise one far leg suffices since the event stops it
+! at the bounce.
+  root_found=.false.
+  do nleg=1,2*max_orbit_steps
+    call vode_integrate_to(rhs_w, vstate, vstate%tn+dt_leg, relerr, atol, &
+                           yout, status, event=return_event, &
+                           event_dir=ODE_EVENT_RISING, t_root=t_root, &
+                           root_found=root_found)
+    if(status%code /= FORTNUM_OK .or. ierrfield /= 0) then
+      call report_left(z_start, yout)
       return
     endif
-!
-    taub=taub+dtau
-    if(nousecut) then
-! check if poloidal distance to the starting point is smaller than
-! sqrt(2) of the poloidal length of the step:
-      dL2_pol=2.d0*((z(1)-r_prev)**2+(z(3)-z_prev)**2)
-      dL2_pol_start=(z(1)-z_start(1))**2+(z(3)-z_start(3))**2
-      if(dL2_pol_start.lt.dL2_pol) exit
-    else
-! check if Poincare cut has been crossed
-!
-      call get_poicut(z(1),Z_tmp,dZ_dR)
-!
-      if(sign_delZ*(z(3)-Z_tmp).lt.0.d0) then
-        if(firstpass) then
-! first crossing (continue integration)
-          firstpass=.false.
-          sign_delZ=-sign_delZ
-        else
-! second crossing (stop integration)
-          exit
-        endif
-      endif
+    Rorb_max=max(Rorb_max, yout(1))
+    if(write_orb) call write_orbit_point(yout)
+    if(root_found) exit
+    if(abs(yout(2)-z_start(2)).gt.fb_max_abs_delphi) then
+      call report_left(z_start, yout)
+      return
     endif
-!
-    if(write_orb) then
-!
-      call velo_ext(dtau,z,vz)
-!
-      write (iunit1,*) z(1:neqm),vz(5)
+    if(.not.write_orb) then
+      if(vstate%tn.ge.tau_max) then
+        call report_left(z_start, yout)
+        return
+      endif
     endif
-    Rorb_max=max(Rorb_max,z(1))
   enddo
+  if(.not.root_found) then
+    call report_left(z_start, yout)
+    return
+  endif
 !
-! End primary search
-!
-! Newton adjustment
-!
-  do iter=1,niter
+  taub = t_root
+  if(next.gt.0) extraset=yout(neqm+1:ndim)
+  delphi = yout(2)-z_start(2)
 !
-    call velo_ext(dtau,z,vz)
+  if(write_orb) write (iunit1,*) 'NaN NaN NaN NaN NaN NaN'
 !
-    vnorm=vz(1)*RNorm+vz(3)*ZNorm
-    dnorm=(z_start(1)-z(1))*RNorm+(z_start(3)-z(3))*ZNorm
-    if(dnorm**2.lt.dL2_pol*relerr) exit
-    dtau_newt=dnorm/vnorm
+  contains
 !
-    call odeint_allroutines(z,ndim,tau0,dtau_newt,relerr,velo_ext)
-    if(ierrfield.ne.0) then
+  subroutine rhs_w(t_, y_, dydt_, ctx_)
+    double precision, intent(in) :: t_
+    double precision, intent(in) :: y_(:)
+    double precision, intent(out) :: dydt_(:)
+    class(*), intent(in), optional :: ctx_
+    if(present(ctx_)) continue
+    call velo_ext(t_, y_, dydt_)
+  end subroutine rhs_w
+!
+  function return_event(t_, y_, ctx_) result(g)
+    double precision, intent(in) :: t_
+    double precision, intent(in) :: y_(:)
+    class(*), intent(in), optional :: ctx_
+    double precision :: g
+    if(present(ctx_)) continue
+    if(t_ < 0.d0) continue
+    g = (y_(1)-z_start(1))*RNorm + (y_(3)-z_start(3))*ZNorm
+  end function return_event
+!
+  subroutine write_orbit_point(y_)
+    double precision, intent(in) :: y_(:)
+    call velo_ext(0.d0, y_, vz)
+    write (iunit1,*) y_(1:neqm), vz(5)
+  end subroutine write_orbit_point
+!
+  subroutine report_left(zs, zc)
+    double precision, intent(in) :: zs(:), zc(:)
+    character(len=256) :: msg
+! Orbits grazing the X-point / trapped-passing separatrix do not close.  During
+! phase-space sampling this is a normal, expected outcome (the class is clipped
+! to the closing region), and there are many such orbits, so the per-orbit
+! string format + stdout write + log_message dominated the runtime (top self-time
+! in perf).  Emit the diagnostic only in explicit orbit-dump mode.
+    if(write_orb) then
       write(msg,'(A,2ES14.6,A,2ES14.6)') &
-        'find_bounce: orbit left domain at R,Z=', &
-        z(1),z(3),' started from R,Z=',z_start(1),z_start(3)
+        'find_bounce: orbit not closed from R,Z=', zs(1), zs(3), &
+        ' last R,Z=', zc(1), zc(3)
       write(*,'(A)') trim(msg)
       call log_message(trim(msg))
-      ierr = 1
-      return
     endif
-!
-    taub=taub+dtau_newt
-  enddo
-!
-  if(next.gt.0) then
-    extraset=z(neqm+1:ndim)
-  endif
-!
-  if(write_orb) then
-!
-    call velo_ext(dtau,z,vz)
-!
-    write (iunit1,*) z(1:neqm),vz(5)
-    write (iunit1,*) 'NaN NaN NaN NaN NaN NaN'
-  endif
-!
-! End Newton adjustment
-!
-  delphi=z(2)-z_start(2)
+    ierr = 1
+  end subroutine report_left
 !
   end subroutine find_bounce
 !
@@ -2395,7 +2386,7 @@ subroutine get_matrix_doublecount
 ! of class parameter x for adaptive refinement of interpolation grid over x
 !
   use sample_matrix_mod, only : n1,n2,x,amat
-  use orbit_dim_mod,     only : neqm,next
+  use orbit_dim_mod,     only : neqm,next,clip_resonance_classes
   use get_matrix_mod,    only : iclass
   use global_invariants, only : dtau,toten,perpinv
   use form_classes_doublecount_mod, only : ifuntype,R_class_beg,R_class_end,sigma_class
@@ -2404,6 +2395,7 @@ subroutine get_matrix_doublecount
 !
   logical :: fullbounce
   integer :: ierr,k
+  double precision, parameter :: twopi_gmd=6.28318530717958648d0
   double precision :: psiast,dpsiast_dRst,taub,delphi,xi,dxi_dx,Rst,sigma,delta_R, &
                       tau_fr,dphi_fr
   double precision, dimension(neqm) :: z
@@ -2421,8 +2413,15 @@ subroutine get_matrix_doublecount
   call starter_doublecount(toten,perpinv,sigma,Rst,   &
                            psiast,dpsiast_dRst,z,ierr)
 !
+! No valid orbit here (start off the cut or negative kinetic energy).  With the
+! class clip on, mark an out-of-cap sentinel so clip_bounds_to_dphi discards the
+! node; otherwise leave the original drop-on-fail behaviour.  No per-node print --
+! it floods the sampler.
   if(ierr.ne.0) then
-    print *,'get_matrix: error in starter'
+    amat(1,1)=0.d0
+    amat(2,1)=0.d0          ! Omega_b=0 : orbit does not close (tau_b -> infinity)
+    amat(3,1)=0.d0
+    if(next.gt.0) amat(4:3+next,1)=0.d0
     return
   endif
 !
@@ -2433,7 +2432,13 @@ subroutine get_matrix_doublecount
     if(fullbounce) then
 !
       call find_bounce(next,velo,dtau,z,taub,delphi,extraset,ierr)
-      if(ierr.ne.0) return
+      if(ierr.ne.0) then
+        amat(1,1)=psiast
+        amat(2,1)=0.d0
+        amat(3,1)=0.d0
+        if(next.gt.0) amat(4:3+next,1)=0.d0
+        return
+      endif
 !
     else
 !
@@ -2447,13 +2452,25 @@ subroutine get_matrix_doublecount
     extraset=0.d0
 !
     call find_bounce(next,velo_pphint,dtau,z,taub,delphi,extraset,ierr)
-    if(ierr.ne.0) return
+    if(ierr.ne.0) then
+      amat(1,1)=psiast
+      amat(2,1)=0.d0
+      amat(3,1)=0.d0
+      if(next.gt.0) amat(4:3+next,1)=0.d0
+      return
+    endif
 !
   endif
 !
+! Interpolate the canonical FREQUENCIES, not the bounce time / toroidal shift:
+! Omega_b = 2 pi / tau_b and Omega_phi = Delta_phi_b / tau_b stay finite and
+! smooth at the separatrix (Omega_b -> 0) where tau_b, Delta_phi_b diverge.  The
+! resonance m*Omega_b + n*Omega_phi = 0 has the same roots as Delta_phi_b + 2 pi
+! m/n = 0 (loss-free), but the polynomial grid converges instead of diverging.
+! get_rescond recovers tau_b = 2 pi / Omega_b.
   amat(1,1)=psiast
-  amat(2,1)=taub
-  amat(3,1)=delphi
+  amat(2,1)=twopi_gmd/taub
+  amat(3,1)=delphi/taub
 !
   if(next.gt.0) then
     amat(4:3+next,1)=extraset
@@ -2471,16 +2488,22 @@ subroutine sample_class_doublecount(iunit,ierr)
 !
   use sample_matrix_mod, only : nlagr,n1,n2,itermax,eps,xbeg,xend,  &
                                 npoi,xarr,amat_arr
-  use orbit_dim_mod,     only : next,numbasef
+  use orbit_dim_mod,     only : next,numbasef,fb_relerr,fb_abserr, &
+                                fb_max_abs_delphi,fb_max_tau, &
+                                clip_resonance_classes
   use get_matrix_mod,    only : relerror,relmargin,iclass
   use global_invariants, only : dtau,toten,perpinv,sigma
   use form_classes_doublecount_mod, only : ifuntype,R_class_beg,R_class_end
+  use potato_input_mod,  only : m_min,m_max,n_tor
   use cc_mod, only : dowrite
 !
   implicit none
 !
+  double precision, parameter :: twopi=6.28318530717958648d0
   integer :: iunit,ierr,i
-  double precision :: psiastbeg,psiastend,widthclass
+  double precision :: psiastbeg,psiastend,widthclass,dphi_cap
+  double precision :: fb_relerr_save,fb_abserr_save
+  double precision :: fb_max_abs_delphi_save,fb_max_tau_save
 !
   external :: get_matrix_doublecount
 !
@@ -2507,16 +2530,169 @@ subroutine sample_class_doublecount(iunit,ierr)
 !
   eps=relerror
 !
-  call sample_matrix(get_matrix_doublecount,ierr)
+! The grid interpolates to eps~1d-3, so the orbit integration needs no more than
+! that; relax the integrator tolerance for the (cost-dominant) sampling, then
+! restore the tight default for the accurate weight integration in pertham.
+  fb_relerr_save=fb_relerr
+  fb_abserr_save=fb_abserr
+  fb_max_abs_delphi_save=fb_max_abs_delphi
+  fb_max_tau_save=fb_max_tau
+  fb_relerr=1.d-4
+  fb_abserr=1.d-6
+  dphi_cap=0.d0   ! legacy, unused
+!
+! The class coordinate runs past the trapped-passing/X-point separatrix into a
+! region where the orbit no longer closes (find_bounce fails, Omega_b=0).  That
+! closing->non-closing step, and the tau_b divergence just before it, are what
+! made the adaptive grid refine without end.  Clip the endpoint to the closing
+! boundary so sample_matrix sees a smooth, all-closing interval.  Done with the
+! clip flag on or off; only the X-point classes (whose endpoint does not close)
+! pay for the boundary search.
+  call clip_class_to_closing(xbeg,xend)
+!
+  if(xend.gt.xbeg) then
+    call sample_matrix(get_matrix_doublecount,ierr)
+  else
+    ierr=1
+  endif
+!
+  fb_relerr=fb_relerr_save
+  fb_abserr=fb_abserr_save
+  fb_max_abs_delphi=fb_max_abs_delphi_save
+  fb_max_tau=fb_max_tau_save
 !
-  if(dowrite) then
+ if(dowrite) then
     print *,'npoi = ',npoi
     do i=1,npoi
       write(iunit,*) xarr(i),amat_arr(:,1,i)
     enddo
   endif
 !
-  end subroutine sample_class_doublecount
+ contains
+!
+ logical function class_closes(xx)
+! True if the orbit at class coordinate xx closes (Omega_b = amat(2,1) > 0).
+   use sample_matrix_mod, only : n1,n2,x,amat
+   implicit none
+   double precision, intent(in) :: xx
+   external :: get_matrix_doublecount
+   if(.not.allocated(amat)) allocate(amat(n1,n2))
+   x=xx
+   call get_matrix_doublecount
+   class_closes=amat(2,1).gt.0.d0
+ end function class_closes
+!
+ subroutine clip_class_to_closing(xb,xe)
+! Clip the separatrix-side endpoint to the closing boundary.  If the endpoint
+! orbit already closes the class is smooth -- nothing to do (cheap, the common
+! case).  Otherwise bisect [closing, non-closing] for the boundary x_sep and set
+! the endpoint just inside it, so the sampled interval is all-closing and smooth.
+   implicit none
+   double precision, intent(inout) :: xb,xe
+   double precision :: a,b,xm
+   integer :: it
+   if(xe.le.xb) return
+   if(class_closes(xe)) return            ! endpoint closes: smooth class, no clip
+   if(.not.class_closes(xb)) then         ! nothing closes: empty class
+     xe=xb
+     return
+   endif
+   a=xb                                   ! closes
+   b=xe                                   ! does not close
+   do it=1,20
+     xm=0.5d0*(a+b)
+     if(class_closes(xm)) then
+       a=xm
+     else
+       b=xm
+     endif
+   enddo
+   xe=a                                   ! last coordinate that still closes
+ end subroutine clip_class_to_closing
+!
+ subroutine preclip_class_to_dphi(xb,xe,dphi_cap)
+   use sample_matrix_mod, only : n1,n2,x,amat
+   implicit none
+   integer, parameter :: nprobe=9
+   integer :: i,ilo,ihi
+   double precision, intent(inout) :: xb,xe
+   double precision, intent(in) :: dphi_cap
+   double precision :: dx
+   double precision :: xprobe(nprobe)
+   double precision :: dphi_probe(nprobe)
+   external :: get_matrix_doublecount
+!
+   if(xe.le.xb) return
+   if(.not.allocated(amat)) allocate(amat(n1,n2))
+!
+   dx=(xe-xb)/dble(nprobe-1)
+   do i=1,nprobe
+     xprobe(i)=xb+dx*dble(i-1)
+     x=xprobe(i)
+     call get_matrix_doublecount
+     dphi_probe(i)=abs(amat(3,1))
+   enddo
+!
+   ilo=0
+   ihi=0
+   do i=1,nprobe
+     if(dphi_probe(i).le.dphi_cap) then
+       if(ilo.eq.0) ilo=i
+       ihi=i
+     endif
+   enddo
+!
+   if(ilo.eq.0) then
+     xe=xb
+     return
+   endif
+!
+   if(ilo.gt.1) ilo=ilo-1
+   if(ihi.lt.nprobe) ihi=ihi+1
+   xb=xprobe(ilo)
+   xe=xprobe(ihi)
+ end subroutine preclip_class_to_dphi
+!
+ end subroutine sample_class_doublecount
+!
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!
+  subroutine clip_bounds_to_dphi(xb,xe,dphi_cap)
+!
+! Shrinks the class search interval [xb,xe] to the span of the (partially
+! sampled) grid where |Delta_phi_b| <= dphi_cap, discarding the separatrix tail
+! where Delta_phi_b diverges.  Uses the grid (xarr) and matrix (amat_arr, row 3
+! = Delta_phi_b) left by the preceding sample_matrix call, so it needs no extra
+! orbit integration.  The physical (resonance-bearing) part is contiguous, so
+! the first and last in-cap grid nodes bracket it.
+!
+  use sample_matrix_mod, only : npoi,xarr,amat_arr
+!
+  implicit none
+!
+  double precision, intent(inout) :: xb,xe
+  double precision, intent(in)    :: dphi_cap
+  integer :: i,ilo,ihi
+!
+  ilo=0
+  ihi=0
+  do i=1,npoi
+    if(abs(amat_arr(3,1,i)).le.dphi_cap) then
+      if(ilo.eq.0) ilo=i
+      ihi=i
+    endif
+  enddo
+!
+  if(ilo.eq.0) then
+! no node within the cap: nothing to resolve, collapse the interval to skip it
+    xe=xb
+    return
+  endif
+!
+  xb=xarr(ilo)
+  xe=xarr(ihi)
+!
+  end subroutine clip_bounds_to_dphi
 !
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 !
diff --git a/POTATO/SRC/theta_rz_mod.f90 b/POTATO/SRC/theta_rz_mod.f90
deleted file mode 120000
index 91d48b83..00000000
--- a/POTATO/SRC/theta_rz_mod.f90
+++ /dev/null
@@ -1 +0,0 @@
-../../../libneo/src/magfie/theta_rz_mod.f90
\ No newline at end of file
diff --git a/POTATO/SRC/tt.f90 b/POTATO/SRC/tt.f90
index 842cd956..050aaada 100644
--- a/POTATO/SRC/tt.f90
+++ b/POTATO/SRC/tt.f90
@@ -2,7 +2,8 @@
   program classify_orbits
 !
   use parmot_mod,                   only : rmu,ro0
-  use orbit_dim_mod,                only : write_orb,iunit1,numbasef
+  use orbit_dim_mod,                only : write_orb,iunit1,numbasef,neqm, &
+                                          orbit_clip_resonance_classes => clip_resonance_classes
   use get_matrix_mod,               only : iclass
   use global_invariants,            only : dtau,toten,perpinv,cE_ref,Phi_eff
   use bounds_fixpoints_mod,         only : nregions
@@ -18,7 +19,13 @@ program classify_orbits
                                           rho_pol, rho_pol_max, scalfac_energy, &
                                           scalfac_efield, Rmax_orbit, ntimstep, &
                                           npoicut, m_min, m_max, n_tor, &
-                                          toten_plot, perpinv_plot, profile_file
+                                          toten_plot, perpinv_plot, profile_file, &
+                                          edge_extension, &
+                                          input_clip_resonance_classes => clip_resonance_classes, &
+                                          orbit_Rstart, orbit_Zstart, orbit_lambda, &
+                                          freq_Rmin, freq_Rmax, freq_n
+  use field_eq_mod,                 only : allow_sol, psi_axis, psi_sep
+  use field_sub,                    only : psif
 !
   implicit none
 !
@@ -30,10 +37,18 @@ program classify_orbits
   double precision,parameter  :: ev=1.6022d-12
 !
   logical :: classes_talk,plot_orbits,compute_equibrium_profiles,compute_resonant_torque,plot_poicut
+  logical :: trace_single_orbit,freq_scan
 !
   integer          :: ifdir_type,ierr,m,iunit,i,log_u
   double precision :: v0,bmod_ref
   double precision, dimension(:,:), allocatable :: resint
+  double precision, dimension(neqm) :: z_single
+  double precision :: taub_single,delphi_single,extraset_single(1)
+  external :: find_bounce, velo, magfie
+! frequency radial scan (itest_type=5):
+  integer          :: ifreq
+  double precision :: Rstep,omega_b,omega_phi,s_norm,rhopol
+  double precision :: bmod_d,sqrtg_d,bder_d(3),hcovar_d(3),hctrvr_d(3),hcurl_d(3),x_d(3)
 
 !
 !
@@ -43,9 +58,15 @@ program classify_orbits
   call read_potato_input('potato.in')
   call print_potato_input(6)
   if (is_logging_enabled()) call print_potato_input(log_u)
+!
+! Allow orbits to cross the separatrix into the scrape-off layer when requested:
+  allow_sol = edge_extension
+  orbit_clip_resonance_classes = input_clip_resonance_classes
 !
   iunit=71
 !
+  trace_single_orbit=.false.
+  freq_scan=.false.
   select case(itest_type)
   case(1)
     call tee_message('Plot orbits')
@@ -62,6 +83,18 @@ program classify_orbits
     plot_orbits=.false.
     compute_equibrium_profiles=.false.
     compute_resonant_torque=.true.
+  case(4)
+    call tee_message('Trace single orbit')
+    plot_orbits=.false.
+    compute_equibrium_profiles=.false.
+    compute_resonant_torque=.false.
+    trace_single_orbit=.true.
+  case(5)
+    call tee_message('Frequency radial scan')
+    plot_orbits=.false.
+    compute_equibrium_profiles=.false.
+    compute_resonant_torque=.false.
+    freq_scan=.true.
   case default
     call tee_message('unknown test')
     call close_logging()
@@ -150,7 +183,70 @@ program classify_orbits
 !
 !  call test_eqmagprofs  !WARNING: contains "stop" inside, comment this test out to continue
 !
-
+!.......................................
+!
+! Trace a single guiding-center orbit over one bounce period and write it.
+! With edge_extension=.true. the orbit may cross the separatrix into the SOL.
+!
+  if(trace_single_orbit) then
+    call tee_message('Tracing single orbit over one bounce')
+    z_single(1)=orbit_Rstart   ! R [cm]
+    z_single(2)=0.d0           ! phi
+    z_single(3)=orbit_Zstart   ! Z [cm]
+    z_single(4)=1.d0           ! p, monoenergetic at E_alpha
+    z_single(5)=orbit_lambda   ! pitch cosine v_par/v
+    write_orb=.true.
+    iunit1=100
+    extraset_single=0.d0
+    call find_bounce(0,velo,dtau,z_single,taub_single,delphi_single, &
+                     extraset_single,ierr)
+    close(iunit1)
+    write_orb=.false.
+    write(*,'(A,ES14.6,A,ES14.6,A,I0)') ' single orbit: taub=',taub_single, &
+          '  delphi=',delphi_single,'  ierr=',ierr
+  endif
+!
+!.......................................
+!
+! Frequency radial scan: at each midplane start radius trace one bounce and
+! record bounce frequency omega_b = 2 pi v0 / taub and toroidal precession
+! frequency omega_phi = delphi v0 / taub against the flux label rho_pol.
+!
+  if(freq_scan) then
+    call tee_message('Frequency radial scan')
+    write_orb=.false.
+    open(iunit,file='freq_scan.dat',status='replace',action='write')
+    write(iunit,'(A)') '# R_start[cm] rho_pol omega_b[1/s] omega_phi[1/s] '// &
+                       'taub delphi ierr'
+    Rstep=(freq_Rmax-freq_Rmin)/dble(max(freq_n-1,1))
+    do ifreq=1,freq_n
+      z_single(1)=freq_Rmin+Rstep*dble(ifreq-1)
+      z_single(2)=0.d0
+      z_single(3)=orbit_Zstart
+      z_single(4)=1.d0
+      z_single(5)=orbit_lambda
+! flux label of the start point: evaluate the field there, then normalize psi.
+      x_d=z_single(1:3)
+      call magfie(x_d,bmod_d,sqrtg_d,bder_d,hcovar_d,hctrvr_d,hcurl_d)
+      s_norm=(psif-psi_axis)/(psi_sep-psi_axis)
+      rhopol=sqrt(max(s_norm,0.d0))
+      extraset_single=0.d0
+      call find_bounce(0,velo,dtau,z_single,taub_single,delphi_single, &
+                       extraset_single,ierr)
+      if(ierr.eq.0 .and. taub_single.gt.0.d0) then
+        omega_b=2.d0*pi*v0/taub_single
+        omega_phi=delphi_single*v0/taub_single
+        write(iunit,'(6ES16.7,I3)') z_single(1),rhopol,omega_b,omega_phi, &
+                                    taub_single,delphi_single,ierr
+      else
+        write(iunit,'(6ES16.7,I3)') z_single(1),rhopol,0.d0,0.d0, &
+                                    taub_single,delphi_single,ierr
+      endif
+    enddo
+    close(iunit)
+    call tee_message('Frequency radial scan done -> freq_scan.dat')
+  endif
+!
 
 
 
diff --git a/test/test_freq_scan.f90 b/test/test_freq_scan.f90
new file mode 100644
index 00000000..be9d7ee3
--- /dev/null
+++ b/test/test_freq_scan.f90
@@ -0,0 +1,49 @@
+program test_freq_scan
+  ! Thin-orbit canonical frequencies vs flux surface for the rung-2 benchmark
+  ! against POTATO. At each s the bounce frequency Om_th and the toroidal
+  ! precession Om_ph are evaluated at a fixed pitch (eta = 0.84/Bmin, i.e.
+  ! v_par/v = 0.4 at the local low-field-side midplane, matching POTATO's
+  ! orbit_lambda) and a fixed energy. No E x B (Om_tE = 0) to match POTATO's
+  ! zero-potential profile, so Om_ph is the magnetic precession alone.
+  use iso_fortran_env, only: dp => real64
+  use do_magfie_mod, only: s, inp_swi, do_magfie_init, iota
+  use neort, only: init
+  use neort_freq, only: Om_th, Om_ph
+  use driftorbit, only: mph, mth, vth, etadt, Om_tE, sign_vpar, magdrift
+
+  implicit none
+
+  integer, parameter :: ns = 60
+  real(dp), parameter :: smin = 0.05_dp, smax = 0.95_dp
+  real(dp), parameter :: v_10keV = 9.787e7_dp   ! 10 keV deuteron, cm/s
+  real(dp), parameter :: lambda0 = 0.4_dp       ! v_par/v at the midplane
+
+  integer :: i, fid
+  real(dp) :: v, eta, Omth, dOmthdv, dOmthdeta, Omph, dOmphdv, dOmphdeta
+
+  inp_swi = 9          ! ASDEX Upgrade Boozer format
+  magdrift = .true.
+  s = smin
+  call do_magfie_init("in_file")
+
+  mth = 1
+  mph = 2
+  vth = v_10keV
+  v = v_10keV
+  Om_tE = 0.0_dp       ! no E x B, magnetic precession only
+  sign_vpar = 1.0_dp
+
+  open(newunit=fid, file='freq_scan_neort.dat')
+  write(fid, '(A)') '# s rho_tor omega_b[1/s] omega_phi[1/s] iota'
+  do i = 1, ns
+    s = smin + (smax - smin)*real(i - 1, dp)/real(ns - 1, dp)
+    call init
+    eta = (1.0_dp - lambda0**2)*etadt   ! = 0.84/Bmin, lambda=0.4 at midplane
+    call Om_th(v, eta, Omth, dOmthdv, dOmthdeta)
+    call Om_ph(v, eta, Omph, dOmphdv, dOmphdeta)
+    write(fid, '(5ES16.7)') s, sqrt(s), Omth, Omph, iota
+  end do
+  close(fid)
+  print *, 'freq_scan_neort.dat written, ns =', ns
+
+end program test_freq_scan