546 compute linear transfer matrix for LINE and RING

.gitignore

@@ -42,4 +42,7 @@ input-files/timing.dat
 build_*
 *.json
 coverage_html/
-.idea
+.idea
+
+# Local study inputs, outputs and exploratory work are not part of the source tree.
+/sandbox/

Doxyfile.in

@@ -1453,14 +1453,16 @@ MATHJAX_FORMAT         = HTML-CSS
 # The default value is: http://cdn.mathjax.org/mathjax/latest.
 # This tag requires that the tag USE_MATHJAX is set to YES.
 
-MATHJAX_RELPATH        = http://cdn.mathjax.org/mathjax/latest
+# The former cdn.mathjax.org endpoint is retired. Pin the MathJax 2 release
+# used by this HTML-CSS configuration rather than following a major-version update.
+MATHJAX_RELPATH        = https://cdn.jsdelivr.net/npm/mathjax@2.7.9
 
 # The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax
 # extension names that should be enabled during MathJax rendering. For example
 # MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols
 # This tag requires that the tag USE_MATHJAX is set to YES.
 
-MATHJAX_EXTENSIONS     = AMSmath AMSsymbols
+MATHJAX_EXTENSIONS     = TeX/AMSmath TeX/AMSsymbols
 
 # The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces
 # of code that will be used on startup of the MathJax code. See the MathJax site

src/AbsBeamline/ElementBase.cpp

@@ -24,6 +24,15 @@
 
 extern Inform* gmsg;
 
