546 compute linear transfer matrix for LINE and RING

[aaadelmann]

Summary and scope

Implements numerical first-order external-field transport maps for the design reference of a LINE or one traversal of a RING, addressing #546 and providing the map infrastructure needed by the future closed-orbit work in #544.

• An OPTION enables map calculation.
• Maps are reconstructed by tracking local (not in PartBunch) perturbed rays through the actual external fields, including supported field-map elements.
• Maps and their reference states are attached to runtime element occurrences. Overlapping fields are handled by unique, path-ordered segments; the complete line/one-turn product is also available and printed to stdout.
• Finite-difference amplitudes and Richardson differentiation order are selectable independently of the integration-method interface.

Note

Some of this documentation will go to the regression tests as soon as we have merged that

Revision status: the implementation and C++ unit tests are included in this PR at b992a3378. The numerical tables below record completed validation runs. Input-file sketches are included directly in this description; the study scripts and output archives are not part of the source diff. Export to the separate regression repository and end-to-end CI integration remain pending.

1. Coordinates and observation planes

The particle momentum used internally is the dimensionless mechanical momentum, not canonical momentum:

$$ \mathbf u=\frac{\mathbf p}{mc}=\boldsymbol\beta\gamma, \qquad \gamma=\sqrt{1+\mathbf u\cdot\mathbf u}. $$

At each reference boundary, let the reference state be $(\mathbf r_0,\mathbf u_0,t_0)$, let $u_0=\lVert\mathbf u_0\rVert$, and define the orthonormal beam frame $(\mathbf e_x,\mathbf e_y,\mathbf e_s)$ with $\mathbf e_s=\mathbf u_0/u_0$. The map coordinates are

$$ \begin{aligned} X&=(x,x',y,y',\zeta,\delta)^T,\\ x&=(\mathbf r-\mathbf r_0)\cdot\mathbf e_x, &x'&=\frac{\mathbf u\cdot\mathbf e_x}{\mathbf u\cdot\mathbf e_s},\\ y&=(\mathbf r-\mathbf r_0)\cdot\mathbf e_y, &y'&=\frac{\mathbf u\cdot\mathbf e_y}{\mathbf u\cdot\mathbf e_s},\\ \zeta&=-\beta_0c(t-t_0), &\delta&=\frac{\lVert\mathbf u\rVert}{u_0}-1, \qquad \beta_0=\frac{u_0}{\sqrt{1+u_0^2}}. \end{aligned} $$

Positions and $\zeta$ are in metres; slopes and $\delta$ are dimensionless. Each ray is observed at its own crossing time of the common reference-normal plane

$$ (\mathbf r-\mathbf r_{0,\mathrm{exit}})\cdot\mathbf e_{s,\mathrm{exit}}=0. $$

The crossing is bracketed and located by repeated numerical integration/bisection, not by sampling all rays at one laboratory time. An entrance perturbation is converted back to a ray by

$$ \begin{aligned} \mathbf r&=\mathbf r_0+x\mathbf e_x+y\mathbf e_y,\\ \mathbf u&=u_0(1+\delta) \frac{\mathbf e_s+x'\mathbf e_x+y'\mathbf e_y}{\sqrt{1+x'^2+y'^2}},\\ t&=t_0-\frac{\zeta}{\beta_0c}. \end{aligned} $$

Minimum-rotation / Bishop frame

Between successive tangents $\mathbf e_s$ and $\mathbf e_s^{\mathrm{new}}$, the transverse axes follow the minimum rotation with

$$ \mathbf a=\frac{\mathbf e_s\times\mathbf e_s^{\mathrm{new}}} {\lVert\mathbf e_s\times\mathbf e_s^{\mathrm{new}}\rVert}, \qquad \theta=\mathbf{atan2}\left( \lVert\mathbf e_s\times\mathbf e_s^{\mathrm{new}}\rVert, \mathbf e_s\cdot\mathbf e_s^{\mathrm{new}}\right). $$

Rodrigues rotation, followed by orthonormalization, transports the frame without adding arbitrary roll about the tangent. Straight sections retain their transverse axes; the degenerate reversal case is handled explicitly. This avoids a Frenet-frame curvature singularity in drifts.

References for the frame construction and rotation:

• R. L. Bishop, There is More than One Way to Frame a Curve, The American Mathematical Monthly 82(3), 246–251 (1975), DOI: 10.2307/2319846. This is the original reference for the continuous relatively parallel (Bishop / rotation-minimizing) frame.
• W. Wang, B. Jüttler, D. Zheng and Y. Liu, Computation of Rotation Minimizing Frames, ACM Transactions on Graphics 27(1), Article 2 (2008), section 2.1, DOI: 10.1145/1330511.1330513. Its discussion of the tangent-aligning rotation method provides the discrete counterpart to the construction above. OPALX uses that minimum-rotation approach with an atan2 angle and explicit degeneracy handling; it does not implement the paper's double-reflection method or an exact finite-step solution of the continuous Bishop equations.
• K. M. Lynch and F. C. Park, Modern Robotics, section 3.2.3: Exponential Coordinates of Rotation, provides a derivation of Rodrigues' axis-angle rotation formula used by rotateRodrigues. The historical source is O. Rodrigues, Des lois géométriques qui régissent les déplacements d'un système solide dans l'espace…, Journal de mathématiques pures et appliquées, series 1, volume 5, 380–440 (1840).

