diff --git a/src/Algorithms/IndexMap.cpp b/src/Algorithms/IndexMap.cpp
index f1bb31b09ec4d4c80b92246e442c9cea0e87db23..330b424c93a2c0054880268a2c8c5e56d9ca8e1a 100644
--- a/src/Algorithms/IndexMap.cpp
+++ b/src/Algorithms/IndexMap.cpp
@@ -6,7 +6,7 @@
 #include "Algorithms/IndexMap.h"
 #include "AbstractObjects/OpalData.h"
 #include "AbsBeamline/Multipole.h"
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "Utilities/Util.h"
 #include "Physics/Physics.h"
 #include "OPALconfig.h"
@@ -24,6 +24,8 @@ IndexMap::IndexMap():
 { }
 
 void IndexMap::print(std::ostream &out) const {
+    if (mapRange2Element_m.size() == 0) return;
+
     out << "Size of map " << mapRange2Element_m.size() << " sections " << std::endl;
     out << std::fixed << std::setprecision(6);
     auto mapIti = mapRange2Element_m.begin();
@@ -252,7 +254,7 @@ void IndexMap::saveSDDS(double startS) const {
             case ElementBase::RBEND:
             case ElementBase::SBEND:
                 {
-                    const Bend* bend = static_cast<const Bend*>(element.get());
+                    const Bend2D* bend = static_cast<const Bend2D*>(element.get());
                     if (bend->getRotationAboutZ() > 0.5 * Physics::pi &&
                         bend->getRotationAboutZ() < 1.5 * Physics::pi) {
                         items[DIPOLE] = -1;
diff --git a/src/Algorithms/OrbitThreader.cpp b/src/Algorithms/OrbitThreader.cpp
index f2ea1d95dc41ceec409588df28cdcfaba4ec800e..a66bbb56b4a358df0446d1bc5ee17c9cb73097ab 100644
--- a/src/Algorithms/OrbitThreader.cpp
+++ b/src/Algorithms/OrbitThreader.cpp
@@ -2,6 +2,7 @@
 #include "Algorithms/CavityAutophaser.h"
 
 #include "AbsBeamline/RFCavity.h"
+#include "AbsBeamline/BendBase.h"
 #include "AbsBeamline/TravelingWave.h"
 
 #include "AbstractObjects/OpalData.h"
@@ -80,12 +81,17 @@ OrbitThreader::OrbitThreader(const PartData &ref,
 
 void OrbitThreader::execute() {
     double initialPathLength = pathLength_m;
-    trackBack();
+    double maxDistance = computeMaximalImplicitDrift();
+
+    auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
+    std::set<std::string> visitedElements;
+
+    trackBack(2 * maxDistance);
+
     Vector_t nextR = r_m / (Physics::c * dt_m);
     integrator_m.push(nextR, p_m, dt_m);
     nextR *= Physics::c * dt_m;
-    auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
-    std::set<std::string> visitedElements;
+
     setDesignEnergy(allElements, visitedElements);
 
     auto elementSet = itsOpalBeamline_m.getElements(nextR);
@@ -98,7 +104,7 @@ void OrbitThreader::execute() {
         double initialS = pathLength_m;
         Vector_t initialR = r_m;
         Vector_t initialP = p_m;
-        integrate(elementSet, maxIntegSteps_m);
+        integrate(elementSet, maxIntegSteps_m, 2 * maxDistance);
 
         registerElement(elementSet, initialS,  initialR, initialP);
 
@@ -123,19 +129,21 @@ void OrbitThreader::execute() {
 
             elementSet = itsOpalBeamline_m.getElements(nextR);
         }
-    } while ((time_m - initialTime_m) / dt_m < maxIntegSteps_m &&
-             errorFlag_m != HITMATERIAL &&
+    } while (errorFlag_m != HITMATERIAL &&
              errorFlag_m != EOL);
 
-    // imap_m.tidyUp();
+    imap_m.tidyUp();
     *gmsg << level1 << "\n" << imap_m << endl;
     imap_m.saveSDDS(initialPathLength);
 
     processElementRegister();
 }
 
-void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxSteps) {
+void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxSteps, double maxDrift) {
     static size_t step = 0;
+    CoordinateSystemTrafo labToBeamline = itsOpalBeamline_m.getCSTrafoLab2Local();
+    const double oldPathLength = pathLength_m;
+    double stepLength;
     Vector_t nextR;
     do {
         errorFlag_m = EVERYTHINGFINE;
@@ -165,7 +173,10 @@ void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxStep
             Bf += itsOpalBeamline_m.rotateFromLocalCS(*it, localB);
         }
 
-        if (step % loggingFrequency_m == 0 && Ippl::myNode() == 0 && !OpalData::getInstance()->isOptimizerRun()) {
+        if (pathLength_m > 0.0 &&
+            pathLength_m < zstop_m &&
+            step % loggingFrequency_m == 0 && Ippl::myNode() == 0 &&
+            !OpalData::getInstance()->isOptimizerRun()) {
             logger_m << std::setw(18) << std::setprecision(8) << pathLength_m + std::copysign(euclidean_norm(r_m - oldR), dt_m)
                      << std::setw(18) << std::setprecision(8) << r_m(0)
                      << std::setw(18) << std::setprecision(8) << r_m(1)
@@ -198,12 +209,15 @@ void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxStep
         integrator_m.push(nextR, p_m, dt_m);
         nextR *= Physics::c * dt_m;
 
-        if (pathLength_m > zstop_m) {
+        if ((activeSet.size() == 0 && std::abs(pathLength_m - oldPathLength) > maxDrift) ||
+            (activeSet.size() > 0  && pathLength_m > zstop_m)) {
             errorFlag_m = EOL;
             return;
         }
 
-    } while (((time_m - initialTime_m) / dt_m < maxSteps) &&
+        stepLength = std::abs(dt_m) * euclidean_norm(p_m) / sqrt(dot(p_m, p_m) + 1) * Physics::c;
+    } while ((dt_m > 0.0 ||
+              euclidean_norm(labToBeamline.transformTo(r_m)) > stepLength)  &&
              (activeSet == itsOpalBeamline_m.getElements(nextR)));
 }
 
@@ -263,7 +277,7 @@ double OrbitThreader::getMaxDesignEnergy(const IndexMap::value_t &elementSet) co
     return designEnergy;
 }
 
-void OrbitThreader::trackBack() {
+void OrbitThreader::trackBack(double maxDrift) {
     dt_m *= -1;
     double initialPathLength = pathLength_m;
 
@@ -271,10 +285,11 @@ void OrbitThreader::trackBack() {
     integrator_m.push(nextR, p_m, dt_m);
     nextR *= Physics::c * dt_m;
 
+    maxDrift = std::min(maxDrift, dZ_m);
     while (initialPathLength - pathLength_m < dZ_m) {
         auto elementSet = itsOpalBeamline_m.getElements(nextR);
 
-        integrate(elementSet, 1000);
+        integrate(elementSet, 1000, maxDrift);
 
         nextR = r_m / (Physics::c * dt_m);
         integrator_m.push(nextR, p_m, dt_m);
@@ -368,3 +383,59 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
         }
     }
 }
+
+double OrbitThreader::computeMaximalImplicitDrift() {
+    FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
+    double maxDrift = 0.0;
+    FieldList::iterator it = allElements.begin();
+    const FieldList::iterator end = allElements.end();
+
+    for (; it != end; ++ it) {
+        auto element1 = it->getElement();
+        const auto &toEdge = element1->getCSTrafoGlobal2Local();
+        auto toEnd = element1->getEdgeToEnd() * toEdge;
+        Vector_t end1 = toEnd.transformFrom(Vector_t(0, 0, 0));
+        Vector_t directionEnd = toEnd.rotateFrom(Vector_t(0, 0, 1));
+        if (element1->getType() == ElementBase::RBEND ||
+            element1->getType() == ElementBase::SBEND ||
+            element1->getType() == ElementBase::RBEND3D ) {
+            auto toBegin = element1->getEdgeToBegin() * toEdge;
+
+            BendBase *bend = static_cast<BendBase*>(element1.get());
+            double angleRelativeToFace = bend->getEntranceAngle() - bend->getBendAngle();
+            directionEnd = toBegin.rotateFrom(Vector_t(sin(angleRelativeToFace), 0, cos(angleRelativeToFace)));
+        }
+
+        double minDistanceLocal = std::numeric_limits<double>::max();
+        FieldList::iterator it2 = allElements.begin();
+        for (; it2 != end; ++ it2) {
+            if (it == it2) continue;
+
+            auto element2 = it2->getElement();
+            const auto &toEdge = element2->getCSTrafoGlobal2Local();
+            auto toBegin = element2->getEdgeToBegin() * toEdge;
+            Vector_t begin2 = toBegin.transformFrom(Vector_t(0, 0, 0));
+            Vector_t directionBegin = toBegin.rotateFrom(Vector_t(0, 0, 1));
+            if (element2->getType() == ElementBase::RBEND ||
+                element2->getType() == ElementBase::SBEND ||
+                element2->getType() == ElementBase::RBEND3D ) {
+                BendBase *bend = static_cast<BendBase*>(element2.get());
+                double E1 = bend->getEntranceAngle();
+                directionBegin = toBegin.rotateFrom(Vector_t(sin(E1), 0, cos(E1)));
+            }
+
+
+            double distance = euclidean_norm(begin2 - end1);
+            double directionProjection = dot(directionEnd, directionBegin);
+            if (directionProjection > 0.999 && minDistanceLocal > distance) {
+                minDistanceLocal = distance;
+            }
+        }
+
+        if (maxDrift < minDistanceLocal) {
+            maxDrift = minDistanceLocal;
+        }
+    }
+
+    return maxDrift;
+}
\ No newline at end of file
diff --git a/src/Algorithms/OrbitThreader.h b/src/Algorithms/OrbitThreader.h
index b9f51da5c412aea9fda250768473fb3e996ecbbe..6a78739ede24c9646cf5cc41ef7d5ff58fe06ab0 100644
--- a/src/Algorithms/OrbitThreader.h
+++ b/src/Algorithms/OrbitThreader.h
@@ -69,8 +69,8 @@ private:
 
     std::multimap<std::string, elementPosition> elementRegistry_m;
 
-    void trackBack();
-    void integrate(const IndexMap::value_t &activeSet, size_t maxSteps);
+    void trackBack(double maxDrift);
+    void integrate(const IndexMap::value_t &activeSet, size_t maxSteps, double maxDrift = 10.0);
     bool containsCavity(const IndexMap::value_t &activeSet);
     void autophaseCavities(const IndexMap::value_t &activeSet, const std::set<std::string> &visitedElements);
     double getMaxDesignEnergy(const IndexMap::value_t &elementSet) const;
@@ -78,6 +78,7 @@ private:
     void registerElement(const IndexMap::value_t &elementSet, double, const Vector_t &r, const Vector_t &p);
     void processElementRegister();
     void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements);
+    double computeMaximalImplicitDrift();
 };
 
 inline
@@ -97,4 +98,4 @@ IndexMap::value_t OrbitThreader::getTouchingElements(const std::pair<double, dou
     return imap_m.getTouchingElements(range);
 }
 
-#endif
+#endif
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/Bend.cpp b/src/Classic/AbsBeamline/Bend2D.cpp
similarity index 96%
rename from src/Classic/AbsBeamline/Bend.cpp
rename to src/Classic/AbsBeamline/Bend2D.cpp
index db872e392d038e8c1a25bbc771f5387ab8b6fe2e..eb85aa78ac1b6017745c8f111e327eae384b146d 100644
--- a/src/Classic/AbsBeamline/Bend.cpp
+++ b/src/Classic/AbsBeamline/Bend2D.cpp
@@ -1,12 +1,12 @@
 // ------------------------------------------------------------------------
-// $RCSfile: Bend.cpp,v $
+// $RCSfile: Bend2D.cpp,v $
 // ------------------------------------------------------------------------
 // $Revision: 1.1.1.1 $
 // ------------------------------------------------------------------------
 // Copyright: see Copyright.readme
 // ------------------------------------------------------------------------
 //
-// Definitions for class: Bend
+// Definitions for class: Bend2D
 //   Defines the abstract interface for a sector bend magnet.
 //
 // ------------------------------------------------------------------------
@@ -18,7 +18,7 @@
 //
 // ------------------------------------------------------------------------
 
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "Algorithms/PartBunchBase.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "Utilities/Options.h"
@@ -34,10 +34,10 @@
 
 extern Inform *gmsg;
 
-// Class Bend
+// Class Bend2D
 // ------------------------------------------------------------------------
 
-Bend::Bend():
+Bend2D::Bend2D():
     BendBase(),
     messageHeader_m(" * "),
     pusher_m(),
@@ -78,7 +78,7 @@ Bend::Bend():
 
 }
 
-Bend::Bend(const Bend &right):
+Bend2D::Bend2D(const Bend2D &right):
     BendBase(right),
     messageHeader_m(" * "),
     pusher_m(right.pusher_m),
@@ -119,7 +119,7 @@ Bend::Bend(const Bend &right):
 
 }
 
-Bend::Bend(const std::string &name):
+Bend2D::Bend2D(const std::string &name):
     BendBase(name),
     messageHeader_m(" * "),
     pusher_m(),
@@ -160,7 +160,7 @@ Bend::Bend(const std::string &name):
 
 }
 
-Bend::~Bend() {
+Bend2D::~Bend2D() {
     if (entryFieldAccel_m != NULL) {
         for (unsigned int i = 0; i < 3u; ++ i) {
             gsl_spline_free(entryFieldValues_m[i]);
@@ -188,7 +188,7 @@ Bend::~Bend() {
  *  This function merely repackages the field arrays as type Vector_t and calls
  *  the equivalent method but with the Vector_t data types.
  */
-bool Bend::apply(const size_t &i,
+bool Bend2D::apply(const size_t &i,
                  const double &t,
                  Vector_t &E,
                  Vector_t &B) {
@@ -212,7 +212,7 @@ bool Bend::apply(const size_t &i,
     return false;
 }
 
-bool Bend::apply(const Vector_t &R,
+bool Bend2D::apply(const Vector_t &R,
                  const Vector_t &P,
                  const double &t,
                  Vector_t &E,
@@ -235,7 +235,7 @@ bool Bend::apply(const Vector_t &R,
 
 }
 
-bool Bend::applyToReferenceParticle(const Vector_t &R,
+bool Bend2D::applyToReferenceParticle(const Vector_t &R,
                                     const Vector_t &P,
                                     const double &t,
                                     Vector_t &E,
@@ -257,11 +257,11 @@ bool Bend::applyToReferenceParticle(const Vector_t &R,
     return false;
 }
 
-void Bend::goOnline(const double &) {
+void Bend2D::goOnline(const double &) {
     online_m = true;
 }
 
-void Bend::initialise(PartBunchBase<double, 3> *bunch,
+void Bend2D::initialise(PartBunchBase<double, 3> *bunch,
                       double &startField,
                       double &endField) {
 
@@ -305,7 +305,7 @@ void Bend::initialise(PartBunchBase<double, 3> *bunch,
     }
 }
 
-void Bend::adjustFringeFields(double ratio) {
+void Bend2D::adjustFringeFields(double ratio) {
     findChordLength(*gmsg, chordLength_m);
 
     double delta = std::abs(entranceParameter1_m - entranceParameter2_m);
@@ -323,7 +323,7 @@ void Bend::adjustFringeFields(double ratio) {
     setupFringeWidths();
 }
 
-double Bend::calculateBendAngle() {
+double Bend2D::calculateBendAngle() {
 
     const double mass = RefPartBunch_m->getM();
     const double gamma = designEnergy_m / mass + 1.0;
@@ -365,7 +365,7 @@ double Bend::calculateBendAngle() {
 
 }
 
-void Bend::calcEngeFunction(double zNormalized,
+void Bend2D::calcEngeFunction(double zNormalized,
                              const std::vector<double> &engeCoeff,
                              int polyOrder,
                              double &engeFunc,
@@ -429,7 +429,7 @@ void Bend::calcEngeFunction(double zNormalized,
     }
 }
 
-Vector_t Bend::calcCentralField(const Vector_t &R,
+Vector_t Bend2D::calcCentralField(const Vector_t &R,
                                 double deltaX) {
 
     Vector_t B(0, 0, 0);
@@ -448,7 +448,7 @@ Vector_t Bend::calcCentralField(const Vector_t &R,
     return B;
 }
 
-Vector_t Bend::calcEntranceFringeField(const Vector_t &R,
+Vector_t Bend2D::calcEntranceFringeField(const Vector_t &R,
                                        double deltaX) {
 
     const CoordinateSystemTrafo toEntranceRegion(Vector_t(0, 0, entranceParameter2_m),
@@ -485,7 +485,7 @@ Vector_t Bend::calcEntranceFringeField(const Vector_t &R,
     return toEntranceRegion.rotateFrom(B);
 }
 
-Vector_t Bend::calcExitFringeField(const Vector_t &R,
+Vector_t Bend2D::calcExitFringeField(const Vector_t &R,
                                    double deltaX) {
 
     const CoordinateSystemTrafo fromEndToExitRegion(Vector_t(0, 0, exitParameter2_m),
@@ -521,7 +521,7 @@ Vector_t Bend::calcExitFringeField(const Vector_t &R,
     return toExitRegion.rotateFrom(B);
 }
 
-bool Bend::calculateMapField(const Vector_t &R, Vector_t &B) {
+bool Bend2D::calculateMapField(const Vector_t &R, Vector_t &B) {
 
     B = Vector_t(0.0);
     bool verticallyInside = (std::abs(R(1)) < 0.5 * gap_m);
@@ -584,7 +584,7 @@ bool Bend::calculateMapField(const Vector_t &R, Vector_t &B) {
     return hitMaterial;
 }
 
-void Bend::calculateRefTrajectory(double &angleX, double &angleY) {
+void Bend2D::calculateRefTrajectory(double &angleX, double &angleY) {
 
     const double mass = RefPartBunch_m->getM();
     const double gamma = designEnergy_m / mass + 1.;
@@ -651,7 +651,7 @@ void Bend::calculateRefTrajectory(double &angleX, double &angleY) {
     angleX = -atan2(P(0), P(2)) + entranceAngle_m;
 }
 
-double Bend::estimateFieldAdjustmentStep(double actualBendAngle,
+double Bend2D::estimateFieldAdjustmentStep(double actualBendAngle,
                                           double mass,
                                           double betaGamma) {
 
@@ -669,7 +669,7 @@ double Bend::estimateFieldAdjustmentStep(double actualBendAngle,
 
 }
 
-void Bend::findBendEffectiveLength(double startField, double endField) {
+void Bend2D::findBendEffectiveLength(double startField, double endField) {
 
     /*
      * Use an iterative procedure to set the width of the
@@ -759,7 +759,7 @@ void Bend::findBendEffectiveLength(double startField, double endField) {
     }
 }
 
-void Bend::findBendStrength(double mass,
+void Bend2D::findBendStrength(double mass,
                             double gamma,
                             double betaGamma,
                             double charge) {
@@ -840,7 +840,7 @@ void Bend::findBendStrength(double mass,
     }
 }
 
-bool Bend::findIdealBendParameters(double chordLength) {
+bool Bend2D::findIdealBendParameters(double chordLength) {
 
     double refMass = RefPartBunch_m->getM();
     double refGamma = designEnergy_m / refMass + 1.0;
@@ -882,7 +882,7 @@ bool Bend::findIdealBendParameters(double chordLength) {
     return reinitialize;
 }
 
-void Bend::findReferenceExitOrigin(double &x, double &z) {
+void Bend2D::findReferenceExitOrigin(double &x, double &z) {
 
     /*
      * Find x,z coordinates of reference trajectory as it passes exit edge
@@ -903,7 +903,7 @@ void Bend::findReferenceExitOrigin(double &x, double &z) {
     }
 }
 
-bool Bend::initializeFieldMap(Inform &msg) {
+bool Bend2D::initializeFieldMap(Inform &msg) {
 
     fieldmap_m = Fieldmap::getFieldmap(fileName_m, fast_m);
 
@@ -918,7 +918,7 @@ bool Bend::initializeFieldMap(Inform &msg) {
 
 }
 
-bool Bend::inMagnetCentralRegion(const Vector_t &R) const {
+bool Bend2D::inMagnetCentralRegion(const Vector_t &R) const {
 
     Vector_t rotationCenter(-designRadius_m * cosEntranceAngle_m, R(1), designRadius_m * sinEntranceAngle_m);
     double distFromRotCenter = euclidean_norm(R - rotationCenter);
@@ -936,14 +936,14 @@ bool Bend::inMagnetCentralRegion(const Vector_t &R) const {
     return false;
 }
 
-bool Bend::inMagnetEntranceRegion(const Vector_t &R) const {
+bool Bend2D::inMagnetEntranceRegion(const Vector_t &R) const {
 
     return (R(2) >= entranceParameter1_m &&
             R(2) < 0.0 &&
             std::abs(R(0)) < aperture_m.second[0]);
 }
 
-bool Bend::inMagnetExitRegion(const Vector_t &R) const {
+bool Bend2D::inMagnetExitRegion(const Vector_t &R) const {
 
     Vector_t Rprime = getBeginToEnd_local().transformTo(R);
 
@@ -952,14 +952,14 @@ bool Bend::inMagnetExitRegion(const Vector_t &R) const {
             std::abs(Rprime(0)) < aperture_m.second[0]);
 }
 
-bool Bend::isPositionInEntranceField(const Vector_t &R) const {
+bool Bend2D::isPositionInEntranceField(const Vector_t &R) const {
 
     return (polyOrderEntry_m >= 0 &&
             R(2) >= entranceParameter1_m &&
             R(2) < entranceParameter3_m);
 }
 
-bool Bend::isPositionInExitField(const Vector_t &R) const {
+bool Bend2D::isPositionInExitField(const Vector_t &R) const {
     Vector_t Rprime = getBeginToEnd_local().transformTo(R);
 
     return (polyOrderExit_m >= 0 &&
@@ -967,7 +967,7 @@ bool Bend::isPositionInExitField(const Vector_t &R) const {
             Rprime(2) < exitParameter3_m);
 }
 
-void Bend::print(Inform &msg, double bendAngleX, double bendAngleY) {
+void Bend2D::print(Inform &msg, double bendAngleX, double bendAngleY) {
     msg << level2 << "\n"
         << "Start of field map:      "
         << startField_m
@@ -1067,7 +1067,7 @@ void Bend::print(Inform &msg, double bendAngleX, double bendAngleY) {
 }
 
 
-void Bend::readFieldMap(Inform &msg) {
+void Bend2D::readFieldMap(Inform &msg) {
     msg << level2 << getName() << " using file ";
     fieldmap_m->getInfo(&msg);
 
@@ -1147,7 +1147,7 @@ void Bend::readFieldMap(Inform &msg) {
     }
 }
 
-bool Bend::reinitialize() {
+bool Bend2D::reinitialize() {
 
     setBendStrength();
     double bendAngleX = 0.0;
@@ -1163,7 +1163,7 @@ bool Bend::reinitialize() {
     return false;
 }
 
-void Bend::setBendEffectiveLength(double startField, double endField) {
+void Bend2D::setBendEffectiveLength(double startField, double endField) {
 
     // Find initial angle.
     double actualBendAngle = calculateBendAngle();
@@ -1175,7 +1175,7 @@ void Bend::setBendEffectiveLength(double startField, double endField) {
 
 }
 
-void Bend::setBendStrength() {
+void Bend2D::setBendStrength() {
 
     // Estimate bend field magnitude.
     double mass = RefPartBunch_m->getM();
@@ -1191,7 +1191,7 @@ void Bend::setBendStrength() {
     findBendStrength(mass, gamma, betaGamma, charge);
 }
 
-void Bend::setEngeOriginDelta(double delta) {
+void Bend2D::setEngeOriginDelta(double delta) {
     /*
      * This function is used to shift the perpendicular distance of the
      * entrance and exit Enge function origins with respect to the entrance
@@ -1213,7 +1213,7 @@ void Bend::setEngeOriginDelta(double delta) {
     setupFringeWidths();
 }
 
-void Bend::setFieldCalcParam() {
+void Bend2D::setFieldCalcParam() {
 
     cosEntranceAngle_m = cos(entranceAngle_m);
     sinEntranceAngle_m = sin(entranceAngle_m);
@@ -1307,7 +1307,7 @@ void Bend::setFieldCalcParam() {
     maxAngleExitAperture -= rotationCenter;
 }
 
-void Bend::setGapFromFieldMap() {
+void Bend2D::setGapFromFieldMap() {
 
     if(gap_m <= 0.0)
         gap_m = fieldmap_m->getFieldGap();
@@ -1316,7 +1316,7 @@ void Bend::setGapFromFieldMap() {
 
 }
 
-bool Bend::setupBendGeometry(Inform &msg, double &startField, double &endField) {
+bool Bend2D::setupBendGeometry(Inform &msg, double &startField, double &endField) {
 
     chordLength_m = 0.0;
     if(!findChordLength(msg, chordLength_m))
@@ -1368,7 +1368,7 @@ bool Bend::setupBendGeometry(Inform &msg, double &startField, double &endField)
 
 }
 
-bool Bend::setupDefaultFieldMap(Inform &msg) {
+bool Bend2D::setupDefaultFieldMap(Inform &msg) {
 
     if(length_m <= 0.0) {
         ERRORMSG("If using \"1DPROFILE1-DEFAULT\" field map you must set the "
@@ -1384,7 +1384,7 @@ bool Bend::setupDefaultFieldMap(Inform &msg) {
 
 }
 
-void Bend::setFieldBoundaries(double startField, double endField) {
+void Bend2D::setFieldBoundaries(double startField, double endField) {
 
     startField_m = startField - deltaBeginEntry_m / cos(entranceAngle_m);
     endField_m = (startField + angle_m * designRadius_m +
@@ -1392,14 +1392,14 @@ void Bend::setFieldBoundaries(double startField, double endField) {
 
 }
 
-void Bend::setupPusher(PartBunchBase<double, 3> *bunch) {
+void Bend2D::setupPusher(PartBunchBase<double, 3> *bunch) {
 
     RefPartBunch_m = bunch;
     pusher_m.initialise(bunch->getReference());
 
 }
 
-bool Bend::treatAsDrift(Inform &msg, double chordLength) {
+bool Bend2D::treatAsDrift(Inform &msg, double chordLength) {
     if(designEnergy_m <= 0.0) {
         WARNMSG("Warning: bend design energy is zero. Treating as drift."
                 << endl);
@@ -1436,7 +1436,7 @@ bool Bend::treatAsDrift(Inform &msg, double chordLength) {
         return false;
 }
 
-std::vector<Vector_t> Bend::getOutline() const {
+std::vector<Vector_t> Bend2D::getOutline() const {
     std::vector<Vector_t> outline;
     Vector_t rotationCenter = Vector_t(-designRadius_m * cosEntranceAngle_m, 0.0, designRadius_m * sinEntranceAngle_m);
     unsigned int numSteps = 2;
@@ -1539,7 +1539,7 @@ std::vector<Vector_t> Bend::getOutline() const {
     return outline;
 }
 
-MeshData Bend::getSurfaceMesh() const {
+MeshData Bend2D::getSurfaceMesh() const {
     MeshData mesh;
     const Vector_t hgap(0, 0.5 * getFullGap(), 0);
     std::vector<Vector_t> outline = getOutline();
@@ -1767,7 +1767,7 @@ MeshData Bend::getSurfaceMesh() const {
     return mesh;
 }
 
-bool Bend::isInside(const Vector_t &R) const {
+bool Bend2D::isInside(const Vector_t &R) const {
     if (std::abs(R(1)) > gap_m) return false;
 
     if (inMagnetEntranceRegion(R)) {
@@ -1781,25 +1781,25 @@ bool Bend::isInside(const Vector_t &R) const {
     return (std::abs(R(1)) < 0.5 * gap_m && inMagnetCentralRegion(R));
 }
 
-void Bend::setupFringeWidths()
+void Bend2D::setupFringeWidths()
 {
     widthEntranceFringe_m = 2 * std::min(entranceParameter3_m - entranceParameter1_m, aperture_m.second[0]) + aperture_m.second[0];
     widthExitFringe_m = 2 * std::min(exitParameter3_m - exitParameter1_m, aperture_m.second[0]) + aperture_m.second[0];
 }
 
 //set the number of slices for map tracking
-void Bend::setNSlices(const std::size_t& nSlices) {
+void Bend2D::setNSlices(const std::size_t& nSlices) {
     nSlices_m = nSlices;
 }
 
 //get the number of slices for map tracking
-std::size_t Bend::getNSlices() const {
+std::size_t Bend2D::getNSlices() const {
     return nSlices_m;
 }
 
 /// Get entrance fringe field length.
 //  Used to create fringe fields in ThickTracker, (before edge[m], after edge[m])
-std::array<double,2> Bend::getEntranceFringeFieldLength() const {
+std::array<double,2> Bend2D::getEntranceFringeFieldLength() const {
     double entrParam1, entrParam2, entrParam3;
 
     fieldmap_m->get1DProfile1EntranceParam(entrParam1,
@@ -1819,7 +1819,7 @@ std::array<double,2> Bend::getEntranceFringeFieldLength() const {
 
 /// Get exit fringe field length.
 //  Used to create fringe fields in ThickTracker, (after edge[m], before edge[m])
-std::array<double,2> Bend::getExitFringeFieldLength() const {
+std::array<double,2> Bend2D::getExitFringeFieldLength() const {
     std::array<double,2> extFFL;
 
     double exitParam1, exitParam2, exitParam3;
diff --git a/src/Classic/AbsBeamline/Bend.h b/src/Classic/AbsBeamline/Bend2D.h
similarity index 88%
rename from src/Classic/AbsBeamline/Bend.h
rename to src/Classic/AbsBeamline/Bend2D.h
index 9f58ebdb322a579cf56ccb571b756719052edb85..aa87c162806ad90b7e80a1f76e5e1bb4970b78db 100644
--- a/src/Classic/AbsBeamline/Bend.h
+++ b/src/Classic/AbsBeamline/Bend2D.h
@@ -9,7 +9,7 @@
 // Copyright: see Copyright.readme
 // ------------------------------------------------------------------------
 //
-// Definitions for class: Bend
+// Definitions for class: Bend2D
 //   Defines the abstract interface for a general bend magnet.
 //
 // ------------------------------------------------------------------------
@@ -39,7 +39,7 @@ class Fieldmap;
 class MeshData;
 
 /*
- * Class Bend
+ * Class Bend2D
  *
  * Interface for general bend magnet.
  *
@@ -47,18 +47,18 @@ class MeshData;
  *
  */
 
-class Bend: public BendBase {
+class Bend2D: public BendBase {
 
 public:
 
     /// Constructor with given name.
-    explicit Bend(const std::string &name);
+    explicit Bend2D(const std::string &name);
 
-    Bend();
-    Bend(const Bend &);
-    virtual ~Bend();
+    Bend2D();
+    Bend2D(const Bend2D &);
+    virtual ~Bend2D();
 
-    /// Apply visitor to Bend.
+    /// Apply visitor to Bend2D.
     virtual void accept(BeamlineVisitor &) const = 0;
 
     /*
@@ -107,21 +107,15 @@ public:
 
     double getStartElement() const;
 
-    void setExitAngle(double exitAngle);
-
-    /*
-     * Set the name of the field map.
-     *
-     * For now this means a file that contains Enge function coefficients
-     * that describe the fringe fields at the entrance and exit.
-     */
     /// Set quadrupole field component.
     void setK1(double k1);
 
     void resetReinitializeFlag();
     void resetRecalcRefTrajFlag();
 
-    double getExitAngle() const;
+    void setExitAngle(double exitAngle);
+    virtual double getExitAngle() const;
+
     double getMapLength() const;
 
     std::vector<Vector_t> getOutline() const;
@@ -158,7 +152,7 @@ private:
     FRIEND_TEST(Maxwell, Zeros);
 #endif
     // Not implemented.
-    void operator=(const Bend &);
+    void operator=(const Bend2D &);
 
     void adjustFringeFields(double ratio);
     double calculateBendAngle();
@@ -314,120 +308,120 @@ private:
 
 
 inline
-void Bend::finalise() {
+void Bend2D::finalise() {
     online_m = false;
 }
 
 inline
-void Bend::getDimensions(double &sBegin, double &sEnd) const {
+void Bend2D::getDimensions(double &sBegin, double &sEnd) const {
     sBegin = startField_m;
     sEnd = endField_m;
 }
 
 inline
-double Bend::getBendRadius() const {
+double Bend2D::getBendRadius() const {
     return designRadius_m;
 }
 
 inline
-double Bend::getEffectiveCenter() const {
+double Bend2D::getEffectiveCenter() const {
     return elementEdge_m + designRadius_m * angle_m / 2.0;
 }
 
 inline
-double Bend::getEffectiveLength() const {
+double Bend2D::getEffectiveLength() const {
     return designRadius_m * angle_m;
 }
 
 inline
-double Bend::getStartElement() const {
+double Bend2D::getStartElement() const {
     return elementEdge_m;
 }
 
 inline
-void Bend::setK1(double k1) {
+void Bend2D::setK1(double k1) {
     if (std::abs(k1) > 0.0) {
-        throw GeneralClassicException("Bend::setK1",
+        throw GeneralClassicException("Bend2D::setK1",
                                       "Quadrupole field temporarily not supported");
     }
     fieldIndex_m = k1;
 }
 
 inline
-void Bend::resetReinitializeFlag() {
+void Bend2D::resetReinitializeFlag() {
     reinitialize_m = true;
 }
 
 inline
-void Bend::resetRecalcRefTrajFlag() {
+void Bend2D::resetRecalcRefTrajFlag() {
     recalcRefTraj_m = true;
 }
 
 inline
-void Bend::setMessageHeader(const std::string & header)
+void Bend2D::setMessageHeader(const std::string & header)
 {
     messageHeader_m = header;
 }
 
 inline
-double Bend::getStartField() const
+double Bend2D::getStartField() const
 {
     return startField_m;
 }
 
 inline
-Fieldmap* Bend::getFieldmap()
+Fieldmap* Bend2D::getFieldmap()
 {
     return fieldmap_m;
 }
 
 inline
-double Bend::getExitAngle() const
+double Bend2D::getExitAngle() const
 {
     return exitAngle_m;
 }
 
 inline
-void Bend::setExitAngle(double angle)
+void Bend2D::setExitAngle(double angle)
 {
     exitAngle_m = angle;
 }
 
 inline
-double Bend::getMapLength() const
+double Bend2D::getMapLength() const
 {
     return exitParameter2_m - entranceParameter2_m;
 }
 
 inline
-CoordinateSystemTrafo Bend::getEdgeToEnd() const
+CoordinateSystemTrafo Bend2D::getEdgeToEnd() const
 {
     return beginToEnd_m;
 }
 
 inline
-CoordinateSystemTrafo Bend::getBeginToEnd_local() const
+CoordinateSystemTrafo Bend2D::getBeginToEnd_local() const
 {
     return beginToEnd_lcs_m;
 }
 
 inline
-void Bend::setCSTrafoToEntranceRegion(const CoordinateSystemTrafo &trafo) {
+void Bend2D::setCSTrafoToEntranceRegion(const CoordinateSystemTrafo &trafo) {
     toEntranceRegion_m = trafo;
 }
 
 inline
-void Bend::setCSTrafoToExitRegion(const CoordinateSystemTrafo &trafo) {
+void Bend2D::setCSTrafoToExitRegion(const CoordinateSystemTrafo &trafo) {
     toExitRegion_m = trafo;
 }
 
 inline
-Vector_t Bend::transformToEntranceRegion(const Vector_t &R) const {
+Vector_t Bend2D::transformToEntranceRegion(const Vector_t &R) const {
     return toEntranceRegion_m.transformTo(R);
 }
 
 inline
-Vector_t Bend::transformToExitRegion(const Vector_t &R) const {
+Vector_t Bend2D::transformToExitRegion(const Vector_t &R) const {
     return toExitRegion_m.transformTo(R);
 }
 
diff --git a/src/Classic/AbsBeamline/BendBase.h b/src/Classic/AbsBeamline/BendBase.h
index c309dc86733bfd742675ca6d1c046c2957b41612..fee02884dab95887884068130b231b6d8d12b46b 100644
--- a/src/Classic/AbsBeamline/BendBase.h
+++ b/src/Classic/AbsBeamline/BendBase.h
@@ -22,7 +22,7 @@ public:
     double getBendAngle() const;
     virtual void setEntranceAngle(double entranceAngle);
     double getEntranceAngle() const;
-
+    virtual double getExitAngle() const = 0;
     void setFullGap(double);
     double getFullGap() const;
 
@@ -144,4 +144,4 @@ std::string BendBase::getFieldMapFN() const {
 }
 
 
-#endif
+#endif
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/CMakeLists.txt b/src/Classic/AbsBeamline/CMakeLists.txt
index a7cb13c1994e9413fbe3ed7f5d00bd4fcb4eb588..67ebd3db127fb8efcced622e6e1bebd824bfde3c 100644
--- a/src/Classic/AbsBeamline/CMakeLists.txt
+++ b/src/Classic/AbsBeamline/CMakeLists.txt
@@ -6,7 +6,7 @@ set (_SRCS
   AttributeSet.cpp
   BeamBeam.cpp
   BeamlineVisitor.cpp
-  Bend.cpp
+  Bend2D.cpp
   BendBase.cpp
   CCollimator.cpp
   Component.cpp
@@ -63,7 +63,7 @@ set (HDRS
     AttributeSet.h
     BeamBeam.h
     BeamlineVisitor.h
-    Bend.h
+    Bend2D.h
     BendBase.h
     CCollimator.h
     Component.h
@@ -115,4 +115,4 @@ install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline"
 
 #Turn on gcov
 #set(CMAKE_CXX_FLAGS "-g -O0 -Wall -fprofile-arcs -ftest-coverage")
-#set(CMAKE_CXX_OUTPUT_EXTENSION_REPLACE 1)
+#set(CMAKE_CXX_OUTPUT_EXTENSION_REPLACE 1)
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/RBend.cpp b/src/Classic/AbsBeamline/RBend.cpp
index c4facb937856835387d43f7ae6d79b0970bdac82..6629b5ddd97439f180fadf54b1a58569a4c2ebde 100644
--- a/src/Classic/AbsBeamline/RBend.cpp
+++ b/src/Classic/AbsBeamline/RBend.cpp
@@ -33,19 +33,19 @@ extern Inform *gmsg;
 // ------------------------------------------------------------------------
 
 RBend::RBend():
-    Bend()
+    Bend2D()
 {
     setMessageHeader("RBend ");
 }
 
 RBend::RBend(const RBend &right):
-    Bend(right)
+    Bend2D(right)
 {
     setMessageHeader("RBend ");
 }
 
 RBend::RBend(const std::string &name):
-    Bend(name)
+    Bend2D(name)
 {
     setMessageHeader("RBend ");
 }
@@ -136,12 +136,12 @@ ElementBase::ElementType RBend::getType() const {
 }
 
 void RBend::setBendAngle(double angle) {
-    Bend::setBendAngle(angle);
+    Bend2D::setBendAngle(angle);
     setExitAngle(angle - getEntranceAngle());
 }
 
 void RBend::setEntranceAngle(double entranceAngle) {
-    Bend::setEntranceAngle(entranceAngle);
+    Bend2D::setEntranceAngle(entranceAngle);
     setExitAngle(getBendAngle() - entranceAngle);
 }
 
diff --git a/src/Classic/AbsBeamline/RBend.h b/src/Classic/AbsBeamline/RBend.h
index 3965a86aac6ac132c8ab208873ce6d827cc7ad36..f35200527736b7de24d64f0becf2e59d77c475fc 100644
--- a/src/Classic/AbsBeamline/RBend.h
+++ b/src/Classic/AbsBeamline/RBend.h
@@ -21,7 +21,7 @@
 //
 // ------------------------------------------------------------------------
 
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "BeamlineGeometry/RBendGeometry.h"
 #include "Fields/BMultipoleField.h"
 
@@ -70,7 +70,7 @@
  * Here we defined the magnet as a field map.
  */
 
-class RBend: public Bend {
+class RBend: public Bend2D {
 
 public:
 
@@ -183,4 +183,4 @@ private:
                                  double &chordLength) override;
 };
 
-#endif // CLASSIC_RBend_HH
+#endif // CLASSIC_RBend_HH
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/RBend3D.h b/src/Classic/AbsBeamline/RBend3D.h
index 96689edfa3f24951537acaa76908b28318256e6b..02a69ac94686ff394b204eebe9c7a68fe4746c21 100644
--- a/src/Classic/AbsBeamline/RBend3D.h
+++ b/src/Classic/AbsBeamline/RBend3D.h
@@ -100,6 +100,7 @@ public:
 
     MeshData getSurfaceMesh() const;
 
+    virtual double getExitAngle() const;
 private:
     double trackRefParticleThrough(double dt, bool print = false);
 
@@ -120,4 +121,9 @@ private:
     void operator=(const RBend3D &);
 };
 
-#endif // CLASSIC_RBend3D_HH
+inline
+double RBend3D::getExitAngle() const {
+    return getBendAngle() - getEntranceAngle();
+}
+
+#endif // CLASSIC_RBend3D_HH
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/SBend.cpp b/src/Classic/AbsBeamline/SBend.cpp
index afb7577dec1955e43e5528ec70784cb9b9e11495..ba04adda1c7a9077d0eed84c5784516d834bf7e5 100644
--- a/src/Classic/AbsBeamline/SBend.cpp
+++ b/src/Classic/AbsBeamline/SBend.cpp
@@ -33,19 +33,19 @@ extern Inform *gmsg;
 // ------------------------------------------------------------------------
 
 SBend::SBend():
-    Bend()
+    Bend2D()
 {
     setMessageHeader("SBend ");
 }
 
 SBend::SBend(const SBend &right):
-    Bend(right)
+    Bend2D(right)
 {
     setMessageHeader("SBend ");
 }
 
 SBend::SBend(const std::string &name):
-    Bend(name)
+    Bend2D(name)
 {
     setMessageHeader("SBend ");
 }
diff --git a/src/Classic/AbsBeamline/SBend.h b/src/Classic/AbsBeamline/SBend.h
index 9f0df63156fb8dafc9aec30bb8b09788d1566e98..61e4669972465a9d24e7101acb2586a087e52d5b 100644
--- a/src/Classic/AbsBeamline/SBend.h
+++ b/src/Classic/AbsBeamline/SBend.h
@@ -21,7 +21,7 @@
 //
 // ------------------------------------------------------------------------
 
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "BeamlineGeometry/PlanarArcGeometry.h"
 #include "Fields/BMultipoleField.h"
 #include <string>
@@ -65,7 +65,7 @@
  * Here we defined the magnet as a field map.
  */
 
-class SBend: public Bend {
+class SBend: public Bend2D {
 
 public:
 
@@ -174,4 +174,4 @@ private:
 
 };
 
-#endif // CLASSIC_SBend_HH
+#endif // CLASSIC_SBend_HH
\ No newline at end of file
diff --git a/src/Classic/Solvers/CSRIGFWakeFunction.cpp b/src/Classic/Solvers/CSRIGFWakeFunction.cpp
index 474f08d84bb1924cbf2da5ca6cb72461cdd7bcbc..bb15837df9d80250354b594931ab2c345a276465 100644
--- a/src/Classic/Solvers/CSRIGFWakeFunction.cpp
+++ b/src/Classic/Solvers/CSRIGFWakeFunction.cpp
@@ -108,27 +108,18 @@ void CSRIGFWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
 }
 
 void CSRIGFWakeFunction::initialize(const ElementBase *ref) {
-    double End;
     if(ref->getType() == ElementBase::RBEND ||
        ref->getType() == ElementBase::SBEND) {
 
-        const Bend *bend = static_cast<const Bend *>(ref);
-        // const RBend *bend = dynamic_cast<const RBend *>(ref);
+        const Bend2D *bend = static_cast<const Bend2D *>(ref);
+        double End;
+
         bendRadius_m = bend->getBendRadius();
         bend->getDimensions(Begin_m, End);
         Length_m = bend->getEffectiveLength();
         FieldBegin_m = bend->getEffectiveCenter() - Length_m / 2.0;
         totalBendAngle_m = std::abs(bend->getBendAngle());
         bendName_m = bend->getName();
-
-    // } else if(dynamic_cast<const SBend *>(ref)) {
-    //     const SBend *bend = dynamic_cast<const SBend *>(ref);
-    //     bendRadius_m = bend->getBendRadius();
-    //     bend->getDimensions(Begin_m, End);
-    //     Length_m = bend->getEffectiveLength();
-    //     FieldBegin_m = bend->getEffectiveCenter() - Length_m / 2.0;
-    //     totalBendAngle_m = bend->getBendAngle();
-    //     bendName_m = bend->getName();
     }
 }
 
diff --git a/src/Classic/Solvers/CSRWakeFunction.cpp b/src/Classic/Solvers/CSRWakeFunction.cpp
index 2d792c55ddee20bb871af35361e772b6a3dfbdfb..4c6c4e06ee39582a51faa932c3a58d8548322df5 100644
--- a/src/Classic/Solvers/CSRWakeFunction.cpp
+++ b/src/Classic/Solvers/CSRWakeFunction.cpp
@@ -113,10 +113,12 @@ void CSRWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
 }
 
 void CSRWakeFunction::initialize(const ElementBase *ref) {
-    double End;
     if(ref->getType() == ElementBase::RBEND ||
        ref->getType() == ElementBase::SBEND) {
-        const Bend *bend = static_cast<const Bend *>(ref);
+
+        const Bend2D *bend = static_cast<const Bend2D *>(ref);
+        double End;
+
         bendRadius_m = bend->getBendRadius();
         bend->getDimensions(Begin_m, End);
         Length_m = bend->getEffectiveLength();
diff --git a/src/Classic/Structure/MeshGenerator.cpp b/src/Classic/Structure/MeshGenerator.cpp
index 4248da385c4b64ce4d19e30252d29e829f53b104..bf70f5eea17414b794c371886137745e7dca0ce8 100644
--- a/src/Classic/Structure/MeshGenerator.cpp
+++ b/src/Classic/Structure/MeshGenerator.cpp
@@ -1,7 +1,7 @@
 #include "Structure/MeshGenerator.h"
 #include "Physics/Physics.h"
 #include "Utilities/Util.h"
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "AbsBeamline/RBend3D.h"
 #include "AbsBeamline/Multipole.h"
 
@@ -22,7 +22,7 @@ void MeshGenerator::add(const ElementBase &element) {
     if (element.getType() == ElementBase::SBEND ||
         element.getType() == ElementBase::RBEND) {
 
-        const Bend* dipole = static_cast<const Bend*>(&element);
+        const Bend2D* dipole = static_cast<const Bend2D*>(&element);
         mesh = dipole->getSurfaceMesh();
         mesh.type_m = DIPOLE;
     } else if (element.getType() == ElementBase::RBEND3D) {
diff --git a/src/Elements/OpalBeamline.cpp b/src/Elements/OpalBeamline.cpp
index 20157a18b84dcc445edb181ec40ee12ca6acd071..216aceb9db11fd3dd6290f5bc426e362acb08991 100644
--- a/src/Elements/OpalBeamline.cpp
+++ b/src/Elements/OpalBeamline.cpp
@@ -3,7 +3,7 @@
 #include "Utilities/Options.h"
 #include "Utilities/Util.h"
 #include "AbstractObjects/OpalData.h"
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "Structure/MeshGenerator.h"
 
 #include <boost/filesystem.hpp>
@@ -400,7 +400,7 @@ void OpalBeamline::plot3DLattice() {
         std::vector<Vector_t> corners;
 
         if (element->getType() == ElementBase::RBEND || element->getType() == ElementBase::SBEND) {
-            std::vector<Vector_t> outline = static_cast<const Bend*>(element.get())->getOutline();
+            std::vector<Vector_t> outline = static_cast<const Bend2D*>(element.get())->getOutline();
 
             for (auto point: outline) {
                 corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(point)));
@@ -495,7 +495,7 @@ void OpalBeamline::save3DLattice() {
         if (element->getType() == ElementBase::SBEND ||
             element->getType() == ElementBase::RBEND) {
 
-            Bend * bendElement = static_cast<Bend*>(element.get());
+            Bend2D * bendElement = static_cast<Bend2D*>(element.get());
             std::vector<Vector_t> designPath = bendElement->getDesignPath();
             CoordinateSystemTrafo toEnd = bendElement->getBeginToEnd_local() * (*it).getCoordTransformationTo();
             Vector_t exit3D = toEnd.getOrigin();
@@ -662,7 +662,7 @@ void OpalBeamline::save3DInput() {
 
         if ((*it).getElement()->getType() == ElementBase::RBEND ||
             (*it).getElement()->getType() == ElementBase::SBEND) {
-            const Bend* dipole = static_cast<const Bend*>((*it).getElement().get());
+            const Bend2D* dipole = static_cast<const Bend2D*>((*it).getElement().get());
             double angle = dipole->getBendAngle();
             double E1 = dipole->getEntranceAngle();
             double E2 = dipole->getExitAngle();
@@ -699,7 +699,7 @@ void OpalBeamline::activateElements() {
         std::shared_ptr<Component> element = (*it).getElement();
         if (element->getType() == ElementBase::SBEND ||
             element->getType() == ElementBase::RBEND) {
-            Bend * bendElement = static_cast<Bend*>(element.get());
+            Bend2D * bendElement = static_cast<Bend2D*>(element.get());
             designEnergy = bendElement->getDesignEnergy() * 1e-6;
         }
         (*it).setOn(designEnergy);
diff --git a/src/Elements/OpalElement.cpp b/src/Elements/OpalElement.cpp
index 76b524d523e42f50fcb0a42232b8d44eb7312a4c..4a1a535c9b3150d2d94b4c239e27a2ae950f6530 100644
--- a/src/Elements/OpalElement.cpp
+++ b/src/Elements/OpalElement.cpp
@@ -21,7 +21,7 @@
 #include "Elements/OpalElement.h"
 #include "AbsBeamline/AlignWrapper.h"
 #include "AbsBeamline/ElementImage.h"
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "AbstractObjects/Attribute.h"
 #include "AbstractObjects/Expressions.h"
 #include "AbstractObjects/OpalData.h"
diff --git a/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp b/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp
index 47110b76fa0957b88ad07568c9c3e738882add16..268c9cf6d9dfaa2fbe5f2ca4f0427007e7017a82 100644
--- a/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp
+++ b/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp
@@ -1,6 +1,6 @@
 #include "gtest/gtest.h"
 #include <gtest/gtest_prod.h>
-#include "AbsBeamline/Bend.h"
+#include "AbsBeamline/Bend2D.h"
 #include "Fields/Fieldmap.h"
 #include "AbsBeamline/BeamlineVisitor.h"
 #include "BeamlineCore/RBendRep.h"
@@ -102,7 +102,7 @@ TEST(Maxwell, Zeros)
     bunch->resetM(0.938);
     bunch->setdT(1.0e-12);//time step
     double startField = 2.0, endField = 10.0 ;
-    myMagnet->Bend::initialise(bunch, startField, endField);
+    myMagnet->Bend2D::initialise(bunch, startField, endField);
     delete partData;
     delete bunch;
     Vector_t R(0. ,0., 0.);
@@ -120,9 +120,9 @@ TEST(Maxwell, Zeros)
 		  //std::cout<<"Step #"<<step<<endl;
 		  counter ++;
 		  Vector_t B(0.0);
-		  R(0) = (myMagnet->Bend::designRadius_m + x) * cos(phi);
+		  R(0) = (myMagnet->Bend2D::designRadius_m + x) * cos(phi);
 		  R(1) = z;
-		  R(2) = (myMagnet->Bend::designRadius_m + x) * sin(phi);
+		  R(2) = (myMagnet->Bend2D::designRadius_m + x) * sin(phi);
 		  double t = 0;
 		  myMagnet->apply(R, P, t , E, B);
 		  //B /= myMagnet->fieldAmplitude_m; //normalisation
@@ -147,14 +147,14 @@ TEST(Maxwell, Zeros)
 
 	        }
     //fout.close();
-    cout<<"bending radius: "<<myMagnet->Bend::designRadius_m<<endl;
+    cout<<"bending radius: "<<myMagnet->Bend2D::designRadius_m<<endl;
     cout<<"field amplitude: "<<myMagnet->fieldAmplitude_m<<endl;
     /**
     Vector_t B(0.0);
     double phi_new = Physics::pi / 20 ;
-    R[0] = (myMagnet->Bend::designRadius_m ) * cos(phi_new);
+    R[0] = (myMagnet->Bend2D::designRadius_m ) * cos(phi_new);
     R[1] = 0.005;
-    R[2] = (myMagnet->Bend::designRadius_m ) * sin(phi_new);
+    R[2] = (myMagnet->Bend2D::designRadius_m ) * sin(phi_new);
     myMagnet->apply(R, P, 0 , E, B);
     cout<<"Derivative vs stepSize:"<<endl;
     cout<<partialsDerivB(R, B,1.0, myMagnet)[1][2]<<' '<<1.0<<endl;
@@ -206,4 +206,4 @@ TEST(Quad, Quadrupole)
     gout.close();
     cout<<"length: "<<quad->getElementLength()<<endl;
 
-}
\ No newline at end of file
+}