Topology (pyna.topo)

The pyna.topo package provides algorithms for analysing the topological structure of Poincaré maps: magnetic islands, X/O cycles, stable/unstable manifolds, and heteroclinic tangles.

---

Poincaré Maps

Multi-section Poincaré map infrastructure.

Provides an extensible framework for collecting field-line crossings on arbitrary 2-D surfaces in (R, Z, φ) space.

Rename history

• PoincareMap → PoincareAccumulator (conflicts with pyna.topo.dynamics.PoincareMap)

Backward-compatible aliases at module level keep old import paths working:

from pyna.topo.poincare import PoincareMap        # still works
from pyna.topo.poincare import ToroidalSection    # still works

Classes

PoincareSection — abstract base for crossing surfaces (poincare-accumulator use) PoincareToroidalSection — φ = φ₀ plane, implements detect_crossing PoincareAccumulator — accumulate crossings on multiple sections

Functions

poincare_from_fieldlines — trace field lines and collect crossings

class pyna.topo.poincare.PoincareAccumulator(sections)

Bases: object

Collect and store trajectory crossings on multiple sections.

Parameters:

sections (list of Section) – The crossing surfaces to monitor.

crossing_array(section_idx)

Return all crossings for a given section as an (N, 3) array.

Parameters:

section_idx (int) – Index into the sections list.

Return type:

ndarray

Returns:

ndarray of shape (N, 3) — columns (R, Z, φ). – Returns shape (0, 3) if no crossings.

record_step(pt_prev, pt_curr)

Check all sections for a crossing between two consecutive points.

Return type:

None

record_trajectory(traj)

Record all crossings in a full trajectory array.

Parameters:

traj (ndarray of shape (N, 3)) – Trajectory with columns (R, Z, φ).

Return type:

None

pyna.topo.poincare.PoincareMap

Alias for PoincareAccumulator kept for backward compatibility.

class pyna.topo.poincare.PoincareSection

Bases: ABC

Abstract crossing surface for use with PoincareAccumulator.

Concrete subclasses define any 2-D surface in 3-D (R, Z, φ) space. The detect_crossing() method uses linear interpolation between consecutive trajectory points to find where the trajectory crosses this surface.

abstractmethod detect_crossing(pt_prev, pt_curr)

Return the interpolated crossing point (R, Z, φ), or None.

Parameters:

• pt_prev (ndarray of shape (3,)) – Consecutive trajectory points (R, Z, φ).

• pt_curr (ndarray of shape (3,)) – Consecutive trajectory points (R, Z, φ).

Return type:

Optional[ndarray]

Returns:

ndarray of shape (3,) or None

abstract property label: str

Human-readable label for this section.

class pyna.topo.poincare.PoincareToroidalSection(phi0=0.0)

Bases: PoincareSection

Toroidal section at φ = φ₀ for use with PoincareAccumulator.

Detects when the toroidal angle φ crosses φ₀ (modulo 2π), using linear interpolation to find the precise crossing point.

Parameters:

phi0 (float) – The toroidal angle of the section (radians). Default 0.0.

detect_crossing(pt_prev, pt_curr)

Detect when φ crosses φ₀ (in the direction of increasing φ).

Return type:

Optional[ndarray]

property label: str

Human-readable label for this section.

pyna.topo.poincare.Section

Backward-compatible alias: Section from poincare module → PoincareSection ABC For the canonical Section ABC, use pyna.topo.section.Section.

pyna.topo.poincare.ToroidalSection

Backward-compatible alias: ToroidalSection from poincare module → PoincareToroidalSection (has detect_crossing for use with PoincareAccumulator) For section-based topology queries, use pyna.topo.section.ToroidalSection.

pyna.topo.poincare.poincare_from_fieldlines(field_func, start_pts, sections, t_max, dt=0.04, backend=None)

Trace field lines and collect Poincaré crossings.

Parameters:

• field_func (callable) – field_func(rzphi) → (dR, dZ, dphi) unit-tangent ODE rhs.

• start_pts (ndarray of shape (N, 3)) – Starting points (R, Z, φ).

• sections (list of Section) – Sections on which to collect crossings.

• t_max (float) – Maximum arc-length parameter for each field line.

• dt (float) – Step size for the RK4 integrator.

• backend (FieldLineTracer or None) – If None, a pyna.flt.FieldLineTracer is created.

Return type:

PoincareAccumulator

Returns:

PoincareAccumulator

pyna.topo.poincare.rotational_transform_from_trajectory(traj, axis_RZ=None, n_turns=None)

Estimate the rotational transform ι (or safety factor q = 1/ι) from a traced field-line trajectory.

The rotational transform is the average poloidal angle advance per toroidal turn. For a closed (periodic) orbit on a flux surface with safety factor q = m/n, we have ι = 1/q = n/m.

Algorithm

1. Compute the cumulative poloidal angle Θ(s) of the trajectory around the magnetic axis axis_RZ.

2. Compute the cumulative toroidal angle Φ(s) from the φ column of the trajectory.

3. Fit a linear relationship Θ ≈ ι · Φ by least squares. The slope ι is the rotational transform.