2. Two passes, field overlap and ownership

OrbitThreader owns the reference pass, records reference states and overlap participation (isOverlapping), and builds the field-support information. LinearTransferMapBuilder then refines the boundaries, transports frames, launches rays, reconstructs segment matrices and composes the result. It does not attach results or print diagnostics; those responsibilities remain with OrbitThreader.

The reference path is partitioned wherever its active-element set changes. There is one calculated matrix per unique segment, not one independent whole-element matrix that could double-count an overlap. Each participating runtime element receives a copy of its segment map. An element can therefore own several consecutive maps.

Every private ray independently evaluates the active set at its own position:

$$ \mathbf E(\mathbf r,t)=\sum_{i\in A(\mathbf r)}\mathbf E_i(\mathbf r,t), \qquad \mathbf B(\mathbf r,t)=\sum_{i\in A(\mathbf r)}\mathbf B_i(\mathbf r,t). $$

Reference-segment ownership never freezes a perturbed ray's field selection. Collective/self fields are excluded; RF cavities and traveling-wave structures are currently rejected.

For segment matrices in increasing reference-path order,

$$ M_{\mathrm{total}}=M_N M_{N-1}\cdots M_2M_1. $$

Every overlap appears exactly once in this product. Field-free gaps are included, and the backward bunch pre-roll is excluded from the requested map interval. A ring result is one-turn transport about the supplied reference; this does not itself find or guarantee a closed orbit.

Storage and access

• Structure/LinearTransferMap.h: the $6\times6$ matrix, entrance/exit reference state and frame, segment/pass IDs, active-element names, overlap flag, differentiation settings, condition estimate and diagnostics.
• ElementBase::getLinearTransferMaps(): maps involving that runtime occurrence. Results are cleared for recalculation and are not inherited by cloned elements.
• OpalBeamline::getLinearTransferMapsInReferenceOrder(): traversal of element-owned maps, deduplicating shared overlap segments. Unowned field-free gaps are not in this element traversal.
• OrbitThreader::getCombinedLinearTransferMap(): the complete product, including those gaps. Rank zero prints it and both diagnostics to stdout; --info 2 captures the study output.

3. Matrix reconstruction and Richardson refinement

For each input coordinate $j$, launch two rays with $X_{\mathrm{in}}^{\pm,j}=\pm\epsilon_j\mathbf e_j$. Encode the actual entrance and exit coordinates and form

$$ (D_{\mathrm{in}})_{:j}=\frac{X_{\mathrm{in}}^{+,j}-X_{\mathrm{in}}^{-,j}}{2}, \qquad (D_{\mathrm{out}})_{:j}=\frac{X_{\mathrm{out}}^{+,j}-X_{\mathrm{out}}^{-,j}}{2}. $$

The matrix satisfies

$$ M D_{\mathrm{in}}=D_{\mathrm{out}}, \qquad M=D_{\mathrm{out}}D_{\mathrm{in}}^{-1}. $$

The implementation uses a pivoted $6\times6$ inversion and rejects singular/ill-conditioned input differences. Re-encoding the entrance states accounts for finite-precision coordinate conversion rather than assuming a perfectly diagonal input matrix. Centered differences have second-order error for a smooth map.

For LINEARTRANSFERMAPRICHARDSON=L, with $0\le L\le4$, repeat the reconstruction at amplitudes $\boldsymbol\epsilon/2^k$ and form

$$ \begin{aligned} R_{k,0}&=M(\boldsymbol\epsilon/2^k),\\ R_{k,j}&=R_{k,j-1} +\frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1}, \qquad 1\le j\le k\le L. \end{aligned} $$

The returned segment matrix is $R_{L,L}$. With a smooth even-power error expansion its formal differentiation error is $O(\epsilon^{2(L+1)})$, at a cost of $12(L+1)$ rays per segment. Extrapolation happens before segment composition. The last change in each matrix column is stored as richardsonCorrection; it is a convergence indicator, not an error bound.

Richardson extrapolates perturbation size, not the integration time step. It does not remove integration, field interpolation, event-location or roundoff errors, nor does it guarantee symplecticity. Time-step and amplitude convergence must be checked separately.

Example supported by the current PR head:

OPTION, ENABLELINEARTRANSFERMAPS = TRUE,
        LINEARTRANSFERMAPINTEGRATOR = "BORIS",
        LINEARTRANSFERMAPRICHARDSON = 1,
        LINEARTRANSFERMAPSTEPS = {1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5};

Meaning of LINEARTRANSFERMAPSTEPS

