diff --git a/.gitignore b/.gitignore index 379271838..2e3035a75 100644 --- a/.gitignore +++ b/.gitignore @@ -42,4 +42,7 @@ input-files/timing.dat build_* *.json coverage_html/ -.idea \ No newline at end of file +.idea + +# Local study inputs, outputs and exploratory work are not part of the source tree. +/sandbox/ diff --git a/Doxyfile.in b/Doxyfile.in index 54d19f34c..9ba3f240d 100644 --- a/Doxyfile.in +++ b/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 diff --git a/src/AbsBeamline/ElementBase.cpp b/src/AbsBeamline/ElementBase.cpp index 05c565970..28dc2b1fc 100644 --- a/src/AbsBeamline/ElementBase.cpp +++ b/src/AbsBeamline/ElementBase.cpp @@ -24,6 +24,15 @@ extern Inform* gmsg; +bool ElementBase::isInsideBody(const Vector_t& 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 ElementBase::defaultAperture_m = std::vector({1e6, 1e6}); const std::map ElementBase::elementTypeToString_s = { @@ -40,7 +49,6 @@ const std::map 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& 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 { diff --git a/src/AbsBeamline/ElementBase.h b/src/AbsBeamline/ElementBase.h index 7bfe95acc..8d0b86600 100644 --- a/src/AbsBeamline/ElementBase.h +++ b/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 { /// @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& 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 { /// @return True if r is inside the field interval. virtual bool isInside(const Vector_t& 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& 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 { std::string elementID; static const std::map 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 linearTransferMaps_m; + /// Reference-pass support-overlap diagnostic, not static lattice membership. + bool isOverlapping_m{false}; + // --- User-defined attributes --- AttributeSet userAttribs; diff --git a/src/AbstractObjects/BeamSequence.h b/src/AbstractObjects/BeamSequence.h index cdec479ef..8eafb0e8f 100644 --- a/src/AbstractObjects/BeamSequence.h +++ b/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); diff --git a/src/Algorithms/CMakeLists.txt b/src/Algorithms/CMakeLists.txt index 4eb8af802..13ca8ebab 100644 --- a/src/Algorithms/CMakeLists.txt +++ b/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 diff --git a/src/Algorithms/CompensatedSum.h b/src/Algorithms/CompensatedSum.h new file mode 100644 index 000000000..427962f51 --- /dev/null +++ b/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 diff --git a/src/Algorithms/ExternalFieldRayTracker.cpp b/src/Algorithms/ExternalFieldRayTracker.cpp new file mode 100644 index 000000000..e02ad2136 --- /dev/null +++ b/src/Algorithms/ExternalFieldRayTracker.cpp @@ -0,0 +1,248 @@ +// Copyright (c) 2026, Paul Scherrer Institute, Villigen PSI, Switzerland +#include "Algorithms/ExternalFieldRayTracker.h" +#include "Algorithms/CompensatedSum.h" +#include "Algorithms/RungeKuttaTableau.h" +#include +#include +#include "Algorithms/PartData.h" +#include "Elements/OpalBeamline.h" +#include "Utilities/OpalException.h" + +namespace { + using State = ExternalFieldRayTracker::State; + using Fields = ExternalFieldRayTracker::FieldEvaluator; + + // Each k holds h*f for (r[3], u[3], s), avoiding sums of large SI derivatives. + template + State rkAdvance(const State& initial, double h, const Fields& fields, + const PartData& reference, const external_field_rk::Tableau& table, + bool& hitMaterial) { + std::array, Stages> k{}; + const auto combine = [&](const auto& weights, unsigned count, double fraction) { + State ray = initial; + std::array increment{}, correction{}; + for (unsigned j = 0; j < count; ++j) + for (unsigned d = 0; d < 7; ++d) + compensated::add(weights[j] * k[j][d], increment[d], correction[d]); + for (unsigned d = 0; d < 3; ++d) { + compensated::add(increment[d] - correction[d], ray.position(d), + ray.positionCorrection(d)); + ray.momentum(d) += increment[d + 3] - correction[d + 3]; + } + compensated::add(increment[6] - correction[6], ray.pathLength, ray.pathLengthCorrection); + compensated::add(fraction * h, ray.time, ray.timeCorrection); + return ray; + }; + const double force = h * reference.getQ() * Physics::c / reference.getM(); + for (unsigned i = 0; i < Stages; ++i) { + const State ray = combine(table.a[i], i, table.c[i]); + Vector_t electric(0.0), magnetic(0.0); + if (fields(ray, electric, magnetic)) { + hitMaterial = true; + return ray; + } + const double gamma = std::sqrt(1.0 + dot(ray.momentum, ray.momentum)); + const Vector_t velocity = (Physics::c / gamma) * ray.momentum; + const Vector_t du = force * (electric + cross(velocity, magnetic)); + for (unsigned d = 0; d < 3; ++d) { + k[i][d] = h * velocity(d); + k[i][d + 3] = du(d); + } + k[i][6] = h * euclidean_norm(velocity); + } + return combine(table.b, Stages, 1.0); + } +} + +ExternalFieldRayTracker::IntegrationMethod ExternalFieldRayTracker::parseIntegrationMethod( + const std::string& name) { + if (name == "BORIS" || name == "LF2") return IntegrationMethod::BORIS; + if (name == "RK4") return IntegrationMethod::RK4; + if (name == "DOP853") return IntegrationMethod::DOP853; + throw OpalException( + "ExternalFieldRayTracker::parseIntegrationMethod", + "Unsupported external-field ray integrator '" + name + "'. Use BORIS (alias LF2), RK4 or DOP853."); +} + +std::string ExternalFieldRayTracker::integrationMethodName(const IntegrationMethod method) { + switch (method) { + case IntegrationMethod::BORIS: + return "BORIS"; + case IntegrationMethod::RK4: + return "RK4"; + case IntegrationMethod::DOP853: + return "DOP853"; + } + throw OpalException( + "ExternalFieldRayTracker::integrationMethodName", "Unsupported ray integrator enum."); +} + +ExternalFieldRayTracker::ExternalFieldRayTracker( + OpalBeamline& beamline, const PartData& reference, const IntegrationMethod method) + : beamline_m(beamline), reference_m(reference), integrationMethod_m(method) { + integrationMethodName(method); // Validate even when no ray will be advanced. + for (const auto& element : beamline_m.getElements()) { + if (element->getType() == ElementType::MARKER || element->getType() == ElementType::MONITOR) + continue; + double begin = 0.0, end = 0.0; + element->getFieldExtent(begin, end); + const double length = std::abs(end - begin); + if (length > 0.0) maximumStep_m = std::min(maximumStep_m, length / (4.0 * Physics::c)); + // Reference samples must also resolve nominal map-ownership intervals when a field + // extends far beyond a short body. This cap does not truncate the field itself. + const double bodyLength = element->getGeometry().getArcLength(); + if (bodyLength > 0.0) + maximumStep_m = std::min(maximumStep_m, bodyLength / (4.0 * Physics::c)); + } +} + +ExternalFieldRayTracker::Step ExternalFieldRayTracker::step( + const State& initial, const double dt, const FieldEvaluator& fields) const { + if (!std::isfinite(dt)) + throw OpalException("ExternalFieldRayTracker::step", "The ray time step must be finite."); + switch (integrationMethod_m) { + case IntegrationMethod::BORIS: + return borisStep(initial, dt, fields); + case IntegrationMethod::RK4: + case IntegrationMethod::DOP853: + return rungeKuttaStep(initial, dt, fields); + } + throw OpalException("ExternalFieldRayTracker::step", "Unsupported ray integrator enum."); +} + +ExternalFieldRayTracker::Step ExternalFieldRayTracker::rungeKuttaStep( + const State& initial, const double dt, const FieldEvaluator& fields) const { + Step result; + result.midpoint = result.end = initial; + result.duration = dt; + if (dt == 0.0) return result; + const auto integrate = [&](double h) { + if (integrationMethod_m == IntegrationMethod::RK4) + return rkAdvance(initial, h, fields, reference_m, external_field_rk::rk4, result.hitMaterial); + return rkAdvance(initial, h, fields, reference_m, external_field_rk::dop853, result.hitMaterial); + }; + result.end = integrate(dt); + if (result.hitMaterial) return result; + // An independent half-duration solve supplies a same-order diagnostic midpoint. + // It is not composed into the endpoint: nominal DT retains its usual meaning. + result.midpoint = integrate(0.5 * dt); + if (result.hitMaterial) return result; + result.hitMaterial = fields(result.midpoint, result.electric, result.magnetic); + return result; +} + +ExternalFieldRayTracker::Step ExternalFieldRayTracker::borisStep( + const State& initial, const double dt, const FieldEvaluator& fields) const { + Step result; + result.midpoint = result.end = initial; + result.duration = dt; + if (dt == 0.0) return result; + const auto halfDrift = [dt](State& ray) { + const double momentum2 = dot(ray.momentum, ray.momentum); + const double factor = (0.5 * Physics::c * dt) / std::sqrt(1.0 + momentum2); + for (unsigned component = 0; component < 3; ++component) + compensated::add(factor * ray.momentum(component), ray.position(component), + ray.positionCorrection(component)); + compensated::add(factor * std::sqrt(momentum2), ray.pathLength, + ray.pathLengthCorrection); + compensated::add(0.5 * dt, ray.time, ray.timeCorrection); + }; + auto& ray = result.midpoint; + halfDrift(ray); + result.hitMaterial = fields(ray, result.electric, result.magnetic); + result.end = ray; + if (result.hitMaterial) return result; + // BorisPusher::kick does not use its position argument. Keep the shared kick; + // only the host drift/bookkeeping changes, not TRACK's particle kernel. + integrator_m.kick( + ray.position, result.end.momentum, result.electric, result.magnetic, dt, reference_m.getM(), + reference_m.getQ()); + halfDrift(result.end); + return result; +} + +ExternalFieldRayTracker::State ExternalFieldRayTracker::advanceToPathLength( + const State& initial, const double dt, const double target) const { + const auto distance = [target](const State& ray) { + return compensated::difference(ray.pathLength, ray.pathLengthCorrection, target, 0.0); + }; + if (!std::isfinite(target)) + throw OpalException("ExternalFieldRayTracker::advanceToPathLength", "Non-finite path target."); + State trial = advance(initial, dt); + const double direction = std::copysign(1.0, dt); + if (direction * distance(initial) > 0.0 || direction * distance(trial) < 0.0) + throw OpalException("ExternalFieldRayTracker::advanceToPathLength", "Path target is not bracketed."); + if (distance(initial) == 0.0) return initial; + if (distance(trial) == 0.0) return trial; + const double timeTolerance = std::max( + 1.e-12 * std::abs(dt), 64.0 * std::numeric_limits::epsilon() + * std::max(1.0, euclidean_norm(initial.position)) / Physics::c); + double lower = 0.0, upper = dt; + for (unsigned iteration = 0; iteration < 64; ++iteration) { + const double middle = 0.5 * (lower + upper); + trial = advance(initial, middle); + if (distance(trial) == 0.0) break; + if (direction * distance(trial) >= 0.0) upper = middle; + else lower = middle; + if (std::abs(upper - lower) <= timeTolerance) break; + } + return trial; +} + +ExternalFieldRayTracker::State ExternalFieldRayTracker::advance( + const State& initial, const double dt, std::vector* accepted) const { + if (!std::isfinite(dt)) { + throw OpalException( + "ExternalFieldRayTracker::advance", "The ray time step must be finite."); + } + if (dt == 0.0) return initial; + return advanceRecursive(initial, dt, 1.0e-12 * std::abs(dt), accepted, 0); +} + +ExternalFieldRayTracker::State ExternalFieldRayTracker::advanceRecursive( + const State& initial, const double dt, const double tolerance, std::vector* accepted, + const unsigned depth) const { + const double timeFloor = 64.0 * std::numeric_limits::epsilon() + * std::max(1.0, euclidean_norm(initial.position)) / Physics::c; + const bool resolved = std::abs(dt) <= std::max(tolerance, timeFloor); + const auto split = [&]() { + if (depth >= 64) { + throw OpalException( + "ExternalFieldRayTracker::advance", + "Field-boundary refinement did not converge."); + } + const auto middle = advanceRecursive(initial, 0.5 * dt, tolerance, accepted, depth + 1); + return advanceRecursive(middle, 0.5 * dt, tolerance, accepted, depth + 1); + }; + if (!resolved && std::abs(dt) > maximumStep_m) return split(); + + const auto initialSet = beamline_m.getElements(initial.position); + bool crossedSupport = false; + const auto trial = step(initial, dt, [&](const State& ray, auto& electric, auto& magnetic) { + const auto stageSet = beamline_m.getElements(ray.position); + crossedSupport = crossedSupport || stageSet != initialSet; + for (const auto& element : stageSet) { + if (element->getType() == ElementType::MARKER + || element->getType() == ElementType::MONITOR) + continue; + const auto localR = beamline_m.transformToLocalCS(element, ray.position); + const auto localP = beamline_m.rotateToLocalCS(element, ray.momentum); + Vector_t localE(0.0), localB(0.0); + if (element->applyToReferenceParticle(localR, localP, ray.time, localE, localB)) + return true; + electric += beamline_m.rotateFromLocalCS(element, localE); + magnetic += beamline_m.rotateFromLocalCS(element, localB); + } + return false; + }); + if (!resolved + && (crossedSupport || initialSet != beamline_m.getElements(trial.end.position))) + return split(); + if (trial.hitMaterial) { + throw OpalException( + "ExternalFieldRayTracker::advance", + "A ray reached material while evaluating external fields."); + } + if (accepted) accepted->push_back(trial); + return trial.end; +} diff --git a/src/Algorithms/ExternalFieldRayTracker.h b/src/Algorithms/ExternalFieldRayTracker.h new file mode 100644 index 000000000..09be7b3a1 --- /dev/null +++ b/src/Algorithms/ExternalFieldRayTracker.h @@ -0,0 +1,146 @@ +// Copyright (c) 2026, Paul Scherrer Institute, Villigen PSI, Switzerland +#ifndef OPAL_EXTERNAL_FIELD_RAY_TRACKER_H +#define OPAL_EXTERNAL_FIELD_RAY_TRACKER_H + +#include +#include +#include +#include +#include "OPALTypes.h" +#include "Steppers/BorisPusher.h" + +class OpalBeamline; +class PartData; + +/** + * @brief Host-side external-field integration shared by reference and transfer-map rays. + * + * Positions are in metres, time and step size in seconds, and mechanical momentum is + * \f$\mathbf u=\mathbf p/(mc)=\vec\beta\gamma\f$. Fields are in V/m and tesla; + * PartData supplies rest energy in eV and charge in elementary-charge units. + * The default is drift--Boris-kick--drift. Map computation can instead select + * classical RK4 or the twelve-stage eighth-order DOP853 formula. No collective + * fields or analytic element transfer matrices enter this calculation. + * + * Boris drifts are evaluated directly in metres: + * \f$\mathbf r_{1/2}=\mathbf r_0+(ch/2)\mathbf u_0/\gamma_0\f$ and + * \f$\mathbf r_1=\mathbf r_{1/2}+(ch/2)\mathbf u_1/\gamma_1\f$, + * where \f$\gamma_i=\sqrt{1+|\mathbf u_i|^2}\f$. The signed path is advanced by + * \f[ + * \Delta s=\frac{ch}{2}\left(\frac{|\mathbf u_0|}{\gamma_0} + * +\frac{|\mathbf u_1|}{\gamma_1}\right). + * \f] + * This integrates speed, not the chord between endpoints. It is exact for constant + * speed (up to momentum/roundoff error), and second-order quadrature with electric + * acceleration. Position, time and path additions retain Kahan correction terms. + * The Boris kick, its order and field-support tolerances are unchanged. Only this + * host-side reference/map tracker is affected, not the production particle pusher. + * + * RK4 and DOP853 solve the same relativistic Lorentz equations, including path: + * \f[ + * \dot{\mathbf r}=c\mathbf u/\gamma,\qquad + * \dot{\mathbf u}=\frac{qc}{\mathcal E_0} + * \left(\mathbf E+\frac{c}{\gamma}\mathbf u\times\mathbf B\right),\qquad + * \dot s=c|\mathbf u|/\gamma, + * \f] + * where \f$\mathcal E_0=mc^2\f$ is the numerical rest energy in eV and q is the + * signed charge in elementary-charge units. Each stage samples its own position, + * momentum and time. For \f$\mathbf y=(\mathbf r,\mathbf u,s)\f$, + * \f[ + * k_i=h f(t_0+c_i h,\mathbf y_0+\sum_{j