diff --git a/gen_OPALrevision b/gen_OPALrevision
index 78bd94b8f711d417c525d7fac0fd3e5e43576da7..9f298fa7937a190611be31e9944591faf263e8e7 100755
--- a/gen_OPALrevision
+++ b/gen_OPALrevision
@@ -4,7 +4,12 @@
 #
 
 print () {
-    echo '#define GIT_VERSION '\"$1\" > src/OPALrevision.h
+    echo '#define GIT_VERSION '\"$1\" > src/OPALrevision-new.h
+    if [[ -e  src/OPALrevision.h ]] && cmp -s src/OPALrevision.h src/OPALrevision-new.h; then
+        rm src/OPALrevision-new.h
+    else
+        mv src/OPALrevision-new.h src/OPALrevision.h
+    fi
 }
 
 # if git(1) is not in $PATH
diff --git a/src/Algorithms/CavityAutophaser.cpp b/src/Algorithms/CavityAutophaser.cpp
index 9fa6b29ecacc0e848295ab56a1fcd8a27e7dc419..4c791f71400bc77f3267baced2efad76627cac83 100644
--- a/src/Algorithms/CavityAutophaser.cpp
+++ b/src/Algorithms/CavityAutophaser.cpp
@@ -72,7 +72,7 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
                    << std::left << std::setw(68) << std::setfill('*') << ss.str()
                    << std::setfill(' ') << endl);
     if (!apVeto) {
-        double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
+        double initialEnergy = Util::getKineticEnergy(P, itsReference_m.getM()) * 1e-6;
         double AstraPhase    = 0.0;
         double designEnergy  = element->getDesignEnergy();
 
@@ -191,7 +191,7 @@ double CavityAutophaser::guessCavityPhase(double t) {
         return orig_phi;
     }
 
-    Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
+    Phimax = element->getAutoPhaseEstimate(Util::getKineticEnergy(refP, itsReference_m.getM()) * 1e-6,
                                            t,
                                            itsReference_m.getQ(),
                                            itsReference_m.getM() * 1e-6);
@@ -288,7 +288,7 @@ double CavityAutophaser::track(Vector_t /*R*/,
                                                             out);
     rfc->setPhasem(initialPhase);
 
-    double finalKineticEnergy = Util::getEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
+    double finalKineticEnergy = Util::getKineticEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
 
     return finalKineticEnergy;
 }
\ No newline at end of file
diff --git a/src/Algorithms/CavityAutophaser.h b/src/Algorithms/CavityAutophaser.h
index 29db048b24349a39a7e9b93b5518d18aafc5571a..a4e49c8da834b0c8cd8d8010aea35439c4315403 100644
--- a/src/Algorithms/CavityAutophaser.h
+++ b/src/Algorithms/CavityAutophaser.h
@@ -28,8 +28,6 @@ private:
                  const double dt,
                  const double phase,
                  std::ofstream *out = NULL) const;
-    double getEnergyMeV(const Vector_t &P);
-    double getMomentum(double kineticEnergyMeV);
 
     const PartData &itsReference_m;
     std::shared_ptr<Component> itsCavity_m;
@@ -39,14 +37,4 @@ private:
 
 };
 
-inline
-double CavityAutophaser::getEnergyMeV(const Vector_t &P) {
-    return itsReference_m.getM() * 1e-6 * (sqrt(dot(P,P) + 1) - 1);
-}
-
-inline
-double CavityAutophaser::getMomentum(double kineticEnergyMeV) {
-    return sqrt(std::pow(kineticEnergyMeV / (itsReference_m.getM() * 1e-6) + 1, 2) - 1);
-}
-
 #endif