LINEARTRANSFERMAPSTEPS supplies the six starting finite-difference perturbation amplitudes, ordered as $(x,x',y,y',\zeta,\delta)$. They are not integration time steps, particle-distribution widths, or error tolerances. At each segment entrance the builder perturbs one coordinate at a time by plus and minus its amplitude, leaving the other five at their reference values. This gives twelve local rays per amplitude level, not one ray displaced in all six coordinates.

Array entry	Coordinate	Unit	Meaning of 1e-5 in the example
1	$x$	m	Horizontal displacement of ±10 micrometres
2	$x'$	dimensionless slope	Horizontal slope perturbation of ±0.00001 (approximately ±10 microradians near the reference)
3	$y$	m	Vertical displacement of ±10 micrometres
4	$y'$	dimensionless slope	Vertical slope perturbation of ±0.00001 (approximately ±10 microradians near the reference)
5	$\zeta$	m	Longitudinal time-of-flight coordinate perturbation of ±10 micrometres
6	$\delta$	dimensionless relative momentum	Relative momentum-magnitude perturbation of ±0.00001, i.e. ±10 ppm

Here $\zeta=-\beta_0c(t-t_0)$, so the fifth entry shifts the ray's launch time by $\Delta t=-\Delta\zeta/(\beta_0c)$ while keeping it on the entrance plane; it is not a displacement along the beam axis. The sixth entry uses $\delta=\lVert\mathbf u\rVert/u_0-1$, with $\mathbf u=\boldsymbol\beta\gamma$. It is a relative change in momentum magnitude, not an absolute beta-gamma increment or a relative kinetic-energy change.

With LINEARTRANSFERMAPRICHARDSON=1, the example is evaluated at amplitudes 1e-5 and 5e-6 and the centered maps are extrapolated before composition. This uses 24 rays per segment. DT independently controls numerical ray integration. Smaller perturbations reduce finite-difference truncation error but eventually amplify roundoff and integration/event-location errors; the example is not a universal optimum.

The starting defaults are {1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3} in the same mixed units, with Richardson level zero. Exactly six finite positive amplitudes are required; they must remain nonzero after refinement, and the sixth amplitude must be less than one so both perturbed momenta remain forward and nonzero.

Default integration method

Boris is the default: LINEARTRANSFERMAPINTEGRATOR = "BORIS". If map calculation is enabled without selecting a method, the initial setting uses symmetric drift–Boris-kick–drift for both the map-enabled reference and the local perturbed rays. The explicit "BORIS" in the example therefore documents the default; it is not required to select it in a fresh run.

LF2 is an alias for BORIS. RK4 and DOP853 require an explicit selection and apply to map computation only; they do not change the production particle pusher. The DOP853 validation runs in section 9 are explicit test settings, not a change of default. Unsupported method names fail explicitly. Later OPTION statements retain previously selected map settings when those attributes are omitted; omission does not reset an earlier RK4/DOP853 choice to Boris.

4. Ray integration and boundary resolution

ExternalFieldRayTracker is the shared host-side integrator for the map-enabled reference and private rays.
The default integrator uses symmetric drift–Boris-kick–drift.

This is exact for constant velocity up to momentum/roundoff error, and second-order quadrature with electric acceleration. Position, time and path additions retain compensated-summation residuals, including when reference/ray differences are formed. This improves small-increment bookkeeping without changing the production particle pusher or the Boris kick order.

For every ray, compare active supports at the trial start, drift midpoint and end. If they differ, bisect the time interval, advance its first half, then recompute the second from that state. Transition-containing substeps terminate when

$$ c\lvert h\rvert\le \max\left(10^{-12}c\lvert\Delta t\rvert, 64\epsilon_{\mathrm{mach}}\max(1,\mathrm m,\lVert\mathbf r\rVert)\right). $$

Here $h$ is the trial duration and $\Delta t$ the nominal step. A cap based on one quarter of the shortest positive longitudinal field-support extent prevents regular longitudinal crossings from skipping a complete thin element. Element isInside() queries and the bend containment chart are authoritative. These operations are numerical subdivision/reintegration, not exact element maps. Grazing geometry, unresolved transverse features and discontinuities inside field tables still require adequate nominal resolution.

The PR also provides classical RK4 and the 12-stage eighth-order DOP853 formula for map computation only. They evaluate the same equations at their own stages; stage support sets are checked as well. Both use a fixed nominal time step, not an adaptive controller. DOP853's embedded estimators/dense output are not used. Neither these methods nor ray-dependent subdivision is a symplecticity guarantee.

Be aware of the costs: 9 field evaluations for RK4 and 25 for DOP853, versus 1 for Boris. The production particle pusher and secondary threading remain Boris.

5. Determinant and symplecticity diagnostics

For each segment and the complete product, report

$$ \eta_{\det}=\lvert\det M-1\rvert,\qquad \eta_J=\max_{i,j}\lvert(M^TJM-J)_{ij}\rvert, $$