type traj:

ndarray

param traj:

Trajectory (R, Z, φ) as returned by pyna.flt.FieldLineTracer.trace().

type traj:

ndarray, shape (N, 3)

type axis_RZ:

Optional[ndarray]

param axis_RZ:

Approximate position [R, Z] of the magnetic axis (or any reference point inside the flux surface). If None, the centroid of the trajectory projection onto the (R, Z) plane is used.

type axis_RZ:

array_like of shape (2,) or None

type n_turns:

Optional[int]

param n_turns:

If provided, restrict the computation to the first n_turns toroidal traversals. If None, use the full trajectory.

type n_turns:

int or None

rtype:

float

returns:

iota (float) – Estimated rotational transform ι = dΘ/dΦ. The safety factor is q = 1 / ι.

Notes

For a regular (KAM) flux surface, the returned value converges as the trajectory length increases. For a chaotic trajectory, the “effective” rotational transform still gives a useful indicator of the local winding rate, though it will not converge to a rational number.

Examples

>>> iota = rotational_transform_from_trajectory(traj, axis_RZ=[1.0, 0.0])
>>> q = 1.0 / iota   # safety factor

---

Magnetic Islands

pyna.topo.toroidal_island — backward-compat shim.

Re-exports utilities from pyna.topo.toroidal._utils.

pyna.topo.toroidal_island.all_rational_q(m_max, n_max, q_min=None, q_max=None)

Enumerate all unique q = m/n rational values (m, n > 0).

Return type:

List[List[List[int]]]

pyna.topo.toroidal_island.island_halfwidth(m, n, S_res, S, q_profile, tilde_b_mn, tilde_b_mn_index=None)

Estimate the half-width of a magnetic island at a rational surface.

Return type:

float

pyna.topo.toroidal_island.locate_all_rational_surfaces(S, q_profile, m_max=12, n_max=3, s=0.01)

Find all rational surfaces q = m/n for |m| ≤ m_max, 1 ≤ n ≤ n_max.

Return type:

Dict[int, Dict[int, List[float]]]

pyna.topo.toroidal_island.locate_rational_surface(S, q_profile, m, n, s=0.01)

Find the S locations where q(S) = m/n.

Return type:

List[float]

---

Manifold Computation

pyna.topo.manifold.accumulate_s_from_RZ_arr(W1d_RZ)

pyna.topo.manifold.create_W1d_interpolator_s_to_RZdRZds(afield, Wivp_bundle, phi, phi_epsilon, mturn=<class 'int'>)

pyna.topo.manifold.grow_manifold_from_Xcycle_eig_interp(afield, Xcycle_RZdiff, Jac_evosol_along_Xcycle, eigind, S_span, S_num, Phi_span, Phi_num, rev_eigvec=False, first_step=5e-05, max_step=0.0001)

pyna.topo.manifold.grow_manifold_from_Xcycle_naive_carousel(afield, Xcycle_RZdiff, Jac_evosol_along_Xcycle, eigind, Phi_span, W_nPhi, Ind_num, first_step=5e-05, max_step=0.0001)

pyna.topo.manifold.grow_manifold_from_Xcycle_naive_init_segment(afield, Xcycle_RZdiff, Jac_evosol_along_Xcycle, eigind, Phi_span, total_deltaPhi, ptsnum_initseg=300, initseg_len=8e-05, *arg_sovle_ivp, **kwarg_solve_ivp)

pyna.topo.manifold.smooth1D(x, window_len=11, window='hanning')

smooth the data using a window with requested size.

This method is based on the convolution of a scaled window with the signal. The signal is prepared by introducing reflected copies of the signal (with the window size) in both ends so that transient parts are minimized in the begining and end part of the output signal. This function is borrowed from https://scipy-cookbook.readthedocs.io/items/SignalSmooth.html

input:

x: the input signal window_len: the dimension of the smoothing window; should be an odd integer window: the type of window from ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’

flat window will produce a moving average smoothing.

output:

the smoothed signal

example:

t=linspace(-2,2,0.1) x=sin(t)+randn(len(t))*0.1 y=smooth(x)

see also:

numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve scipy.signal.lfilter

TODO: the window parameter could be the window itself if an array instead of a string NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.

---

Periodic Orbits / Cycles

Toroidal periodic-orbit traces and Newton solvers.

This module is explicitly for toroidal field-line dynamics parameterised by phi. It stores sampled toroidal periodic-orbit traces and provides Newton search utilities in (R, Z, phi).

class pyna.topo.toroidal_cycle.ToroidalPeriodicOrbitTrace(rzphi0, period_m, trajectory, DPm)

Bases: object

A periodic field-line orbit.

rzphi0

Starting point (R, Z, phi0).

Type:

ndarray, shape (3,)

period_m

Number of toroidal turns (period of Poincaré map). Corresponds to m in q=m/n notation: P^m(x) = x.

Type:

int

trajectory

Full orbit trajectory (R, Z, phi).

Type:

ndarray, shape (N, 3)

DPm

Monodromy matrix DPm = DX_pol(phi_end), phi_end = 2π·m. Eigenvalues characterize stability.