\ No newline at end of file
diff --git a/src/Algorithms/ParallelTTracker.cpp b/src/Algorithms/ParallelTTracker.cpp
index f38da15e80bd5e36a201c9a9637a7aeaee16c12c..23e2a6c2fdfb97e3d3d8cec97a629ab3756b0452 100644
--- a/src/Algorithms/ParallelTTracker.cpp
+++ b/src/Algorithms/ParallelTTracker.cpp
@@ -1162,7 +1162,7 @@ void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
     selectDT();
 
     if ((back_track && itsOpalBeamline_m.containsSource()) ||
-        Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
+        Util::getKineticEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
         double gamma = 0.1 / itsBunch_m->getM() + 1.0;
         double beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
         itsBunch_m->RefPartP_m = itsBunch_m->toLabTrafo_m.rotateTo(beta * gamma * Vector_t(0, 0, 1));
diff --git a/src/Classic/AbsBeamline/RFCavity.cpp b/src/Classic/AbsBeamline/RFCavity.cpp
index b8ebb79459043ae993e9c3b3b05712ad5183e541..33cecafa90dc8f8fee454843f7be5d0b6e1894e3 100644
--- a/src/Classic/AbsBeamline/RFCavity.cpp
+++ b/src/Classic/AbsBeamline/RFCavity.cpp
@@ -494,22 +494,22 @@ ElementBase::ElementType RFCavity::getType() const {
 double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
 
     const double dt = 1e-13;
-    const double p0 = Util::getP(E0, mass);
+    const double p0 = Util::getBetaGamma(E0, mass);
     const double origPhase =getPhasem();
     double dphi = Physics::pi / 18;
 
     double phi = 0.0;
     setPhasem(phi);
-    std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
+    std::pair<double, double> ret = trackOnAxisParticle(p0, t0, dt, q, mass);
     double phimax = 0.0;
-    double Emax = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
+    double Emax = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
     phi += dphi;
 
     for (unsigned int j = 0; j < 2; ++ j) {
         for (unsigned int i = 0; i < 36; ++ i, phi += dphi) {
             setPhasem(phi);
             ret = trackOnAxisParticle(p0, t0, dt, q, mass);
-            double Ekin = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
+            double Ekin = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
             if (Ekin > Emax) {
                 Emax = Ekin;
                 phimax = phi;
@@ -635,7 +635,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
         }
 
         double cosine_part = 0.0, sine_part = 0.0;
-        double p0 = std::sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
+        double p0 = Util::getBetaGamma(E0, mass);
         cosine_part += scale_m * std::cos(frequency_m * t0) * F[0];
         sine_part += scale_m * std::sin(frequency_m * t0) * F[0];
 
@@ -675,11 +675,11 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
     const double zend = length_m + startField_m;
 
     Vector_t z(0.0, 0.0, zbegin);
-    double dz = 0.5 * p(2) / std::sqrt(1.0 + dot(p, p)) * cdt;
+    double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;
     Vector_t Ef(0.0), Bf(0.0);
 
     if (out) *out << std::setw(18) << z[2]
-                  << std::setw(18) << Util::getEnergy(p, mass)
+                  << std::setw(18) << Util::getKineticEnergy(p, mass)
                   << std::endl;
     while (z(2) + dz < zend && z(2) + dz > zbegin) {
         z /= cdt;
@@ -700,7 +700,7 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
         t += dt;
 
         if (out) *out << std::setw(18) << z[2]
-                      << std::setw(18) << Util::getEnergy(p, mass)
+                      << std::setw(18) << Util::getKineticEnergy(p, mass)
                       << std::endl;
     }
 
diff --git a/src/Classic/Utilities/Util.h b/src/Classic/Utilities/Util.h
index edf7178517c527c010c2bfcfd183278465f0fd1a..42479de4b024055ced8e0e5cc28a80a702a933e2 100644
--- a/src/Classic/Utilities/Util.h
+++ b/src/Classic/Utilities/Util.h
@@ -6,6 +6,7 @@
 
 #include <string>
 #include <cstring>
+#include <limits>
 #include <sstream>
 #include <type_traits>
 #include <functional>
@@ -28,16 +29,23 @@ namespace Util {
     }
 
     inline
-    double getEnergy(Vector_t p, double mass) {
+    double getKineticEnergy(Vector_t p, double mass) {
         return (getGamma(p) - 1.0) * mass;
     }
 
     inline
-    double getP(double E, double mass) {
-        double gamma = E / mass + 1;
-        return std::sqrt(std::pow(gamma, 2.0) - 1.0);
+    double getBetaGamma(double Ekin, double mass) {
+        double value = std::sqrt(std::pow(Ekin / mass + 1.0, 2) - 1.0);
+        if (value < std::numeric_limits<double>::epsilon())
+            value = std::sqrt(2 * Ekin / mass);
+        return value;
     }
 
+    inline
+    double convertMomentumeVToBetaGamma(double p, double mass) {
+        return p / mass;
+    }
+  
     inline
     std::string getTimeString(double time, unsigned int precision = 3) {
         std::string timeUnit(" [ps]");
@@ -162,7 +170,7 @@ namespace Util {
     }
 
     Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName = "");
-
+  
     std::string toUpper(const std::string &str);
 
     std::string combineFilePath(std::initializer_list<std::string>);
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 6e39d527c85de989fe58527fb1a4fa7c00f4fec1..89976348f672e09d09d2c3a1450bbf6fe475da37 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -17,7 +17,6 @@
 #include <string>
 #include <vector>
 #include <numeric>
-#include <limits>
 
 // IPPL
 #include "DataSource/DataConnect.h"
@@ -856,21 +855,18 @@ void Distribution::checkParticleNumber(size_t &numberOfParticles) {
     }
 
     if (numberOfDistParticles != numberOfParticles) {
-        *gmsg << "\n--------------------------------------------------" << endl
-              << "Warning!! The number of particles in the initial" << endl
-              << "distribution is " << numberOfDistParticles << "." << endl << endl
-              << "This is different from the number of particles" << endl
-              << "defined by the BEAM command: " << numberOfParticles << endl << endl
-              << "This is often happens when using a FROMFILE type" << endl
-              << "distribution and not matching the number of" << endl
-              << "particles in the particles file(s) with the number" << endl
-              << "given in the BEAM command." << endl << endl
-              << "The number of particles in the initial distribution" << endl
-              << "(" << numberOfDistParticles << ") "
-              << "will take precedence." << endl
-              << "---------------------------------------------------\n" << endl;
+        throw OpalException("Distribution::checkParticleNumber",
+                            "The number of particles in the initial\n"
+                            "distribution " +
+                            std::to_string(numberOfDistParticles) + "\n"
+                            "is different from the number of particles\n"
+                            "defined by the BEAM command " +
+                            std::to_string(numberOfParticles) + ".\n"
+                            "This often happens when using a FROMFILE type\n"
+                            "distribution and not matching the number of\n"
+                            "particles in the input distribution file(s) with\n"
+                            "the number given in the BEAM command.");
     }
-    numberOfParticles = numberOfDistParticles;
 }
 
 void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits) {
@@ -888,14 +884,6 @@ void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUn
 
 }
 
-double Distribution::converteVToBetaGamma(double valueIneV, double massIneV) {
-    double value = std::copysign(std::sqrt(std::pow(std::abs(valueIneV) / massIneV + 1.0, 2) - 1.0), valueIneV);
-    if (std::abs(value) < std::numeric_limits<double>::epsilon())
-        value = std::copysign(std::sqrt(2 * std::abs(valueIneV) / massIneV), valueIneV);
-
-    return value;
-}
-
 void Distribution::createDistributionBinomial(size_t numberOfParticles, double massIneV) {
 
     setDistParametersBinomial(massIneV);
@@ -1082,9 +1070,9 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
                 saveProcessor = 0;
 
             if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-                P(0) = converteVToBetaGamma(P(0), massIneV);
-                P(1) = converteVToBetaGamma(P(1), massIneV);
-                P(2) = converteVToBetaGamma(P(2), massIneV);
+                P(0) = Util::convertMomentumeVToBetaGamma(P(0), massIneV);
+                P(1) = Util::convertMomentumeVToBetaGamma(P(1), massIneV);
+                P(2) = Util::convertMomentumeVToBetaGamma(P(2), massIneV);
             }
 
             pmean_m += P;
@@ -1483,7 +1471,7 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
     // This is PC from BEAM
     double deltaP = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETP]);
     if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-        deltaP = converteVToBetaGamma(deltaP, beam->getM());
+        deltaP = Util::convertMomentumeVToBetaGamma(deltaP, beam->getM());
     }
 
     avrgpz_m = beam->getP()/beam->getM() + deltaP;
@@ -3580,9 +3568,9 @@ void Distribution::setSigmaP_m(double massIneV) {
 
     // Check what input units we are using for momentum.
     if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-        sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
-        sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
-        sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
+        sigmaP_m[0] = Util::convertMomentumeVToBetaGamma(sigmaP_m[0], massIneV);
+        sigmaP_m[1] = Util::convertMomentumeVToBetaGamma(sigmaP_m[1], massIneV);
+        sigmaP_m[2] = Util::convertMomentumeVToBetaGamma(sigmaP_m[2], massIneV);
     }
 }
 
