Closed orbit finder for OPALX

[aaadelmann]

[aaadelmann]

Proposal: Add closed orbit finder for OPALX

Summary

Add a new command block with the same outer-command/sub-command semantics as TRACK:

COF, LINE=myRing, BEAM=beam1,
     DT=1.0e-12, TURNS=200000, TIMEINTEGRATOR=LF2;
    RUN, METHOD=NEWTON, DIMENSION=4,
         X=0.0, PX=0.0, Y=0.0, PY=0.0,
         MAXIT=20, XTOL=1.0e-10, PTOL=1.0e-10;
ENDCOF;

LINE deliberately remains the attribute name for compatibility with TRACK, but the selected BeamSequence must have the RING topology proposed in #543. COF does not make a LINE cyclic, and RING itself remains topology metadata only.

The first implementation should find the four-dimensional, fixed-energy, single-particle closed orbit of a static magnetic ring
representing a Cyclotron. We then can check easy with OLD OPAL.

RF/6D closure and collective closed orbits should be separate extensions.

Relationship to #543

Issue #543 gives a RING the same list grammar and beamline representation as a LINE, while attaching topology and canonical ring-name metadata to its occurrences. It explicitly does not add repeated tracking, an end-to-start connection, or geometric closure.

The division of responsibility should be:

• RING declares topology and supplies the elements.
• COF verifies the additional geometric and physics preconditions needed by a solver.
• A one-turn tracking advances a test particle exactly once and expresses the exit state in the entrance frame.
• The nonlinear solver finds a fixed point of that map.

This avoids making parsing depend on a tracking algorithm and avoids overloading ElementType::RING, which already denotes the cyclotron ring-field element.

Note

ElementType::RING maybe need to revisited

What is retained from the old OPAL cyclotron COF

The old Distribution/ClosedOrbitFinder.h implements Gordon's method for the special cyclotron state (r, p_r). It integrates the orbit and variational equations, forms the end-minus-start residual, applies a Newton correction using the one-turn matrix, optionally continues from low energy to the requested energy, and obtains tunes from the resulting linear motion.

Retain these ideas:

1. Define closure as a one-turn fixed point.
2. Compute a one-turn linear matrix.
3. Use the one-turn matrix in a Newton correction.
4. Support energy continuation and extrapolate the next initial guess from prior solutions.
5. Report the residual, orbit samples, one-turn matrix, stability, and tunes.

At the end, we have to revisit the Cyclotron specifics w.r.t. OPALX general coordinate system:

• polar coordinates and the (r, p_r)-only correction;
• direct calls to Cyclotron::apply() and field-map sector averaging;
• dependence on FMLOWE, FMHIGHE, cyclotron harmonic number, or the first RF frequency;
• the independent Odeint RK4 integration;
• tune calculation by zero-crossing counts.

The new COF must use the same element placement, field evaluation, reference-particle conventions, and pusher as TRACK, otherwise a mathematically closed COF orbit may not be closed when subsequently tracked.

Mathematical definition (c.c. from old OPAL)

Note this is a bit explicite because in old OPAL the cyclotron flavour uses different conventions and units.

At the entrance Poincare section use local coordinates

u = (x, px, y, py)

Positions are in metres and momenta use OPALX's normalized kinetic momentum beta*gamma. The beam supplies rest mass, charge, and fixed reference momentum p0. For a trial state, place the longitudinal coordinate on the entrance section and construct forward momentum as

pz = sqrt(p0^2 - px^2 - py^2).

Let T_E(u) be the state after one ordered traversal of the selected ring, transformed from the nominal exit frame back into the entrance frame. The closed orbit satisfies

F(u) = T_E(u) - u = 0.

If M = dT_E/du, the direct generalization of the old Gordon correction is

(M - I) delta_u = -F(u),
u_next = u + lambda * delta_u,

where lambda is selected by a backtracking line search. Solve with rank-revealing QR or SVD, not an explicit inverse. Report a near-singular (M-I) and its condition estimate; it commonly signals an underconstrained coordinate or a tune near an integer.

Convergence should be componentwise:

max(|delta x|, |delta y|, |Fx|, |Fy|)       <= XTOL
max(|delta px|, |delta py|, |Fpx|, |Fpy|)   <= PTOL

It must also require a successful, non-lost one-turn integration. Verify the final result with one fresh, unperturbed traversal.

One-turn map and TRACK consistency

Introduce a small OneTurnMap/ReferenceParticleStepper service under src/Algorithms. It should be shared by OrbitThreader, COF, and eventually other optics tools. It should:

• build the runtime OpalBeamline through the normal visitor path;
• use the same ELEMEDGE and explicit-6D placement transforms as TRACK;
• call ElementBase::applyToReferenceParticle() for external fields;
• use the same Boris/LF2 reference-particle stepping path as current ParallelTracker;
• honor the selected TRACK time step, maximum-step guard, and element step-size checks;
• stop on the ring exit Poincare section rather than ZSTOP or approximate elapsed time;
• return the exit state in the ring entrance frame;
• optionally retain orbit samples and active-element names for diagnostics.

Do not implement an iteration by constructing a PartBunch, field solver, H5 sink, and ParallelTracker for every trial. That would couple deterministic single-particle optics to emission, space charge, output, restart, and global OpalData state.

The current TIMEINTEGRATOR interface and actual reference-particle pusher must be reconciled in the shared service. Initially it is acceptable to support only the pusher actually used by TRACK and explicitly reject other values; silently using a COF-only RK4 would repeat the legacy inconsistency.

One-turn Map

For v1 compute M using central finite differences of the identical one-turn map:

M[:,i] = (T_E(u + h_i e_i) - T_E(u - h_i e_i)) / (2 h_i).

This costs nine one-turn traversals per Newton iteration in 4D, but is simple, testable, and guarantees the derivative corresponds to the actual OPALX pusher. Select h_i from coordinate-specific absolute floors and relative scales; expose expert overrides as FDSTEP={hx,hpx,hy,hpy}.

A later tangent-map or automatic-differentiation implementation can replace finite differences behind the same interface after comparison with the finite-difference reference. Obtain tunes from eigenvalues of the final M. Stable complex-conjugate pairs give nu = arg(lambda)/(2*pi); report unstable real/reciprocal pairs rather than forcing a real tune.

Command model

COF outer command

Reuse common TRACK meanings and validation:

• LINE (required): selected sequence; it must be a RING and its owner name must match.
• BEAM (required in v1): species and design momentum/energy. Reject BEAMS because a fixed point belongs to one species and one reference rigidity.
• DT, MAXSTEPS, T0, TIMEINTEGRATOR: identical units and numerical meanings to TRACK. V1 requires exactly one DT and one MAXSTEPS; MAXSTEPS is a failure guard while successful traversal stops at the ring boundary.

ZSTART, ZSTOP, SOURCES, DTSCINIT, DTAU, and MAP_ORDER should not be accepted in v1. Reject user-set values precisely rather than ignoring them.

COF-local RUN

Use a separate CofRun registered as RUN only in CofParser; do not overload TrackRun with solver attributes.

Attribute	Default	Meaning
METHOD	NEWTON	Fixed-point solver
DIMENSION	4	Only 4 accepted initially
X, PX, Y, PY	0	Entrance-frame initial guess
MAXIT	20	Maximum accepted Newton iterations
XTOL	1e-10 m	Position correction/residual tolerance
PTOL	1e-10	Normalized-momentum tolerance
JACOBIAN	CENTRAL	One-turn map construction
FDSTEP	automatic	Four coordinate perturbations
DAMPING	TRUE	Backtracking line search
ENERGYSTART	unset	Optional continuation start
DENERGY	automatic	Maximum continuation increment
GEOMTOL	1e-9 m	Exit-to-entrance position tolerance
ANGLETOL	1e-10 rad	Exit-to-entrance orientation tolerance
OUTPUT	input stem + _cof	Diagnostic output stem

ENDCOF stops the local parser exactly as ENDTRACK stops TRACK.