$$ J=\mathbf{diag}(J_2,J_2,J_2),\qquad J_2=\begin{pmatrix}0&1\-1&0\end{pmatrix}. $$

The unchanged stdout diagnostic threshold is $10^{-6}$. Determinant one is necessary but not sufficient for a canonical symplectic map. Expressing mechanical momenta in $\beta\gamma$ does not make them canonical; slopes, momentum normalization and nonzero endpoint fields also matter. Consequently the canonical-$J$ residual is explicitly labelled a diagnostic, not a universal invariant or a substitute for accuracy against an independent reference.

6. DBA validation: hard-edge baseline in this PR

The hard-edge DBA baseline is the symmetric sequence B1–D1–QACH–D2–B2, with

$$ \rho=2,\mathrm m,\quad \theta=\pi/6,\quad L_b=\rho\theta, \quad L_d=1,\mathrm m,\quad L_q=0.2,\mathrm m, $$

and $K_1=-6.371966681365967,\mathrm m^{-2}$ in OPALX's quadrupole convention. HGAP=FINT=0. The thick-lens hard-edge solution has $R_{16}=R_{26}=0$ and is obtained by the ordered product $B D Q D B$, including longitudinal transport, not merely the transverse blocks.

The saved single-rank Boris case at $\Delta t=10^{-13},\mathrm s$, starting $\epsilon_j=10^{-5}$ and one Richardson level gives:

Quantity	Result
Maximum scaled error over all 36 entries	8.36648e-9
$R_{16}$ [m]	-5.338758153958e-9
$R_{26}$	-2.974797943206e-10
$\eta_{\det}$	2.9031169e-9
$\eta_J$	2.9001920e-9

The lattice can be summarized by this input-file sketch. It omits the beam, source, field-solver and tracking blocks; survey coordinates are rounded for readability, so this is not a precision-replay deck. The reference particle and numerical settings must be supplied separately as stated for each benchmark.

REAL LB = 2 * PI / 6;
REAL KACH = -6.371966681365967;
B1: SBEND, L=LB, ANGLE=PI/6, HGAP=0, FINT=0,
    X=0, Y=0, Z=0, THETA=0;
D1: DRIFT, L=1,
    X=-0.267949192431123, Y=0, Z=1, THETA=-PI/6;
QACH: QUADRUPOLE, L=0.2, K1=KACH,
    X=-0.767949192431123, Y=0, Z=1.866025403784439, THETA=-PI/6;
D2: DRIFT, L=1,
    X=-0.867949192431123, Y=0, Z=2.039230484541327, THETA=-PI/6;
B2: SBEND, L=LB, ANGLE=PI/6, HGAP=0, FINT=0,
    X=-1.367949192431123, Y=0, Z=2.905255888325765, THETA=-PI/6;
LATTICE: LINE = (B1,D1,QACH,D2,B2);

With zero gap the FINT value has no effect. The numbered validation cases use the rigidity-matched momentum and explicit DOP853 settings in section 9; they are distinct from the historical Boris run reported here.

7. DBA with finite fringes: wider-spacing benchmark

This retained three-integrator benchmark has wider spacing to separate the fringe supports; it is distinct from the fixed-nominal-geometry in section 9.

The same bend bodies and quadrupole seed are retained, with HGAP=0.1 m, FINT=0.1, full gap $g=0.2,\mathrm m$ and support extension $w=5g=1,\mathrm m$ outside each body face. Explicitly field-free D1/D2 remain 1 m long; supports are separated:

Element	Full support interval [m]	Nominal body interval [m]
B1	0.000000–3.047198	1.000000–2.047198
D1	3.047198–4.047198	same
QACH	4.047198–4.247198	same
D2	4.247198–5.247198	same
B2	5.247198–8.294395	6.247198–7.294395

The total requested reference path is 8.294395102393196 m. The old quadrupole strength is now an unmatched seed: zero dispersion must not be assumed for this larger-spacing, finite-fringe lattice. No dipole-strength renormalization or rematching was performed. All five elements were visited, with no recorded support overlaps and zero logged fields in D1/D2. Field checks use the actual global reference positions, not the nominal design-line coordinate alone.

The native profile is

$$ F(d)=\frac{1}{1+\exp!\left(\sum_{i=0}^{5}c_i(d/g)^i\right)}, \qquad A(z)=\min{F(-z),F(z-L_b)}. $$

The six coefficients are $(0.478959,1.911289,-1.185953,1.630554,-1.082657,0.318111)$. In local bend coordinates, the implemented near-axis law is

$$ B_x=\sigma C A'(z)y,\qquad B_y=B_0\left(A(z)-\tfrac12 A''(z)y^2\right),\qquad B_s=B_0A'(z)y, $$

where $\sigma=+1$ for the entrance and $-1$ for the exit. The full definition of $C$, units, support mask and limitations is in the DBA manual chapter. FINT is invisible on the mid-plane field trace, so off-axis derivatives are included in the reference validation.

Independent reference and full-map errors