@@ -3915,14 +3903,14 @@ void Distribution::setupEmissionModel(PartBunchBase<double, 3> *beam) {
 void Distribution::setupEmissionModelAstra(PartBunchBase<double, 3> *beam) {
 
     double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
-    pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
+    pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
     pmean_m = Vector_t(0.0, 0.0, 0.5 * pTotThermal_m);
 }
 
 void Distribution::setupEmissionModelNone(PartBunchBase<double, 3> *beam) {
 
     double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
-    pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
+    pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
     double avgPz = std::accumulate(pzDist_m.begin(), pzDist_m.end(), 0.0);
     size_t numParticles = pzDist_m.size();
     reduce(avgPz, avgPz, OpAddAssign());
@@ -4076,9 +4064,9 @@ void Distribution::shiftDistCoordinates(double massIneV) {
 
         // Check input momentum units.
         if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-            deltaPx = converteVToBetaGamma(deltaPx, massIneV);
-            deltaPy = converteVToBetaGamma(deltaPy, massIneV);
-            deltaPz = converteVToBetaGamma(deltaPz, massIneV);
+            deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
+            deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
+            deltaPz = Util::convertMomentumeVToBetaGamma(deltaPz, massIneV);
         }
 
         size_t endIdx = startIdx + particlesPerDist_m[i];
@@ -4359,8 +4347,8 @@ void Distribution::adjustPhaseSpace(double massIneV) {
     double deltaPy = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETPY]);
     // Check input momentum units.
     if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-        deltaPx = converteVToBetaGamma(deltaPx, massIneV);
-        deltaPy = converteVToBetaGamma(deltaPy, massIneV);
+        deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
+        deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
     }
 
     double avrg[6];
