Restore liberal Newton continuation outside SPECTRE

DOC/config.md

@@ -12,20 +12,22 @@
   values to refine event location independently of the nonlinear solve.
 
 * `symplectic_newton_warning_mode` defaults to `.true.`. When an implicit solve
-  reaches its iteration limit, SIMPLE continues only if the final Newton
-  correction is finite and no more than ten times the requested relative
-  tolerance. Each occurrence is reported by the corresponding `*_maxit`
-  diagnostic. The production RK, symplectic, and full-orbit paths then use the
-  same terminal convention: if bounded recovery cannot resolve a numerical
-  microstep, SIMPLE retains any contiguous accepted prefix and holds only the
-  unresolved remainder, records `warning_step_skip`, advances the clock for
-  the complete interval, and retries from that valid state on the next
-  microstep. A warning hold does not terminate or
-  numerically disqualify the marker: a marker that reaches the requested end
-  time remains a resolved survivor, while a later physical boundary event is
-  still a loss. Set the option to `.false.` for strict diagnostic runs that end
-  only the affected marker at the first exhausted recovery and report a
-  101--105 `orbit_exit_code` with `NaN` in `times_lost`.
+  reaches its iteration limit, SIMPLE commits any finite final Newton iterate
+  and continues the orbit. This default path uses the historical undamped
+  Newton correction; the radial solver-basin limiter remains active in strict
+  mode. Each occurrence is reported by the corresponding `*_maxit` diagnostic.
+  Recursive recovery is reserved for an unusable step, such as a failed linear
+  solve or non-finite iterate. The production RK, symplectic, and full-orbit
+  paths use the same terminal convention: if that recovery cannot resolve a
+  numerical microstep, SIMPLE retains any contiguous accepted prefix and holds
+  only the unresolved remainder, records
+  `warning_step_skip`, advances the clock for the complete interval, and
+  retries from that valid state on the next microstep. A warning hold does not
+  terminate or numerically disqualify the marker: a marker that reaches the
+  requested end time remains a resolved survivor, while a later physical
+  boundary event is still a loss. Set the option to `.false.` for strict
+  diagnostic runs that end only the affected marker at the first exhausted
+  recovery and report a 101--105 `orbit_exit_code` with `NaN` in `times_lost`.
 
 * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi`
   control the Meiss or Albert canonical-map grid. Their defaults are 62, 63,

DOC/coordinates-and-fields.md

@@ -1200,7 +1200,29 @@ because they have no reference-coordinate geometry.
 | `src/magfie.f90` | Unified field evaluation interface |
 | `src/magfie_can_boozer.f90` | Boozer/Canflux implementations |
 
-### 10.5 Integration
+### 10.5 Canonical Frequency API
+
+`src/orbit_frequencies.f90` provides the source-level Fortran interface
+`compute_canonical_frequencies(tracer, initial_state, options, result)`.
+`initial_state` uses SIMPLE integration coordinates in the standard ordering
+`[s, theta, phi, v/v0, lambda]`; all angles are radians and remain unwrapped.
+The caller must initialize `tracer_t` and set `dtaumin`, `v0`, `relerr`, and
+`integmode`. The routine copies the tracer before integration and therefore
+does not advance caller-owned orbit state.
+
+The result contains the positive bounce/transit frequency
+`omega_b = 2*pi/period` and the signed mean toroidal frequency
+`omega_phi = delta_phi/period`. Periods are seconds, frequencies are rad/s,
+and toroidal displacement is radians. Named status and orbit-class constants
+are exported by the module. A downstream Fortran target can link the CMake
+target `simple` and `use orbit_frequencies`; the target publishes SIMPLE's
+build-tree module directory.
+
+The flat wrapper in `simple_main` accepts the public reference-coordinate
+state, performs `ref_to_integ`, and is the boundary used by f90wrap and
+`pysimple.compute_canonical_frequencies`.
+
+### 10.6 Integration
 
 | File | Purpose |
 |------|---------|

python/CMakeLists.txt

@@ -73,6 +73,7 @@ message(STATUS "Python binary output dir: ${CMAKE_CURRENT_BINARY_DIR}")
 set(FILES_TO_WRAP
     simple_main.f90
     simple.f90
+    orbit_frequencies.f90
     samplers.f90
     cut_detector.f90
     callback.f90

python/README.md

@@ -52,6 +52,51 @@ First verify the bindings import cleanly:
 Complete examples are in ``examples/simple_api.py``, ``examples/classify_fast.py``,
 and ``examples/classify_fractal.py``.
 
+Canonical orbit frequencies
+---------------------------
+
+``compute_canonical_frequencies`` traces one particle until it completes the
+requested number of poloidal cycles:
+
+.. code-block:: python
+
+    position = [0.6, 0.0, 0.25 * np.pi, 1.0, 0.3]
+    frequency = pysimple.compute_canonical_frequencies(
+        position, integrator="midpoint", n_periods=4
+    )
+    print(frequency["omega_b"], frequency["omega_phi"])
+
+For trapped particles, a cycle is measured between same-direction
+``v_parallel = 0`` crossings. For passing particles, it is measured between
+successive signed ``2*pi`` advances of the unwrapped poloidal angle. The
+toroidal displacement is also unwrapped. Periods are seconds, displacements
+are radians, and both frequencies are angular frequencies in rad/s.
+
+Use ``n_periods=1`` for the inexpensive single-cycle value in an axisymmetric
+field. Larger values return the mean and sample spread over several cycles,
+which is useful in asymmetric fields. Always inspect ``status`` before using a
+result; losses, integrator errors, and the maximum-step limit remain visible
+instead of being silently discarded. An orbit that leaves the plasma through
+the outer radial boundary reports ``FREQ_ORBIT_LOST``.
+
+The particle species and energy follow the same defaults as ``simple.in``
+(3.5 MeV alphas). A frequency computed for a different particle needs the
+matching overrides in ``init``, for example a 5 keV deuteron:
+
+.. code-block:: python
+
+    pysimple.init("field.nc", deterministic=True, n_e=1, n_d=2, facE_al=700.0)
+
+A fast particle in a small equilibrium can be genuinely unconfined; a
+``FREQ_ORBIT_LOST`` result for a seed that should be confined usually means
+the species/energy overrides are missing.
+
+The same event machinery is available as ``pysimple.trace_to_cut``. Select
+``cut="tip"`` for a same-direction ``v_parallel=0`` section or
+``cut="toroidal"`` for a field-period section. The returned state is converted
+back to public reference coordinates. This driver currently requires one of
+SIMPLE's symplectic integrators.
+
 Legacy script note
 ------------------