The reference integrates the analytic Lorentz variational equations, including analytic field derivatives, crossing-time sensitivity at hard quadrupole faces, and projection to the common exit plane. It uses neither OPALX nor shadow-ray differencing. This is a numerically evaluated distributed-field reference, not a closed-form solution. Its last DOP853 refinement changes the full matrix by 2.0233e-11; an independent RK45 calculation differs by 1.7714e-11. These are convergence checks, not rigorous error bounds.

The actual distribution has $u_0=490.23677597553325$, hence momentum 0.2505104777748827 GeV/c with the current electron mass. The field-normalization input is P0 = 0.2505104781131461 GeV/c. Both are preserved as actually run; their small difference is not silently removed.

Complete independent distributed-field reference matrix (rounded for display)

Coordinates: (x, x', y, y', zeta, delta); positions in metres.

-2.085824445604  -8.697818516437   0               0               0  -1.976820978592
-0.385368056578  -2.086398703073   0               0               0  -0.247747623263
 0               0               6.936040964462  32.072898829273  0   0
 0               0               1.467539017179   6.930211437947  0   0
 0.245045609781   1.969572860924   0               0               1   0.692183661897
 0               0               0               0               0   1

Entries below 1e-14 are shown as zero. Reference determinant and canonical-J residuals are 1.24e-12 and 9.97e-13 before this display rounding.

Settings are the observed best full-map points from the incomplete previous map-2 study, not guaranteed optima for map-3:

Integrator	DT [s]	Starting epsilon	Richardson levels	$R_{16}$ [m]	$R_{26}$
Boris	1e-13	1e-5	1	-1.976821004850	-0.247747632004
RK4	3e-11	3e-5	0	-1.976824228889	-0.247749986316
DOP853	1e-11	0.003	1	-1.976825802611	-0.247748090546

Define $\Delta M=M_{\mathrm{OPALX}}-M_{\mathrm{ref}}$. To compare all 36 entries without mixing units, use $S=\mathbf{diag}(L_,1,L_,1,L_,1)$ with $L_=1,\mathrm m$ and measure

$$ \eta_M=\max_{i,j}\lvert(S^{-1}\Delta M S)_{ij}\rvert. $$

The 1 m value is a coordinate scale, not an element length.

Integrator	$\eta_M$	Largest entry	$\Delta R_{16}$ [m]	$\Delta R_{26}$
Boris	9.64778e-8	$R_{12}$	-2.62585e-8	-8.74138e-9
RK4	4.05225e-5	$R_{34}$	-3.25030e-6	-2.36305e-6
DOP853	9.96550e-6	$R_{12}$	-4.82402e-6	-4.67283e-7

Integrator	Signed $\det M-1$ (approximately)	$\eta_{\det}$ from stdout	$\eta_J$ from stdout	OrbitThreader wall time [s]
Boris	-7.76e-9	7.7564677e-9	9.6620733e-9	41.60
RK4	-1.03e-6	1.0262537e-6	1.0222106e-6	2.898
DOP853	-1.74e-7	1.7433327e-7	5.6822864e-7	19.49

The signed values are recomputed from the saved text-rounded matrices; stdout evaluates the absolute residuals before that rounding. RK4 is slightly outside the existing 1e-6 thresholds in both diagnostics. No tolerances were loosened. These are different parameter choices, not an equal-cost or equal-accuracy integrator ranking; the timing includes reference threading and map construction.

Distinguish numerical error from the CERN thin-edge approximation

Analytic fringe edge maps exist, as derived in the CERN lecture, slides 23–29. In the full-gap convention, for zero face rotation, $\psi=h(2,\mathrm{HGAP})\mathrm{FINT}$. The current native helper instead uses $h,\mathrm{HGAP},\mathrm{FINT}$: 0.010 rad versus 0.005 rad for this case. Its caller does not compensate for the factor of two. This concerns the explicit added edge coefficient, not all focusing from the distributed field, and is not corrected by this PR description.

The fixed Enge profile is also not constrained to have the input FINT as its measured fringe integral. Accordingly a thin-edge model on the nominal design orbit is not identical to the tracked distributed model. The separate CAS thin-edge DBA gives $R_{16}=-1.929713679677,\mathrm m$ and $R_{26}=-0.207853671267$, with a largest map difference of about 0.393 from the three OPALX runs. That is a model/reference-orbit difference, not an integration error. The zero-rotation FINT coefficient affects the vertical block only and cannot explain the horizontal/dispersion discrepancy by itself. Resolving the fringe-model convention is separate follow-up work.

8. Tests and planned regression coverage

The pushed TestOrbitThreader coverage includes map attachment, independent builder ownership, Bishop frames, overlap summation/deduplication, thin-support crossings, reverse integration, one-turn traversal independent of production step count, Richardson order/tableau/settings, compensated increments, magnetic path-length bookkeeping and accelerated path-target localization.

