current status

src/Classic/AbsBeamline/Bend2D.cpp

@@ -1680,4 +1680,37 @@ std::array<double,2> Bend2D::getExitFringeFieldLength() const {
     extFFL[1] = ( exitParameter2_m-exitParameter1_m ); //before edge
 
     return extFFL;
-}
\ No newline at end of file
+}
+
+BoundaryBox Bend2D::getBoundaryBoxInLabCoords() const {
+    CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
+    CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
+
+    double &x = aperture_m.second[0];
+    double &y = aperture_m.second[1];
+
+    std::vector<Vector_t> corners(8);
+    for (int i = -1; i < 2; i += 2) {
+        for (int j = -1; j < 2; j += 2) {
+            unsigned int idx = (i + 1)/2 + (j + 1);
+            corners[idx] = toBegin.transformFrom(Vector_t(i * x, j * y, 0.0));
+            corners[idx + 4] = toEnd.transformFrom(Vector_t(i * x, j * y, 0.0));
+        }
+    }
+
+    BoundingBox bb;
+    bb.lowerLeftCorner = bb.upperRightCorner = corners[0];
+
+    for (unsigned int i = 1; i < 8u; ++ i) {
+        for (unsigned int d; d < 3u; ++ d) {
+            if (bb.lowerLeftCorner(d) > corners[i](d)) {
+                bb.lowerLeftCorner(d) = corners[i](d);
+            }
+            if (bb.upperLeftCorner(d) < corners[i](d)) {
+                bb.upperLeftCorner(d) = corners[i](d);
+            }
+        }
+    }
+
+    return bb;
+}

src/Classic/AbsBeamline/Bend2D.h

@@ -129,6 +129,7 @@ public:
    //  Used to create fringe fields in ThickTracker, (before edge[m], after edge[m])
    std::array<double,2> getExitFringeFieldLength() const;
 
+    BoundaryBox getBoundaryBoxInLabCoords() const override;
 protected:
     void setMessageHeader(const std::string & header);
     double getStartField() const;
@@ -375,5 +376,4 @@ Vector_t Bend2D::transformToExitRegion(const Vector_t &R) const {
     return toExitRegion_m.transformTo(R);
 }
 
-#endif // CLASSIC_BEND_H
-
+#endif // CLASSIC_BEND_H
\ No newline at end of file

src/Classic/AbsBeamline/BendBase.h

@@ -37,7 +37,6 @@ public:
 
     void setFieldMapFN(std::string fileName);
     std::string getFieldMapFN() const;
-
 protected:
     /// Calculate design radius from design energy and field amplitude
     double calcDesignRadius(double fieldAmplitude) const;

src/Classic/AbsBeamline/ElementBase.cpp

@@ -304,4 +304,64 @@ void ElementBase::setCurrentSCoordinate(double s) {
             elementEdge_m = actionRange_m.front().first;
         }
     }