Avoid a static Cof::block if practical. CofCmd::execute() should own a scoped context and give the local parser a non-owning reference. If infrastructure forces the legacy static pattern, contain it behind RAII so exceptions cannot leave a dangling global context.

Preconditions and unsupported physics in v1

Before iterating, check:

• the selected root sequence has BeamlineTopology::RING;
• the ring is non-empty and has one unambiguous entrance/exit boundary;
• nominal exit position and orientation close onto the entrance within separately reported tolerances;
• the beam has nonzero charge, positive mass, and defined momentum/energy;
• the ring contains no time-dependent RF/electric element;
• all field maps initialize successfully;
• the initial transverse momentum permits a real forward pz.

Space charge, wakes, particle-matter interactions, emission, stochastic processes, decays, and self-consistent collective closed orbits are disabled. Static errors, misalignments, fringe fields, and overlapping external fields are included because they belong to the one-turn map. Aperture/material interception invalidates the trial and must identify the responsible element.

Geometry checking belongs in COF, not RING: #543 intentionally permits a topological ring without asserting geometry.

Result and output

On success print and write:

• entrance coordinates and beam energy/momentum;
• iterations, damping factors, residuals, corrections, and condition estimates;
• verified closure residual from a fresh turn;
• the 4x4 one-turn matrix, eigenvalues, stability, and transverse tunes;
• sampled s, lab/local positions, normalized momentum, time, and active element;
• integration settings and ring/beam names for reproducibility.

Use a lightweight text/SDDS orbit file plus a machine-readable summary. COF must not silently mutate BEAM, DISTRIBUTION, or later TRACK initial conditions. A separately designed named CLOSEDORBIT result can later provide explicit reuse.

Proposed source layout

src/COF/CofCmd.{h,cpp}
src/COF/CofContext.{h,cpp}
src/COF/CofParser.{h,cpp}
src/COF/CofRun.{h,cpp}
src/COF/CofEnd.{h,cpp}
src/Algorithms/OneTurnMap.{h,cpp}
src/Algorithms/ClosedOrbitSolver.{h,cpp}
src/Algorithms/ClosedOrbitResult.h

Verification plan

1. Pure solver tests: affine maps with known fixed points, coupled stable 4D maps, singular (M-I), damping, non-convergence, and finite-difference accuracy.
2. Geometry tests: reject LINE; accept a closed RING; reject position/orientation mismatch with both errors reported.
3. Physics sanity: ideal symmetric magnetic ring converges; a fresh turn closes; tune agrees with a known linear model.
4. Legacy benchmark: compare a representable cyclotron configuration against old OPAL at the same energy.
5. TRACK consistency: seed one particle with the COF result and advance it through the shared TRACK physics path; closure meets XTOL/PTOL.
6. Perturbations: static misalignment/field errors shift the orbit as expected; lost trials and near-integer tunes fail cleanly.
7. Convergence study: repeat with DT, DT/2, DT/4 and demonstrate expected pusher order. Solver tolerances must be justified against discretization error.

Delivery sequence

1. Land RING Modelling Capabilities #543 RING topology and clone-propagation tests.
2. Extract/test the shared reference-particle one-turn service without changing TRACK results.
3. Add the pure 4D Newton solver and analytic-map tests.
4. Add COF/RUN/ENDCOF, validation, diagnostics, and a magnetic-ring integration test.
5. Add legacy cyclotron and MPI comparisons, then update the user and physics manuals.
6. Consider 6D RF closure only after the 4D interface is stable.

Deferred 6D extension

For RF rings, extend the state with longitudinal phase/time and momentum deviation, define time residual modulo a declared master RF period, include energy gain, and solve a 6D fixed point. Multiple incommensurate RF frequencies, radiation loss, and stochastic effects require explicit semantics and should not be approximated by quietly adding two coordinates.

References

• RING Modelling Capabilities #543, New RING Line
• Legacy OPAL implementation: https://github.com/OPALX-project/OPAL/blob/master/src/Distribution/ClosedOrbitFinder.h
• M. M. Gordon, Computation of closed orbits and basic focusing properties for sector-focused cyclotrons and the design of CYCLOPS (1983), cited by the legacy implementation.