The map-3 analysis passed 25 Python tests, both in the working tree and an isolated export prepared for reproducibility without private campaign directories. These cover exact full-matrix drift/quadrupole/uniform relativistic bend maps, analytic field/Lorentz Jacobians, independently checked quadrupole-face sensitivity, reference cross-integration, archived settings/artifact integrity and thin-edge-model sanity checks. Four shared SDDS parser tests also pass. The rebuilt C++ TestOrbitThreader and TestBendFieldModel targets pass on one rank; all 25 OrbitThreader cases also pass in a dedicated two-rank MPI run. These are performed checks, not a claim that end-to-end study replay is already wired into PR CI.

The DBA cell and these results will become an automated regression test. The planned fixture will retain the input, reference matrix, numerical settings and provenance, and check all 36 entries, especially $R_{16}$/$R_{26}$, both determinant and canonical-J diagnostics, separated supports, and the complete line-map product. The analytic hard-edge baseline and distributed-fringe reference must remain distinct. Any future fringe-physics correction requires an explicitly reviewed reference update rather than freezing the present convention as a universal physics truth.

Regression integration is still pending; this description does not claim it has already been added to CI. Acceptance tolerances will be justified from reference convergence and reproducibility, not chosen to hide the RK4 diagnostic failure. Study runs remain single-rank; dedicated MPI regression checks are a separate concern.

9. Fixed-geometry linear-map validation: seven tests

This is an additional validation set, distinct from the support-separated, wider-spacing map-3 benchmark in section 7. It tests nominal-body map ownership, summed fringe fields inside nominal drifts, and same-section ring return handling from 86830cacf. The source/unit-test diff recorded with these runs matches that commit's changes relative to e15681b58.

Regression status: the new Python harness, numbered inputs, recorded outputs and regression export are prepared locally for the separate regression-test repository; they are not yet published there or wired into CI. The results below are completed local runs. The wider-spacing benchmark above remains unchanged.

Common setup and independent references

All seven cases use one electron, one MPI rank, one OpenMP thread, no collective fields, map-only DOP853, DT=1e-11 s, starting perturbations 1e-4 in each coordinate's units, and LINEARTRANSFERMAPRICHARDSON=1. Production tracking uses MAXSTEPS=1; OrbitThreader still computes the entire requested map. Stdout is retained as lin-map-N.out with --info 2.

The launch momentum is $u_0=\beta\gamma=490.23677597553325$. The field-normalization parameter is set consistently to $p_0=0.2505104777748827~\mathrm{GeV}/c$. This removes the small field/particle rigidity mismatch retained in the historical map-3 fixture; these are new tests, not bitwise replays of that benchmark.

Python calculates the independent combined map before reading OPALX output. Hard-edge references use the full $6\times6$ ordered product $M_N\cdots M_1$, including longitudinal transport and dispersion. The finite-fringe reference is an independently integrated variational ODE for the analytic native field, not a closed-form fringe map and not a fit to OPALX rays. OPALX itself continues to obtain every map numerically from perturbed-ray tracking.

What each test checks

1. lin-map-1 — 1 m drift. Checks the complete relativistic drift map, including $R_{12}=R_{34}=L$ and $R_{56}=L/\gamma_0^2$, with zero transverse dispersion. This is the simplest check of coordinate encoding, observation-plane crossing and longitudinal time-of-flight sensitivity.

2. lin-map-2 — single thick quadrupole. A 0.4 m quadrupole with $K_1=1~\mathrm{m}^{-2}$. For this electron and the native sign convention it defocuses horizontally and focuses vertically. The reference uses the exact hyperbolic/trigonometric thick-lens blocks and the longitudinal drift term. It checks focusing signs and finite-length transport, not a thin-lens approximation.

3. lin-map-3 — single sector bend. A 30-degree bend with $\rho=2~\mathrm m$, $L_b=\pi/3~\mathrm m$ and HGAP=0. Checks the full sector-bend map, including $R_{16}=\rho(1-\cos\theta)$, $R_{26}=\sin\theta$, the transverse-to-longitudinal entries, and the relativistic $R_{56}$ term. These dispersion entries are intentionally nonzero.

4. lin-map-4 — thick-lens FODO. The sequence is D(0.4 m)–QF(0.2 m)–D(0.8 m)–QD(0.2 m)–D(0.4 m), with $K_{1F}=-1~\mathrm{m}^{-2}$ and $K_{1D}=+1~\mathrm{m}^{-2}$. Checks ordered composition across multiple hard boundaries against the analytic product. The full matrix is checked; $R_{16}=R_{26}=0$ because there is no bending.

5. lin-map-5 — hard-edge DBA. B1(30 degrees, radius 2 m)–D1(1 m)–Q(0.2 m)–D2(1 m)–B2(30 degrees, radius 2 m), with $K_1=-6.371966681365967~\mathrm{m}^{-2}$ and HGAP=0. Its nominal length is 4.294395102393195 m. Checks the complete thick-lens DBA product and dispersion cancellation: analytic $R_{16}$ and $R_{26}$ are zero to roundoff. FINT has no effect in this zero-gap limit.