+bool ElementBase::isInsideBody(const Vector_t<double, 3>& r) const {
+    const auto& geometry = getGeometry();
+    const auto local = geometry.kind() == GeometryKind::SBend
+            ? GeometryHelper::toBendArcCoords(r, geometry.getCurvature(), geometry.getElementLength())
+            : r;
+    return local(2) >= 0.0 && local(2) < geometry.getElementLength()
+            && ApertureHelper::isInsideAperture(local, aperture_m);
+}
+
 const std::vector<double> ElementBase::defaultAperture_m = std::vector<double>({1e6, 1e6});
 
 const std::map<ElementType, std::string> ElementBase::elementTypeToString_s = {
@@ -40,7 +49,6 @@ const std::map<ElementType, std::string> ElementBase::elementTypeToString_s = {
         {ElementType::SBEND, "SBEND"},
         {ElementType::RBEND, "RBEND"},
         {ElementType::RBEND3D, "RBEND3D"},
-        {ElementType::RING, "Ring"},
         {ElementType::SOURCE, "SOURCE"},
         {ElementType::SOLENOID, "SOLENOID"},
         {ElementType::PROBE, "Probe"},
@@ -60,6 +68,7 @@ ElementBase::ElementBase(const ElementBase& right)
       RefPartBunch_m(nullptr),
       online_m(right.online_m),
       elementID(right.elementID),
+      beamlineMembership_m(right.beamlineMembership_m),
       userAttribs(right.userAttribs),
       positionIsFixed(right.positionIsFixed),
       elementPosition_m(right.elementPosition_m),
@@ -75,6 +84,7 @@ ElementBase::ElementBase(const std::string& name)
       RefPartBunch_m(nullptr),
       online_m(false),
       elementID(name),
+      beamlineMembership_m(),
       userAttribs(),
       positionIsFixed(false),
       elementPosition_m(0.0),
@@ -89,6 +99,49 @@ const std::string& ElementBase::getName() const { return elementID; }
 
 void ElementBase::setName(const std::string& name) { elementID = name; }
 
+const BeamlineMembership& ElementBase::getBeamlineMembership() const {
+    return beamlineMembership_m;
+}
+
+BeamlineTopology ElementBase::getBeamlineTopology() const { return beamlineMembership_m.topology; }
+
+const std::string& ElementBase::getBeamlineOwnerName() const {
+    return beamlineMembership_m.ownerName;
+}
+
+bool ElementBase::hasLinearTransferMaps() const { return !linearTransferMaps_m.empty(); }
+
+const std::vector<LinearTransferMap>& ElementBase::getLinearTransferMaps() const {
+    return linearTransferMaps_m;
+}
+
+void ElementBase::addLinearTransferMap(LinearTransferMap map) {
+    linearTransferMaps_m.push_back(std::move(map));
+}
+
+void ElementBase::clearLinearTransferMaps() { linearTransferMaps_m.clear(); }
+
+bool ElementBase::isOverlapping() const { return isOverlapping_m; }
+
+void ElementBase::setOverlapping(const bool overlapping) { isOverlapping_m = overlapping; }
+
+void ElementBase::setBeamlineMembership(BeamlineTopology topology, std::string ownerName) {
+    if (topology == BeamlineTopology::RING && ownerName.empty()) {
+        throw GeneralOpalException(
+                "ElementBase::setBeamlineMembership()",
+                "RING membership requires a non-empty owner name");
+    }
+    if (topology == BeamlineTopology::LINE && !ownerName.empty()) {
+        throw GeneralOpalException(
+                "ElementBase::setBeamlineMembership()",
+                "LINE membership cannot have an owner name");
+    }
+
+    beamlineMembership_m = {topology, std::move(ownerName)};
+}
+
+void ElementBase::clearBeamlineMembership() { beamlineMembership_m = {}; }
+
 void ElementBase::setOutputFN(const std::string fn) { outputfn_m = fn; }
 
 std::string ElementBase::getOutputFN() const {

src/AbsBeamline/ElementBase.h

@@ -17,6 +17,7 @@
 
 #include "AbsBeamline/AttributeSet.h"
 #include "Algorithms/CoordinateSystemTrafo.h"
+#include "Structure/LinearTransferMap.h"
 #include "Algorithms/Quaternion.hpp"
 #include "BeamlineGeometry/Geometry.h"
 #include "OPALTypes.h"
@@ -49,7 +50,6 @@ enum class ElementType : unsigned short {
     SBEND,
     RBEND,
     RBEND3D,
-    RING,
     PROBE,
     VACUUM,
     SOLENOID,
@@ -58,6 +58,37 @@ enum class ElementType : unsigned short {
     CONSTANTFOCUSING
 };
 
+/// Logical sequence topology, independent of ElementType and field-support overlap.
+/// RING is a sequence declaration, not a new physical element type or a closure test.
+enum class BeamlineTopology : unsigned short { LINE, RING };
+
+/**
+ * @brief Logical RING membership carried by a sequence member or runtime occurrence.
+ *
+ * An occurrence is one use of an element in a sequence, rather than its reusable input
+ * definition. Ring::prepareForTracking() clones occurrences recursively and assigns the
+ * enclosing RING name to each, including members inside nested LINEs. For example, two
+ * uses of a quadrupole in ring R1 have distinct element objects but both carry
+ * {BeamlineTopology::RING, "R1"}; neither the quadrupole name nor an inner LINE name
+ * is the ownerName. The original quadrupole definition is not retagged.
+ *
+ * Ordinary LINE membership is {BeamlineTopology::LINE, ""}: LINE names and a full
+ * parent hierarchy are deliberately not recorded. ElementBase's copy constructor
+ * preserves this metadata, including when OpalBeamline clones tracking elements.
+ * This is host-side value metadata, not a pointer to a sequence, C++ memory ownership,
+ * a unique occurrence identifier, a turn counter, or evidence of orbit closure.
+ *
+ * Do not use membership to select fields or assign transfer-map intervals. Those use,
+ * respectively, ElementBase::isInside() and ElementBase::isInsideBody(). Field tails
+ * can act in another element's nominal body without changing either element's membership.
+ */
+struct BeamlineMembership {
+    /// LINE by default; RING after assignment by the enclosing Ring.
+    BeamlineTopology topology = BeamlineTopology::LINE;
+    /// Enclosing RING's name (required for RING), empty for LINE; stored by value.
+    std::string ownerName;
+};
+
 enum class ApertureType : unsigned short { RECTANGULAR, ELLIPTICAL };
 
 /**
@@ -184,6 +215,70 @@ class ElementBase : public std::enable_shared_from_this<ElementBase> {
     /// @return The element type string.
     static std::string getTypeString(ElementType type);
 
+    /* ======================= Beamline membership =========================== */
+
+    /// Read the logical sequence tag; the returned reference is owned by this element.
+    /// See BeamlineMembership for its distinction from map ownership and field support.
+    const BeamlineMembership& getBeamlineMembership() const;
+
+    /// Return the topology of the owning beam sequence.
+    BeamlineTopology getBeamlineTopology() const;
+
+    /// Return the owning RING name, or an empty string for LINE membership.
+    const std::string& getBeamlineOwnerName() const;
+
+    /// @brief Whether this runtime occurrence owns at least one calculated linear transfer map.
+    bool hasLinearTransferMaps() const;
+
+    /**
+     * @brief Maps attached to this nominal body by the design-orbit threader.
+     *
+     * OrbitThreader appends in reference-path order. One element may receive several
+     * segments, and a segment with several nominal owners is copied to each owner.
+     * Field contribution alone does not imply ownership: a magnet's fringe can act
+     * in a drift-owned map. These are runtime results, not maps of isolated fields.
+     *
+     * The vector belongs to this element. References/pointers to its entries can be
+     * invalidated by addLinearTransferMap(), clearLinearTransferMaps(), or a new
+     * map-enabled design pass. Clones start with an empty vector.
+     * @see OpalBeamline::getLinearTransferMapsInReferenceOrder()
+     */
+    const std::vector<LinearTransferMap>& getLinearTransferMaps() const;
+
+    /// Append a map by value; does not validate, sort, deduplicate, or update overlap flags.
+    void addLinearTransferMap(LinearTransferMap map);
+
+    /// Remove attached maps only; membership and isOverlapping() are unchanged.
+    void clearLinearTransferMaps();
+
+    /**
+     * @brief Whether the map-enabled design pass observed shared field support.
+     *
+     * True if this occurrence participated in a sampled support set containing more
+     * than one element along the requested reference interval. This is an accumulated
+     * flag for that pass, not an overlap partner list or a full 3D intersection test.
+     * Support selection can include a field-free drift; it does not test whether each
+     * selected element contributes a nonzero field. Thus true need not mean overlapping
+     * nominal bodies or an invalid lattice.
+     *
+     * Initially false, not inherited by clones, and reset before map-enabled design
+     * threading. Secondary-species threading leaves the design result unchanged.
+     * @see LinearTransferMap::includesOverlappingFields
+     */
+    bool isOverlapping() const;
+
+    /// @brief Set overlap participation discovered by the design-orbit threader.
+    void setOverlapping(bool overlapping);
+
+    /// Assign the logical sequence tag only; does not change placement, maps or overlap flags.
+    /// @param topology LINE or RING.
+    /// @param ownerName Enclosing RING name, or empty for LINE.
+    /// @throws GeneralOpalException If RING has an empty name or LINE a nonempty name.
+    void setBeamlineMembership(BeamlineTopology topology, std::string ownerName = {});
+
+    /// Restore {BeamlineTopology::LINE, ""} only; maps and overlap flags are unchanged.
+    void clearBeamlineMembership();
+
     /// @brief Apply a visitor.
     /// @note This method must be overridden by derived classes. It should call
     ///       the method of the visitor corresponding to the element class. If
@@ -207,6 +302,19 @@ class ElementBase : public std::enable_shared_from_this<ElementBase> {
     /// @return True if r is inside the field interval.
     virtual bool isInside(const Vector_t<double, 3>& r) const;
 
+    /**
+     * @brief Nominal body containment, independent of longitudinal field support.
+     *
+     * The local longitudinal interval is [0,L), with arc coordinates for a sector bend.
+     * The transverse aperture is respected. A drift can own this interval even when a
+     * neighbouring magnet contributes a fringe field there. Zero-length bodies own no interval.
+     * This defines transfer-map ownership, not the region used for field evaluation.
+     * @param r Position [m] in the element entrance frame, not the lab frame.
+     * @return Whether the point belongs to the nominal body and transverse aperture.
+     * @see OpalBeamline::getBodyElements()
+     */
+    bool isInsideBody(const Vector_t<double, 3>& r) const;
+
     /// @brief Get the bounding box.
     /// @return The bounding box in lab coordinates.
     virtual BoundingBox getBoundingBoxInLabCoords() const;
@@ -425,6 +533,15 @@ class ElementBase : public std::enable_shared_from_this<ElementBase> {
     std::string elementID;
     static const std::map<ElementType, std::string> elementTypeToString_s;
 
+    // --- Beamline membership ---
+    /// Logical sequence tag, preserved by the copy constructor; see BeamlineMembership.
+    BeamlineMembership beamlineMembership_m;
+
+    /// Host-owned map copies. Runtime results below are deliberately not copied.
+    std::vector<LinearTransferMap> linearTransferMaps_m;
+    /// Reference-pass support-overlap diagnostic, not static lattice membership.
+    bool isOverlapping_m{false};
+
     // --- User-defined attributes ---
     AttributeSet userAttribs;
 

src/AbstractObjects/BeamSequence.h

@@ -47,6 +47,10 @@ class BeamSequence : public Element {
     //  The result it the ideal line.
     virtual Beamline* fetchLine() const = 0;
 
+    /// Prepare occurrence-specific state before TrackRun constructs the tracking lattice.
+    /// The default is a no-op (LINE); Ring overrides it to clone/tag member occurrences.
+    virtual void prepareForTracking() {}
+
 protected:
     /// Constructor for exemplars.
     BeamSequence(int size, const char* name, const char* help);

src/Algorithms/CMakeLists.txt

@@ -6,6 +6,8 @@ set (_SRCS
     DistributionMoments.cpp
     Flagger.cpp
     IndexMap.cpp
+    ExternalFieldRayTracker.cpp
+    LinearTransferMapBuilder.cpp
     OrbitThreader.cpp
     ParallelTracker.cpp
     PartData.cpp
@@ -22,10 +24,15 @@ set (HDRS
     AbstractTimeDependence.h
     AbstractTracker.h
     CavityAutophaser.h
+    CompensatedSum.h
+    RungeKuttaTableau.h
     DefaultVisitor.h
     DistributionMoments.h
     Flagger.h
     IndexMap.h
+    ExternalFieldRayTracker.h
+    LinearTransferMapBuilder.h
+    LinearTransferMap.h
     OrbitThreader.h
     ParallelTracker.h
     PartData.h

src/Algorithms/CompensatedSum.h

@@ -0,0 +1,24 @@
+// Copyright (c) 2026, Paul Scherrer Institute, Villigen PSI, Switzerland
+#ifndef OPAL_COMPENSATED_SUM_H
+#define OPAL_COMPENSATED_SUM_H
+
+namespace compensated {
+    /**
+     * @brief Kahan addition retaining increments below the current sum's spacing.
+     * The represented value is sum - correction. Do not reassociate these operations
+     * (e.g. with fast-math). Both values must be retained when copying a tracked state.
+     */
+    inline void add(double increment, double& sum, double& correction) {
+        const double adjusted = increment - correction;
+        const double next = sum + adjusted;
+        correction = (next - sum) - adjusted;
+        sum = next;
+    }
+
+    /// Difference of two compensated values, preserving their small residuals.
+    inline double difference(double left, double leftCorrection, double right,
+                             double rightCorrection) {
+        return (left - right) - (leftCorrection - rightCorrection);
+    }
+}
+#endif