diff --git a/src/Distribution/Distribution.h b/src/Distribution/Distribution.h
index 159669467fcb19f09bfd6ec40acee7ad0a76cef0..7eda65809a5b2f0e968abf155fb6b4a5c64833bb 100644
--- a/src/Distribution/Distribution.h
+++ b/src/Distribution/Distribution.h
@@ -291,7 +291,6 @@ private:
     void checkIfEmitted();
     void checkParticleNumber(size_t &numberOfParticles);
     void chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits);
-    double converteVToBetaGamma(double valueIneV, double massIneV);
     size_t getNumberOfParticlesInFile(std::ifstream &inputFile);
 
     class BinomialBehaviorSplitter {
diff --git a/src/Elements/OpalElement.cpp b/src/Elements/OpalElement.cpp
index 9d31e27baaa5d34f63e89113bad5a8f7b5937c9a..1188416ba8a98e4cadd0f57adf39d467676f24c7 100644
--- a/src/Elements/OpalElement.cpp
+++ b/src/Elements/OpalElement.cpp
@@ -454,7 +454,8 @@ void OpalElement::update() {
             rotation = rotTheta * (rotPhi * rotPsi);
         } else {
             if (itsAttr[ORIENTATION]) {
-                throw OpalException("Line::parse","Parameter orientation is array of 3 values (theta, phi, psi);\n" +
+                throw OpalException("OpalElement::update",
+                                    "Parameter orientation is array of 3 values (theta, phi, psi);\n" +
                                     std::to_string(dir.size()) + " values provided");
             }
         }
@@ -463,7 +464,8 @@ void OpalElement::update() {
             origin = Vector_t(ori[0], ori[1], ori[2]);
         } else {
             if (itsAttr[ORIGIN]) {
-                throw OpalException("Line::parse","Parameter origin is array of 3 values (x, y, z);\n" +
+                throw OpalException("OpalElement::update",
+                                    "Parameter origin is array of 3 values (x, y, z);\n" +
                                     std::to_string(ori.size()) + " values provided");
             }
         }
diff --git a/src/Elements/OpalRBend.cpp b/src/Elements/OpalRBend.cpp
index 0b453b0d22086aec5b51b533e94875445e4db90d..27437e7b5f89ff9534f099f8b90b93a3cc5ae40b 100644
--- a/src/Elements/OpalRBend.cpp
+++ b/src/Elements/OpalRBend.cpp
@@ -24,7 +24,6 @@
 #include "Structure/OpalWake.h"
 #include "Structure/ParticleMatterInteraction.h"
 #include "Utilities/OpalException.h"
-#include <cmath>
 
 OpalRBend::OpalRBend():
     OpalBend("RBEND",
@@ -47,10 +46,8 @@ OpalRBend::OpalRBend(const std::string &name, OpalRBend *parent):
 
 
 OpalRBend::~OpalRBend() {
-    if(owk_m)
-        delete owk_m;
-    if(parmatint_m)
-        delete parmatint_m;
+    delete owk_m;
+    delete parmatint_m;
 }
 
 
@@ -97,7 +94,7 @@ void OpalRBend::update() {
     double k0s = itsAttr[K0S] ? Attributes::getReal(itsAttr[K0S]) : 0.0;
     //JMJ 4/10/2000: above line replaced
     //    length ? angle / length : angle;
-    // to avoid closed orbit created by RBEND with defalt K0.
+    // to avoid closed orbit created by RBEND with default K0.
     field.setNormalComponent(1, factor * k0);
     field.setSkewComponent(1, factor * Attributes::getReal(itsAttr[K0S]));
     field.setNormalComponent(2, factor * Attributes::getReal(itsAttr[K1]));
@@ -143,8 +140,11 @@ void OpalRBend::update() {
     }
 
     // Energy in eV.
-    if(itsAttr[DESIGNENERGY]) {
+    if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
         bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
+    } else if (bend->getName() != "RBEND") {
+        throw OpalException("OpalRBend::update",
+                            "RBend requires non-zero DESIGNENERGY");
     }
 
     bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
diff --git a/src/Elements/OpalRBend3D.cpp b/src/Elements/OpalRBend3D.cpp
index 59aaf0162b0881c550babd959664b3eff2d6d0aa..8c6773e62346acf38fae67ba44d4048ae577cf21 100644
--- a/src/Elements/OpalRBend3D.cpp
+++ b/src/Elements/OpalRBend3D.cpp
@@ -62,10 +62,8 @@ OpalRBend3D::OpalRBend3D(const std::string &name, OpalRBend3D *parent):
 }
 
 OpalRBend3D::~OpalRBend3D() {
-    if(owk_m)
-        delete owk_m;
-    if(parmatint_m)
-        delete parmatint_m;
+    delete owk_m;
+    delete parmatint_m;
 }
 
 OpalRBend3D *OpalRBend3D::clone(const std::string &name) {
@@ -115,8 +113,11 @@ void OpalRBend3D::update() {
     bend->setEntranceAngle(e1);
 
     // Energy in eV.
-    if(itsAttr[DESIGNENERGY]) {
+    if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
         bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
+    } else if (bend->getName() != "RBEND3D") {
+        throw OpalException("OpalRBend3D::update",
+                            "RBend3D requires non-zero DESIGNENERGY");
     }
 
     bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
diff --git a/src/Elements/OpalSBend.cpp b/src/Elements/OpalSBend.cpp
index 02c4275ad4f856d4fc736058cc0eaa0466e2602c..ccc3c2da8f03e81b3659e54e58e9e0e7a438a3ed 100644
--- a/src/Elements/OpalSBend.cpp
+++ b/src/Elements/OpalSBend.cpp
@@ -24,7 +24,6 @@
 #include "Structure/OpalWake.h"
 #include "Structure/ParticleMatterInteraction.h"
 #include "Utilities/OpalException.h"
-#include <cmath>
 
 OpalSBend::OpalSBend():
     OpalBend("SBEND",
@@ -47,10 +46,8 @@ OpalSBend::OpalSBend(const std::string &name, OpalSBend *parent):
 
 
 OpalSBend::~OpalSBend() {
-    if(owk_m)
-        delete owk_m;
-    if(parmatint_m)
-        delete parmatint_m;
+    delete owk_m;
+    delete parmatint_m;
 }
 
 
@@ -96,7 +93,7 @@ void OpalSBend::update() {
     double k0s = itsAttr[K0S] ? Attributes::getReal(itsAttr[K0S]) : 0.0;
     //JMJ 4/10/2000: above line replaced
     //    length ? angle / length : angle;
-    // to avoid closed orbit created by SBEND with defalt K0.
+    // to avoid closed orbit created by SBEND with default K0.
     field.setNormalComponent(1, factor * k0);
     field.setSkewComponent(1, factor * Attributes::getReal(itsAttr[K0S]));
     field.setNormalComponent(2, factor * Attributes::getReal(itsAttr[K1]));
@@ -128,11 +125,11 @@ void OpalSBend::update() {
 
     if(itsAttr[GREATERTHANPI])
         throw OpalException("OpalSBend::update",
-                            "GREATERTHANPI not supportet any more");
+                            "GREATERTHANPI not supported anymore");
 
     if(itsAttr[ROTATION])
         throw OpalException("OpalSBend::update",
-                            "ROTATION not supportet any more; use PSI instead");
+                            "ROTATION not supported anymore; use PSI instead");
 
     if(itsAttr[FMAPFN])
         bend->setFieldMapFN(Attributes::getString(itsAttr[FMAPFN]));
@@ -148,15 +145,18 @@ void OpalSBend::update() {
     bend->setExitAngle(e2);
 
     // Units are eV.
-    if(itsAttr[DESIGNENERGY]) {
+    if(itsAttr[DESIGNENERGY] && Attributes::getReal(itsAttr[DESIGNENERGY]) != 0.0) {
         bend->setDesignEnergy(Attributes::getReal(itsAttr[DESIGNENERGY]), false);
+    } else if (bend->getName() != "SBEND") {
+        throw OpalException("OpalSBend::update",
+                            "SBend requires non-zero DESIGNENERGY");
     }
 
     bend->setFullGap(Attributes::getReal(itsAttr[GAP]));
 
     if(itsAttr[APERT])
-        throw OpalException("OpalRBend::update",
-                            "APERTURE in RBEND not supported; use GAP and HAPERT instead");
+        throw OpalException("OpalSBend::update",
+                            "APERTURE in SBEND not supported; use GAP and HAPERT instead");
 
     if(itsAttr[HAPERT]) {
         double hapert = Attributes::getReal(itsAttr[HAPERT]);
diff --git a/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp b/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp
index e6b4a98ddd711a46cac0ae49fd5612973f4bba94..b839e4d8009c6adc1e2a0670f7400b738bfe66c8 100644
--- a/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp
+++ b/tests/classic_src/Solvers/GreenWakeFunctionTest.cpp
@@ -25,10 +25,10 @@ TEST(GreenWakeFunctionTest, TestApply)
     bool const_length = true;
     std::string fname = "";
 
-    std::vector<double> finalWakeValues   = {1.61757e+10,  2.78525e+19};
-    std::vector<double> finalEnergyValues = {-0.00125303, -2.15739e+06};
+    std::vector<double> finalWakeValues   = { 1.61757e+10,  2.78525e+19};
+    std::vector<double> finalEnergyValues = {-2.446072e-5, -42118.09};
     std::vector<double> relativeErrorWake   = {1e+6, 1e+15};
-    std::vector<double> relativeErrorEnergy = {1e-7, 1e+2};
+    std::vector<double> relativeErrorEnergy = {1e-9, 1};
 
     for (int acmode : acmodes) {
         GreenWakeFunction gwf("opal", nullptr, filters, nbin, Z0, radius, sigma, acmode, tau, 0, const_length, fname);
@@ -37,7 +37,7 @@ TEST(GreenWakeFunctionTest, TestApply)
         // determine K and charge
         double charge = 0.8e-9; // nC
         double K = 0.20536314319923724e-9; //K normalizes nC data in lambda.h?
-        gwf.NBin_m = 294;
+        gwf.NBin_m = 10;
 
         std::cout << "# Z0 = "        << gwf.Z0_m        << std::endl
                   << "# radius = "    << gwf.radius_m    << std::endl