6. lin-map-6 — finite-fringe DBA at fixed nominal geometry. Retains exactly the bodies, placements, drift lengths and quadrupole strength of lin-map-5, but uses HGAP=0.1 m, FINT=0.1. Bend tails may occupy nominal drifts; they do not add lattice length. The independent reference sums all field sources, integrates analytic field sensitivities, applies hard-quadrupole crossing sensitivity, and projects onto the exit observation plane. The launch is at the first body entrance, where the field is nonzero; the exit plane is normal to the reference at a travelled path of 4.294395102393195 m. There is no support padding or rematching, and zero dispersion is not expected. This case currently fails the map-accuracy criterion.

7. lin-map-7 — six-cell hard-edge DBA ring. Six copies of lin-map-5, represented by 30 uniquely named elements in a RING, give 360 degrees of nominal bending and a design circumference of 25.766370614359165 m. The independent map is $M_{\mathrm{DBA}}^6$. Checks same-section return handling and the complete one-turn composition; it does not invoke a closed-orbit finder. The reference returns within $7.3905\times10^{-13}~\mathrm m$ of its launch position, with relative momentum mismatch $5.6603\times10^{-13}$. This seed lattice is optically unstable; it is a composition test, not a proposed matched operating ring.

Input-file sketches for the numbered cases

Each fragment below belongs to a separate input file; it is not a complete standalone deck. Common beam/source/solver blocks are omitted. The straight examples use the global +Z direction. For both DBA cases use the five-body placement sketch in section 6, with the changes described below. No deprecated longitudinal placement attribute is needed.

// lin-map-1: drift
D: DRIFT, L=1, X=0, Y=0, Z=0;
LATTICE: LINE = (D);

// lin-map-2: thick quadrupole
Q: QUADRUPOLE, L=0.4, K1=1, X=0, Y=0, Z=0;
LATTICE: LINE = (Q);

// lin-map-3: single sector bend
B: SBEND, L=2*PI/6, ANGLE=PI/6, HGAP=0, X=0, Y=0, Z=0;
LATTICE: LINE = (B);

// lin-map-4: FODO
D1: DRIFT,      L=0.4,        X=0, Y=0, Z=0;
QF: QUADRUPOLE, L=0.2, K1=-1, X=0, Y=0, Z=0.4;
D2: DRIFT,      L=0.8,        X=0, Y=0, Z=0.6;
QD: QUADRUPOLE, L=0.2, K1=1,  X=0, Y=0, Z=1.4;
D3: DRIFT,      L=0.4,        X=0, Y=0, Z=1.6;
LATTICE: LINE = (D1,QF,D2,QD,D3);

For lin-map-5, retain the section-6 DBA body positions and use HGAP=0, FINT=0.1. For lin-map-6, change only the two bend field definitions to HGAP=0.1, FINT=0.1, with HAPERT=0.5; do not move the bodies or increase the drift lengths. The actual numbered fixtures set APERTURE="RECTANGLE(1.0,0.2)" on the drift and quadrupole elements. The nominal drifts can contain the summed fringe fields. This differs from the wider-spacing benchmark in section 7.

For lin-map-7, define and explicitly survey six successive copies of the hard-edge cell, with unique names for all 30 occurrences, and use this sequence declaration. C1_* through C6_* are separate placed elements, not six copies left at the same global coordinates.

LATTICE: RING = (
  C1_B1,C1_D1,C1_Q,C1_D2,C1_B2,
  C2_B1,C2_D1,C2_Q,C2_D2,C2_B2,
  C3_B1,C3_D1,C3_Q,C3_D2,C3_B2,
  C4_B1,C4_D1,C4_Q,C4_D2,C4_B2,
  C5_B1,C5_D1,C5_Q,C5_D2,C5_B2,
  C6_B1,C6_D1,C6_Q,C6_D2,C6_B2);

The shared numerical settings for reproducing the section-9 table are the following explicit overrides; Boris remains the default when no method has been selected.

OPTION, ENABLELINEARTRANSFERMAPS=TRUE,
        LINEARTRANSFERMAPINTEGRATOR="DOP853",
        LINEARTRANSFERMAPRICHARDSON=1,
        LINEARTRANSFERMAPSTEPS={1e-4,1e-4,1e-4,1e-4,1e-4,1e-4};
REAL P0 = 0.2505104777748827;
// Supply a one-electron source with u0=490.23677597553325 and no collective fields.
// TRACK uses MAXSTEPS=1, DT=1e-11 and the case's stated nominal path target.
// Run on one MPI rank with one OpenMP thread and save stdout using --info 2.

Manual-ready comparison table

Let $\Delta M=M_{\mathrm{OPALX}}-M_{\mathrm{ref}}$. Use the coordinate scale $L_{\mathrm{scale}}=1~\mathrm m$ and define

$$ \begin{aligned} S &= \mathrm{diag}(L_{\mathrm{scale}},1,L_{\mathrm{scale}},1,L_{\mathrm{scale}},1),\\ \widehat M &= S^{-1}MS. \end{aligned} $$

