Merge branch 'master' of gitlab.psi.ch:OPAL/src into 285-matched-distribution-not-matched

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

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

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));

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;
     }
 

src/Classic/Algorithms/PartBunch.cpp

@@ -123,8 +123,7 @@ void PartBunch::computeSelfFields(int binNumber) {
 
     if(fs_m->hasValidSolver()) {
         /// Mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         /// Scatter charge onto space charge grid.
         this->Q *= this->dt;
@@ -328,12 +327,14 @@ void PartBunch::computeSelfFields(int binNumber) {
 }
 
 void PartBunch::resizeMesh() {
+    if (fs_m->getFieldSolverType() != "SAAMG") {
+        return;
+    }
+
     double xmin = fs_m->solver_m->getXRangeMin();
     double xmax = fs_m->solver_m->getXRangeMax();
     double ymin = fs_m->solver_m->getYRangeMin();
     double ymax = fs_m->solver_m->getYRangeMax();
-    double zmin = fs_m->solver_m->getZRangeMin();
-    double zmax = fs_m->solver_m->getZRangeMax();
 
     if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
        ymin > rmin_m[1] || ymax < rmax_m[1]) {
@@ -354,14 +355,14 @@ void PartBunch::resizeMesh() {
         boundp();
         get_bounds(rmin_m, rmax_m);
     }
-    Vector_t mymin = Vector_t(xmin, ymin , zmin);
-    Vector_t mymax = Vector_t(xmax, ymax , zmax);
 
-    for (int i = 0; i < 3; i++)
-        hr_m[i]   = (mymax[i] - mymin[i])/nr_m[i];
+    Vector_t origin = Vector_t(0.0, 0.0, 0.0);
+
+    // update the mesh origin and mesh spacing hr_m
+    fs_m->solver_m->resizeMesh(origin, hr_m, rmin_m, rmax_m, dh_m);
 
     getMesh().set_meshSpacing(&(hr_m[0]));
-    getMesh().set_origin(mymin);
+    getMesh().set_origin(origin);
 
     rho_m.initialize(getMesh(),
                      getFieldLayout(),
@@ -384,8 +385,7 @@ void PartBunch::computeSelfFields() {
 
     if(fs_m->hasValidSolver()) {
         //mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         //scatter charges onto grid
         this->Q *= this->dt;
@@ -510,8 +510,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
 
     if(fs_m->hasValidSolver()) {
         /// mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         /// scatter particles charge onto grid.
         this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
@@ -639,8 +638,7 @@ void PartBunch::computeSelfFields_cycl(int bin) {
 
     if(fs_m->hasValidSolver()) {
         /// mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         /// scatter particles charge onto grid.
         this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());

src/Classic/Algorithms/PartBunchBase.h

@@ -411,6 +411,8 @@ public:
 //     virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
     virtual FieldLayout_t &getFieldLayout() = 0;
 
+    virtual void resizeMesh() { };
+
     /*
      * Wrapped member functions of IpplParticleBase
      */

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>);

src/Distribution/Distribution.cpp

@@ -17,7 +17,6 @@
 #include <string>
 #include <vector>
 #include <numeric>
-#include <limits>
 
 // IPPL
 #include "DataSource/DataConnect.h"
@@ -875,14 +874,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);
@@ -1069,9 +1060,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;
@@ -1422,7 +1413,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;
@@ -3650,9 +3641,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);
     }
 }
 
@@ -3985,14 +3976,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());
@@ -4146,9 +4137,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];
@@ -4429,8 +4420,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];

src/Distribution/Distribution.h

@@ -292,7 +292,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 {

src/Elements/OpalElement.cpp

@@ -568,7 +568,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");
             }
         }
@@ -577,7 +578,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");
             }
         }

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;
 }
 
 
@@ -131,7 +128,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]));
@@ -177,8 +174,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]));

src/Elements/OpalRBend3D.cpp

@@ -72,10 +72,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) {
@@ -130,8 +128,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]));

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;
 }
 
 
@@ -131,7 +128,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]));
@@ -163,11 +160,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]));
@@ -183,15 +180,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::fillRegisteredAttributes",
-                            "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]);

src/Solvers/EllipticDomain.h

 //
 // Class EllipticDomain
-//   :FIXME: add brief description
+//   This class provides an elliptic beam pipe. The mesh adapts to the bunch size
+//   in the longitudinal direction. At the intersection of the mesh with the
+//   beam pipe, three stencil interpolation methods are available.
 //
 // Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
 //               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
@@ -32,100 +34,127 @@
 #include <cmath>
 #include "IrregularDomain.h"
 #include "Structure/BoundaryGeometry.h"
+#include "Utilities/OpalException.h"
 
 class EllipticDomain : public IrregularDomain {
 
 public:
 
     EllipticDomain(Vector_t nr, Vector_t hr);
-    EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl);
-    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
+
+    EllipticDomain(double semimajor, double semiminor, Vector_t nr,
+                   Vector_t hr, std::string interpl);
+
+    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr,
+                   Vector_t hr, std::string interpl);
+
     EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
 
     ~EllipticDomain();
 
     /// returns discretization at (x,y,z)
-    void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
+    void getBoundaryStencil(int x, int y, int z,
+                            double &W, double &E, double &S,
+                            double &N, double &F, double &B,
+                            double &C, double &scaleFactor);
+
     /// returns discretization at 3D index
-    void getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
+    void getBoundaryStencil(int idx, double &W, double &E,
+                            double &S, double &N, double &F,
+                            double &B, double &C, double &scaleFactor);
+
     /// returns index of neighbours at (x,y,z)
-    void getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B);
+    void getNeighbours(int x, int y, int z, int &W, int &E,
+                       int &S, int &N, int &F, int &B);
+
     /// returns index of neighbours at 3D index
     void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
+
     /// returns type of boundary condition
     std::string getType() {return "Elliptic";}
+
     /// queries if a given (x,y,z) coordinate lies inside the domain
     inline bool isInside(int x, int y, int z) {
-        double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
-        double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
-        return ((xx * xx / (SemiMajor * SemiMajor) + yy * yy / (SemiMinor * SemiMinor) < 1) && z != 0 && z != nr[2] - 1);
+        double xx = - semiMajor_m + hr[0] * (x + 0.5);
+        double yy = - semiMinor_m + hr[1] * (y + 0.5);
+
+        bool isInsideEllipse = (xx * xx / (semiMajor_m * semiMajor_m) +
+                                yy * yy / (semiMinor_m * semiMinor_m) < 1);
+
+        return (isInsideEllipse && z > 0 && z < nr[2] - 1);
     }
 
     int getNumXY(int /*z*/) { return nxy_m; }
-    /// set semi-minor
-    void setSemiMinor(double sm) {SemiMinor = sm;}
-    /// set semi-major
-    void setSemiMajor(double sm) {SemiMajor = sm;}
 
-    /// calculates intersection 
-    void compute(Vector_t hr);
+
+    /// calculates intersection
+    void compute(Vector_t /*hr*/) {
+        throw OpalException("EllipticDomain::compute()", "This function is not available.");
+    }