make sure that the initial position is inside the bounding box

src/Algorithms/IndexMap.cpp

+//
+// Class IndexMap
+//
+// This class stores and prints the sequence of elements that the referenc particle passes.
+// Each time the reference particle enters or leaves an element an entry is added to the map.
+// With help of this map one can determine which element can be found at a given position.
+//
+// Copyright (c) 2016,       Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
+//               2017 - 2020 Christof Metzger-Kraus
+//
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
 #include <map>
 #include <limits>
 #include <iostream>
@@ -128,8 +150,7 @@ void IndexMap::tidyUp(double zstop) {
 
     if (rit != mapRange2Element_m.rend() &&
         (*rit).second.size() == 0 &&
-        zstop > (*rit).first.first &&
-        zstop < (*rit).first.second) {
+        zstop > (*rit).first.first) {
 
         key_t key((*rit).first.first, zstop);
         value_t val;

src/Algorithms/IndexMap.h

+//
+// Class IndexMap
+//
+// This class stores and prints the sequence of elements that the referenc particle passes.
+// Each time the reference particle enters or leaves an element an entry is added to the map.
+// With help of this map one can determine which element can be found at a given position.
+//
+// Copyright (c) 2016,       Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
+//               2017 - 2020 Christof Metzger-Kraus
+//
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
 #ifndef OPAL_INDEXMAP_H
 #define OPAL_INDEXMAP_H
 

src/Algorithms/OrbitThreader.cpp

@@ -112,6 +112,7 @@ void OrbitThreader::execute() {
     std::set<std::string> visitedElements;
 
     trackBack();
+    updateBoundingBoxWithCurrentPosition();
 
     Vector_t nextR = r_m / (Physics::c * dt_m);
     integrator_m.push(nextR, p_m, dt_m);
@@ -414,6 +415,7 @@ void OrbitThreader::computeBoundingBox() {
     FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
     FieldList::iterator it = allElements.begin();
     const FieldList::iterator end = allElements.end();
+
     globalBoundingBox_m.lowerLeftCorner = Vector_t(std::numeric_limits<double>::max());
     globalBoundingBox_m.upperRightCorner = Vector_t(-std::numeric_limits<double>::max());
 
@@ -424,6 +426,15 @@ void OrbitThreader::computeBoundingBox() {
         ElementBase::BoundingBox other = it->getBoundingBoxInLabCoords();
         globalBoundingBox_m.getCombinedBoundingBox(other);
     }
+
+    updateBoundingBoxWithCurrentPosition();
+}
+
+
+void OrbitThreader::updateBoundingBoxWithCurrentPosition() {
+    Vector_t dR = Physics::c * dt_m * p_m / Util::getGamma(p_m);
+    ElementBase::BoundingBox pos = ElementBase::BoundingBox::getBoundingBox({r_m - 10 * dR, r_m + 10 * dR});
+    globalBoundingBox_m.getCombinedBoundingBox(pos);
 }
 
 double OrbitThreader::computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,

src/Algorithms/OrbitThreader.h

@@ -108,7 +108,7 @@ private:
     void processElementRegister();
     void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements);
     void computeBoundingBox();
-    // double computeMaximalImplicitDrift();
+    void updateBoundingBoxWithCurrentPosition();
     double computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
                                            const Vector_t & position,
                                            const Vector_t & direction) const;

src/Classic/AbsBeamline/ElementBase.cpp

@@ -309,6 +309,17 @@ bool ElementBase::isInsideTransverse(const Vector_t &r) const
     }
 }
 
+ElementBase::BoundingBox ElementBase::BoundingBox::getBoundingBox(const std::vector<Vector_t> & points) {
+    const Vector_t & point = points.front();
+    BoundingBox result{point, point};
+    for (const Vector_t & point: points) {
+        BoundingBox tmp{point, point};
+        result.getCombinedBoundingBox(tmp);
+    }
+
+    return result;
+}
+
 bool ElementBase::BoundingBox::isInside(const Vector_t & position) const {
     Vector_t relativePosition = position - lowerLeftCorner;
     Vector_t diagonal = upperRightCorner - lowerLeftCorner;
@@ -392,4 +403,26 @@ ElementBase::BoundingBox ElementBase::getBoundingBoxInLabCoords() const {
     }
 
     return bb;
+}
+
+void ElementBase::BoundingBox::print(std::ostream & out) const {
+    const Vector_t & ll = lowerLeftCorner;
+    const Vector_t & ur = upperRightCorner;
+    Vector_t diagonal = ur - ll;
+    Vector_t dX(diagonal(0), 0, 0), dY(0, diagonal(1), 0), dZ(0, 0, diagonal(2));
+
+    std::vector<Vector_t> corners{ll, ll + dX, ll + dX + dY, ll + dY, ur, ur - dX, ur - dX - dY, ur - dY};
+    std::vector<std::vector<unsigned int>> paths{{0, 1, 2, 3}, {0, 1, 7, 6}, {1, 2, 4, 7}, {2, 3, 5, 4}, {3, 0, 6, 5}, {4, 5, 6, 7}};
+
+    out << std::setprecision(8);
+    for (const std::vector<unsigned int> path: paths) {
+        for (unsigned int i : {0, 1, 2, 3, 0}) {
+            const Vector_t & corner = corners[path[i]];
+            out << std::setw(16) << corner(0)
+                << std::setw(16) << corner(1)
+                << std::setw(16) << corner(2)
+                << std::endl;
+        }
+        out << std::endl;
+    }
 }
\ No newline at end of file

src/Classic/AbsBeamline/ElementBase.h

@@ -345,6 +345,8 @@ public:
         Vector_t lowerLeftCorner;
         Vector_t upperRightCorner;
 
+        static BoundingBox getBoundingBox(const std::vector<Vector_t> & points);
+
         void getCombinedBoundingBox(const BoundingBox & other) {
             for (unsigned int d = 0; d < 3; ++ d) {
                 lowerLeftCorner[d] = std::min(lowerLeftCorner[d], other.lowerLeftCorner[d]);
@@ -354,6 +356,8 @@ public:
 
         bool isInside(const Vector_t &) const;
 
+        void print(std::ostream &) const;
+
         /*! Computes the intersection point between a bounding box and the ray which
          *  has the direction 'direction' and starts at the position 'position'. If
          *  the position is inside the box then the algorithm should find an inter-
@@ -609,4 +613,4 @@ inline
 bool ElementBase::isElementPositionSet() const
 { return elemedgeSet_m; }
 
-#endif // CLASSIC_ElementBase_HH
+#endif // CLASSIC_ElementBase_HH
\ No newline at end of file