$$ \begin{aligned} \eta_M &= \max_{i,j}\left|\widehat M_{ij}-\widehat M^{\mathrm{ref}}_{ij}\right|,\\ q_{\mathrm{tol}} &= \max_{i,j} \frac{\left|\widehat M_{ij}-\widehat M^{\mathrm{ref}}_{ij}\right|} {a+r\left|\widehat M^{\mathrm{ref}}_{ij}\right|}. \end{aligned} $$

The 1 m scale makes the mixed-unit map comparison dimensionless; it is not an element length. PASS means $q_{\mathrm{tol}}\le1$ across all 36 entries, including expected zeros. The fixed absolute/relative tolerances $(a,r)$ are (1e-9, 1e-9) for lin-map-1, (2e-7, 2e-8) for lin-map-2 through lin-map-5, and (5e-7, 5e-8) for lin-map-6 and lin-map-7. They were selected before the runs and were not relaxed after a failure. The relative term accommodates large ring-map entries; these bounds are not uncertainty estimates or cross-platform guarantees.

Test	Model / reference	Max map error	ΔR16 [m]	ΔR26	Tolerance ratio	Time [s]	Result
lin-map-1	1 m drift (analytic)	5.062e-12	0.000e+00	0.000e+00	5.062e-03	2.39	PASS
lin-map-2	0.4 m quadrupole (analytic)	8.886e-12	0.000e+00	0.000e+00	4.443e-05	2.78	PASS
lin-map-3	30 degree sector bend (analytic)	5.922e-11	-5.423e-12	-2.790e-11	2.884e-04	2.68	PASS
lin-map-4	Thick-lens FODO (analytic)	1.010e-11	0.000e+00	0.000e+00	5.051e-05	16.72	PASS
lin-map-5	Hard-edge DBA (analytic)	4.755e-11	-9.164e-12	3.534e-11	2.378e-04	18.39	PASS
lin-map-6	DBA with fringes in nominal drifts (ODE)	1.898e-05	7.118e-06	2.563e-06	3.359e+01	18.71	FAIL
lin-map-7	Six-cell hard-edge DBA ring (analytic)	2.817e-08	1.514e-09	1.113e-09	8.241e-03	309.87	PASS

Time is elapsed wall time of the complete OPALX process, not just an integrator kernel; the Python reference calculation is excluded. These are different lattices at common numerical settings, not an integrator-efficiency ranking.

Known discrepancy and diagnostic caveats

• The finite-fringe reference changes by at most 6.3221e-12 under the recorded DOP853 refinement/RK45 cross-check and also passes an independent event-located nonlinear-ray check. OPALX's vertical block agrees within 3.4e-12, but horizontal/longitudinal entries differ. Reducing epsilon from 1e-4 to 1e-5 and 1e-6 gives maximum errors 2.1770e-5 and 2.2052e-5; halving DT to 5e-12 s at epsilon 1e-6 gives 2.2044e-5. The cause remains unresolved. Both calculations use the same native FINT convention; this failure is distinct from the separate CAS/MAD-X convention question.
• PASS above is the entrywise accuracy test, not a claim that every printed invariant diagnostic passes. The ring reports abs(det(M)-1)=3.4098e-7 and a canonical-J residual of 1.6625e-6; the latter exceeds the unchanged 1e-6 stdout threshold. Large vertical amplification makes these subtractions roundoff-sensitive. Both diagnostics remain available in the recorded CSV/stdout.

Regression-suite preparation

The local checkout prepared for regression-tests-x contains six passing fixtures under RegressionTests/lin-map-N; the finite-fringe lin-map-6 is retained under disabledTests with its strict independent reference and unchanged tolerance. The checker never substitutes an observed OPALX matrix for the independent golden matrix and never converts a known failure into success.

The standard-library runner checks all 36 entries, input/reference hashes, missing/duplicate/nonfinite maps, process failures and timeouts. It preserves stdout and returns nonzero on a mismatch. Twenty local Python verification tests pass, including corrupted-entry tests and replay of the exported independent reference generator; the packaged .local launcher also passed a fresh single-rank drift run.

These fixtures will be used for the OPALX regression suite after review and push; CI wiring is still pending. The existing .rt particle-statistics checks do not validate transfer matrices, so CI must invoke the dedicated full-map runner. The known-failing fringe case must be resolved before activation. The existing manual links below provide the implementation background; the table and descriptions in this section are prepared for later inclusion in the manual.

Manual

• OPALX manual
• Numerical linear transfer maps — implementation mathematics and limitations
• map-3 fringe DBA — field model, both references, results and convention caveats

The chapter links point to the verified current manual source because the corresponding new published HTML chapter URLs were not yet available when this description was prepared. The implementation is included in this PR; the numerical studies are documented here through input sketches, reference definitions and results. The separate regression-suite integration status is stated explicitly above.

[aaadelmann]

sandbox will not be merged in master !