+}
+
+bool ElementBase::isInsideTransverse(const Vector_t &r) const
+{
+    double x = aperture_m.second[0];
+    double y = aperture_m.second[1];
+    double f = 1.0;
+    if (aperture_m.first == CONIC_RECTANGULAR ||
+        aperture_m.first == CONIC_ELLIPTIC) {
+        Vector_t rRelativeToBegin = getEdgeToBegin().transformTo(r);
+        double fractionLength = rRelativeToBegin(2) / getElementLength();
+        f = fractionLength * aperture_m.second[2];
+    }
+
+    switch(aperture_m.first) {
+    case RECTANGULAR:
+        return (std::abs(r[0]) < aperture_m.second[0] && std::abs(r[1]) < aperture_m.second[1]);
+    case ELLIPTICAL:
+        return (std::pow(r[0] / aperture_m.second[0], 2) + std::pow(r[1] / aperture_m.second[1], 2) < 1.0);
+    case CONIC_RECTANGULAR:
+        return (std::abs(r[0]) < f * aperture_m.second[0] && std::abs(r[1]) < f * aperture_m.second[1]);
+    case CONIC_ELLIPTICAL:
+        return (std::pow(r[0] / (f * aperture_m.second[0]), 2) + std::pow(r[1] / (f * aperture_m.second[1]), 2) < 1.0);
+    default:
+        return false;
+    }
+}
+
+BoundaryBox ElementBase::getBoundaryBoxInLabCoords() const {
+    CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
+    CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
+
+    double &x = aperture_m.second[0];
+    double &y = aperture_m.second[1];
+    double &f = aperture_m.second[2];
+
+    std::vector<Vector_t> corners(8);
+    for (int i = -1; i < 2; i += 2) {
+        for (int j = -1; j < 2; j += 2) {
+            unsigned int idx = (i + 1)/2 + (j + 1);
+            corners[idx] = toBegin.transformFrom(Vector_t(i * x, j * y, 0.0));
+            corners[idx + 4] = toEnd.transformFrom(Vector_t(i * f * x, j * f * y, 0.0));
+        }
+    }
+
+    BoundingBox bb;
+    bb.lowerLeftCorner = bb.upperRightCorner = corners[0];
+
+    for (unsigned int i = 1; i < 8u; ++ i) {
+        for (unsigned int d; d < 3u; ++ d) {
+            if (bb.lowerLeftCorner(d) > corners[i](d)) {
+                bb.lowerLeftCorner(d) = corners[i](d);
+            }
+            if (bb.upperLeftCorner(d) < corners[i](d)) {
+                bb.upperLeftCorner(d) = corners[i](d);
+            }
+        }
+    }
+
+    return bb;
 }
\ No newline at end of file

src/Classic/AbsBeamline/ElementBase.h

@@ -359,6 +359,8 @@ public:
     void setRotationAboutZ(double rotation);
     double getRotationAboutZ() const;
 
+    virtual BoundaryBox getBoundaryBoxInLabCoords() const;
+
 protected:
     bool isInsideTransverse(const Vector_t &r, double f = 1) const;
 
@@ -531,24 +533,7 @@ inline
 bool ElementBase::isInside(const Vector_t &r) const
 {
     const double length = getElementLength();
-    return isInsideTransverse(r, r(2) / length * aperture_m.second[2]) && r(2) >= 0.0 && r(2) < length;
-}
-
-inline
-bool ElementBase::isInsideTransverse(const Vector_t &r, double f) const
-{
-    switch(aperture_m.first) {
-    case RECTANGULAR:
-        return (std::abs(r[0]) < aperture_m.second[0] && std::abs(r[1]) < aperture_m.second[1]);
-    case ELLIPTICAL:
-        return (std::pow(r[0] / aperture_m.second[0], 2) + std::pow(r[1] / aperture_m.second[1], 2) < 1.0);
-    case CONIC_RECTANGULAR:
-        return (std::abs(r[0]) < f * aperture_m.second[0] && std::abs(r[1]) < f * aperture_m.second[1]);
-    case CONIC_ELLIPTICAL:
-        return (std::pow(r[0] / (f * aperture_m.second[0]), 2) + std::pow(r[1] / (f * aperture_m.second[1]), 2) < 1.0);
-    default:
-        return false;
-    }
+    return r(2) >= 0.0 && r(2) < length && isInsideTransverse(r);
 }
 
 inline

src/Elements/OpalRBend.cpp

@@ -189,7 +189,11 @@ void OpalRBend::update() {
 
     if(itsAttr[HAPERT]) {
         double hapert = Attributes::getReal(itsAttr[HAPERT]);
-        bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({hapert, hapert, 1.0}));
+        double gap = Attributes::getReal(itsAttr[GAP]);
+        bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({hapert, gap, 1.0}));
+    } else {
+        double gap = Attributes::getReal(itsAttr[GAP]);
+        bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({0.5, gap, 1.0});
     }
 
     if(itsAttr[LENGTH])
@@ -216,4 +220,4 @@ void OpalRBend::update() {
 
     // Transmit "unknown" attributes.
     OpalElement::updateUnknown(bend);
-}
+}
\ No newline at end of file