Type:

ndarray, shape (2, 2)

DPm: ndarray

property Jac: ndarray

Deprecated alias for DPm.

diagnostics()

Return type:

dict

property eigenvalues: ndarray

property is_stable: bool

True if |eigenvalues| ≤ 1 (elliptic, O-point type).

property kind: str

‘O’ if elliptic (stable), ‘X’ if hyperbolic.

period_m: int

property period_n: int

Alias for period_m (backward compatibility).

rzphi0: ndarray

section_cut(section)

Return the point(s) at which this orbit crosses a toroidal section.

Parameters:

section (float | ToroidalSection) – Toroidal section angle (radians) or a concrete ToroidalSection.

Return type:

list

Returns:

• list of dict, each with keys 'R', 'Z', 'phi', 'kind'.

• Returns an empty list if the trajectory is not available or no

• crossing is found within tolerance.

property stability_index: float

Tr(DPm)/2 for a 2x2 symplectic map. |k|<1 → elliptic, |k|>1 → hyperbolic.

trajectory: ndarray

pyna.topo.toroidal_cycle.find_all_cycles_near_resonance(field_func, equilibrium, m, n, n_seeds=8, dt=0.05, RZlimit=None)

Find all O- and X-point cycles near the q=m/n resonant surface.

For a q = m/n resonance, there are m O-points and m X-points at equally spaced angular positions around the resonant surface. This function seeds the Newton-Raphson solver at 2·m·n_seeds angular positions around the resonant flux surface and deduplicates the resulting orbits.

Parameters:

• field_func (callable)

• equilibrium (StellaratorSimple or similar) – Must have resonant_psi(m, n) and R0, r0 attributes.

• m (int) – Mode numbers defining the resonance q = m/n.

• n (int) – Mode numbers defining the resonance q = m/n.

• n_seeds (int) – Number of angular seeds per expected fixed point.

• dt (float) – Integration step.

• RZlimit (tuple or None)

Return type:

list

Returns:

list of ToroidalPeriodicOrbitTrace

pyna.topo.toroidal_cycle.find_cycle(field_func, init_rzphi, n_turns=1, dt=0.05, RZlimit=None, max_iter=50, tol=1e-08, n_fallback_seeds=12, fallback_radius=0.05)

Find a periodic orbit starting from init_rzphi using Newton-Raphson.

G(x0) = P^n(x0) - x0 = 0

If Newton diverges or leaves domain, automatically tries n_fallback_seeds alternative starting points distributed on a circle of radius fallback_radius around init_rzphi.

Parameters:

• field_func (callable) – Field function f(rzphi) → (dR/dl, dZ/dl, dphi/dl).

• init_rzphi (array (3,)) – Initial guess (R0, Z0, phi0).

• n_turns (int) – Period (number of toroidal turns).

• dt (float) – Integration step size in φ.

• RZlimit (tuple or None) – Domain limits (R_min, R_max, Z_min, Z_max).

• max_iter (int) – Maximum Newton iterations.

• tol (float) – Convergence tolerance on |G(x0)|.

• n_fallback_seeds (int) – Number of fallback seeds to try if primary Newton fails.

• fallback_radius (float) – Radius around init_rzphi for fallback seeds.

Return type:

Optional[ToroidalPeriodicOrbitTrace]

Returns:

ToroidalPeriodicOrbitTrace or None if not found.

pyna.topo.toroidal_cycle.jacobian_of_poincare_map(field_func, rzphi0, n_turns, dt=0.05, eps=1e-05)

Finite-difference Jacobian ∂(R_f, Z_f)/∂(R_0, Z_0) of n-turn Poincaré map.

Parameters:

• field_func (callable)

• rzphi0 (array-like (3,))

• n_turns (int)

• dt (float) – Integration step size.

• eps (float) – Finite-difference perturbation.

Return type:

ndarray

Returns:

DPm (ndarray, shape (2, 2)) – Monodromy matrix (Jacobian of the n-turn Poincaré map). det(DPm) ≈ 1 for area-preserving.

pyna.topo.toroidal_cycle.poincare_map_n(field_func, rzphi0, n_turns, dt=0.05, RZlimit=None)

Integrate field line for n toroidal turns, return final (R, Z).

The integration uses the toroidal angle φ as the independent variable.

Parameters:

• field_func (callable) – field_func(rzphi) → (dR/dl, dZ/dl, dphi/dl).

• rzphi0 (array-like (3,)) – Starting point (R, Z, phi).

• n_turns (int) – Number of toroidal turns to integrate (φ increases by 2π·n_turns).

• dt (float) – Step size in φ (radians) for the integrator.

• RZlimit (tuple or None) – Optional (R_min, R_max, Z_min, Z_max) domain boundary.

Return type:

Tuple[float, float]

Returns:

(R_final, Z_final) or (nan, nan) if field left domain.

pyna.topo.toroidal_cycle.poincare_map_n_trajectory(field_func, rzphi0, n_turns, dt=0.05, RZlimit=None)

Integrate field line for n_turns, return full (R, Z, phi) trajectory.

Return type:

ndarray