diff --git a/ippl/src/AppTypes/Vektor.h b/ippl/src/AppTypes/Vektor.h index 92f1fcd9a6f314db9f85d655d93c87c4f337fd6a..7f92412ccaaa370542d42311b4593044958bb655 100644 --- a/ippl/src/AppTypes/Vektor.h +++ b/ippl/src/AppTypes/Vektor.h @@ -40,10 +40,6 @@ public: TSV_MetaAssignScalar<Vektor<T,D>,T,OpAssign>::apply(*this,T(0)); } - // A noninitializing ctor. - class DontInitialize {}; - Vektor(DontInitialize) {} - // Copy Constructor Vektor(const Vektor<T,D> &rhs) { TSV_MetaAssign< Vektor<T,D> , Vektor<T,D> ,OpAssign >::apply(*this,rhs); diff --git a/ippl/src/Message/Message.h b/ippl/src/Message/Message.h index 392540111fe8b07fde4c4554ae7bdcdcf543e41f..4505bf58a9bd479e5125392d3d85b03473fac835 100644 --- a/ippl/src/Message/Message.h +++ b/ippl/src/Message/Message.h @@ -211,7 +211,7 @@ public: } // copy constructor - MsgItem(const MsgItem &m) : item(&defbuf), elements(m.elements), bytesize(m.bytesize), + MsgItem(const MsgItem &m) : item(&defbuf), elements(m.elements), bytesize(m.bytesize), defbuf(m.defbuf), needDelete(m.needDelete) { // either we just copy the 'item' pointer, or we copy the default buf diff --git a/optimizer/Optimizer/EA/Variator.h b/optimizer/Optimizer/EA/Variator.h index 8cb72b35bbba93566c1724ce03d30d39a6bf287c..77183e0cb370c19365fc48594c6a1dd2a00e6640 100644 --- a/optimizer/Optimizer/EA/Variator.h +++ b/optimizer/Optimizer/EA/Variator.h @@ -213,7 +213,7 @@ public: if (viableA == true && viableB == true) { break; } - std::cout << "Individual not viable, I try again: iter= " << iter << std::endl; + std::cout << "Individuals not viable, I try again: iter= " << iter << std::endl; iter++; // if maximum number of tries then create new individual(s) if (iter > 100) { diff --git a/src/BasicActions/DumpFields.cpp b/src/BasicActions/DumpFields.cpp index 2bbcf654de568897cdb622dfc477b0712874c50b..b9c029e1f84d22442ef1f2a2b769164d1bc3a453 100644 --- a/src/BasicActions/DumpFields.cpp +++ b/src/BasicActions/DumpFields.cpp @@ -34,8 +34,6 @@ #include "Attributes/Attributes.h" #include "BasicActions/DumpFields.h" -// extern Inform *gmsg; - std::unordered_set<DumpFields*> DumpFields::dumpsSet_m; std::string DumpFields::dumpfields_docstring = @@ -116,9 +114,7 @@ void DumpFields::buildGrid() { checkInt(nx, "X_STEPS"); checkInt(ny, "Y_STEPS"); checkInt(nz, "Z_STEPS"); - if (grid_m != NULL) { - delete grid_m; - } + delete grid_m; grid_m = new interpolation::ThreeDGrid(dx, dy, dz, x0, y0, z0, @@ -148,11 +144,11 @@ void DumpFields::checkInt(double real, std::string name, double tolerance) { void DumpFields::writeFieldThis(Component* field) { if (grid_m == NULL) { - throw OpalException("DumpFields::writeField", + throw OpalException("DumpFields::writeFieldThis", "The grid was NULL; there was a problem with the DumpFields initialisation."); } if (field == NULL) { - throw OpalException("DumpFields::writeField", + throw OpalException("DumpFields::writeFieldThis", "The field to be written was NULL."); } double time = 0.; @@ -160,7 +156,7 @@ void DumpFields::writeFieldThis(Component* field) { Vector_t centroid(0., 0., 0.); std::ofstream fout(filename_m.c_str(), std::ofstream::out); if (!fout.good()) { - throw OpalException("DumpFields::writeField", + throw OpalException("DumpFields::writeFieldThis", "Failed to open DumpFields file "+filename_m); } // set precision @@ -183,7 +179,7 @@ void DumpFields::writeFieldThis(Component* field) { fout << B[0] << " " << B[1] << " " << B[2] << "\n"; } if (!fout.good()) { - throw OpalException("DumpFields::writeField", + throw OpalException("DumpFields::writeFieldThis", "Something went wrong during writing "+filename_m); } fout.close(); diff --git a/src/Classic/AbsBeamline/Bend2D.cpp b/src/Classic/AbsBeamline/Bend2D.cpp index 14554e0712988f510edff40cace87220ed66eb0c..53b138bdc9a6b9b5352fe3d03347821295b872c3 100644 --- a/src/Classic/AbsBeamline/Bend2D.cpp +++ b/src/Classic/AbsBeamline/Bend2D.cpp @@ -46,20 +46,22 @@ Bend2D::Bend2D(const Bend2D &right): messageHeader_m(" * "), pusher_m(right.pusher_m), fieldmap_m(right.fieldmap_m), - fast_m(right.fast_m), designRadius_m(right.designRadius_m), exitAngle_m(right.exitAngle_m), fieldIndex_m(right.fieldIndex_m), startField_m(right.startField_m), endField_m(right.endField_m), + widthEntranceFringe_m(right.widthEntranceFringe_m), + widthExitFringe_m(right.widthExitFringe_m), reinitialize_m(right.reinitialize_m), - recalcRefTraj_m(right.recalcRefTraj_m), entranceParameter1_m(right.entranceParameter1_m), entranceParameter2_m(right.entranceParameter2_m), entranceParameter3_m(right.entranceParameter3_m), exitParameter1_m(right.exitParameter1_m), exitParameter2_m(right.exitParameter2_m), exitParameter3_m(right.exitParameter3_m), + engeCoeffsEntry_m(right.engeCoeffsEntry_m), + engeCoeffsExit_m(right.engeCoeffsExit_m), entryFieldValues_m(NULL), exitFieldValues_m(NULL), entryFieldAccel_m(NULL), @@ -67,8 +69,6 @@ Bend2D::Bend2D(const Bend2D &right): deltaBeginEntry_m(right.deltaBeginEntry_m), deltaEndEntry_m(right.deltaEndEntry_m), polyOrderEntry_m(right.polyOrderEntry_m), - xExit_m(right.xExit_m), - zExit_m(right.zExit_m), deltaBeginExit_m(right.deltaBeginExit_m), deltaEndExit_m(right.deltaEndExit_m), polyOrderExit_m(right.polyOrderExit_m), @@ -76,6 +76,12 @@ Bend2D::Bend2D(const Bend2D &right): sinEntranceAngle_m(right.sinEntranceAngle_m), tanEntranceAngle_m(right.tanEntranceAngle_m), tanExitAngle_m(right.tanExitAngle_m), + beginToEnd_m(right.beginToEnd_m), + beginToEnd_lcs_m(right.beginToEnd_lcs_m), + toEntranceRegion_m(right.toEntranceRegion_m), + toExitRegion_m(right.toExitRegion_m), + computeAngleTrafo_m(right.computeAngleTrafo_m), + maxAngle_m(right.maxAngle_m), nSlices_m(right.nSlices_m){ setElType(isDipole); @@ -87,14 +93,14 @@ Bend2D::Bend2D(const std::string &name): messageHeader_m(" * "), pusher_m(), fieldmap_m(NULL), - fast_m(false), designRadius_m(0.0), exitAngle_m(0.0), fieldIndex_m(0.0), startField_m(0.0), endField_m(0.0), + widthEntranceFringe_m(0.0), + widthExitFringe_m(0.0), reinitialize_m(false), - recalcRefTraj_m(false), entranceParameter1_m(0.0), entranceParameter2_m(0.0), entranceParameter3_m(0.0), @@ -108,8 +114,6 @@ Bend2D::Bend2D(const std::string &name): deltaBeginEntry_m(0.0), deltaEndEntry_m(0.0), polyOrderEntry_m(0), - xExit_m(0.0), - zExit_m(0.0), deltaBeginExit_m(0.0), deltaEndExit_m(0.0), polyOrderExit_m(0), @@ -117,6 +121,7 @@ Bend2D::Bend2D(const std::string &name): sinEntranceAngle_m(0.0), tanEntranceAngle_m(0.0), tanExitAngle_m(0.0), + maxAngle_m(0.0), nSlices_m(1){ setElType(isDipole); @@ -133,13 +138,9 @@ Bend2D::~Bend2D() { } delete[] entryFieldValues_m; delete[] exitFieldValues_m; - entryFieldValues_m = NULL; - exitFieldValues_m = NULL; gsl_interp_accel_free(entryFieldAccel_m); gsl_interp_accel_free(exitFieldAccel_m); - entryFieldAccel_m = NULL; - exitFieldAccel_m = NULL; } } /* @@ -155,24 +156,8 @@ bool Bend2D::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) { - - if(designRadius_m > 0.0) { - - // Shift position to magnet frame. - Vector_t X = RefPartBunch_m->R[i]; - - Vector_t bField(0.0); - if (!calculateMapField(X, bField)) { - - B += fieldAmplitude_m * bField; - - return false; - } - - return true; - } - - return false; + Vector_t P; + return apply(RefPartBunch_m->R[i], P, t, E, B); } bool Bend2D::apply(const Vector_t &R, @@ -203,21 +188,7 @@ bool Bend2D::applyToReferenceParticle(const Vector_t &R, const double &t, Vector_t &E, Vector_t &B) { - if(designRadius_m > 0.0) { - - // Get field from field map. - Vector_t bField(0.0); - Vector_t X = R;// + Vector_t(0, 0, startField_m/* - elementEdge_m*/); - if (!calculateMapField(X, bField)) { - - B += fieldAmplitude_m * bField; - - return false; - } - - return true; - } - return false; + return apply(R, P, t, E, B); } void Bend2D::goOnline(const double &) { @@ -244,7 +215,6 @@ void Bend2D::initialise(PartBunchBase<double, 3> *bunch, double bendAngleX = 0.0; double bendAngleY = 0.0; calculateRefTrajectory(bendAngleX, bendAngleY); - recalcRefTraj_m = true; print(msg, bendAngleX, bendAngleY); // Pass start and end of field to calling function. @@ -816,7 +786,7 @@ bool Bend2D::findIdealBendParameters(double chordLength) { if(angle_m < 0.0) { // Negative angle is a positive bend rotated 180 degrees. entranceAngle_m = copysign(1, angle_m) * entranceAngle_m; - exitAngle_m = copysign(1, angle_m) * exitAngle_m; + setExitAngle(copysign(1, angle_m) * exitAngle_m); angle_m = std::abs(angle_m); rotationZAxis_m += Physics::pi; } @@ -845,27 +815,6 @@ bool Bend2D::findIdealBendParameters(double chordLength) { return reinitialize; } -void Bend2D::findReferenceExitOrigin(double &x, double &z) { - - /* - * Find x,z coordinates of reference trajectory as it passes exit edge - * of the bend magnet. This assumes an entrance position of (x,z) = (0,0). - */ - if(angle_m <= Physics::pi / 2.0) { - x = - designRadius_m * (1.0 - std::cos(angle_m)); - z = designRadius_m * std::sin(angle_m); - } else if(angle_m <= Physics::pi) { - x = -designRadius_m * (1.0 + std::sin(angle_m - Physics::pi / 2.0)); - z = designRadius_m * std::cos(angle_m - Physics::pi / 2.0); - } else if(angle_m <= 3.0 * Physics::pi / 2.0) { - x = -designRadius_m * (2.0 - std::cos(angle_m - Physics::pi)); - z = -designRadius_m * std::sin(angle_m - Physics::pi); - } else { - x = -designRadius_m * (1.0 - std::cos(angle_m - 3.0 * Physics::pi / 2.0)); - z = -designRadius_m * std::sin(angle_m - 3.0 * Physics::pi / 2.0); - } -} - bool Bend2D::initializeFieldMap(Inform &msg) { fieldmap_m = Fieldmap::getFieldmap(fileName_m, fast_m); @@ -1110,22 +1059,6 @@ void Bend2D::readFieldMap(Inform &msg) { } } -bool Bend2D::reinitialize() { - - setBendStrength(); - double bendAngleX = 0.0; - double bendAngleY = 0.0; - calculateRefTrajectory(bendAngleX, bendAngleY); - - // INFOMSG(level2 << "Bend design energy changed to: " - // << designEnergy_m * 1.0e-6 - // << " MeV" - // << endl); - print(*Ippl::Info, bendAngleX, bendAngleY); - - return false; -} - void Bend2D::setBendEffectiveLength(double startField, double endField) { // Find initial angle. @@ -1183,15 +1116,13 @@ void Bend2D::setFieldCalcParam() { tanEntranceAngle_m = tan(entranceAngle_m); deltaBeginEntry_m = std::abs(entranceParameter1_m - entranceParameter2_m); - deltaEndEntry_m = std::abs(entranceParameter2_m - entranceParameter3_m); + deltaEndEntry_m = std::abs(entranceParameter2_m - entranceParameter3_m); deltaBeginExit_m = std::abs(exitParameter1_m - exitParameter2_m); - deltaEndExit_m = std::abs(exitParameter2_m - exitParameter3_m); + deltaEndExit_m = std::abs(exitParameter2_m - exitParameter3_m); double bendAngle = getBendAngle(); double entranceAngle = getEntranceAngle(); - double exitAngle = getExitAngle(); - tanExitAngle_m = tan(exitAngle); double rotationAngleAboutZ = getRotationAboutZ(); Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ), @@ -1200,8 +1131,8 @@ void Bend2D::setFieldCalcParam() { Vector_t rotationAxis(0, -1, 0); Quaternion_t halfRotationAboutAxis(cos(0.5 * (0.5 * bendAngle - entranceAngle)), sin(0.5 * (0.5 * bendAngle - entranceAngle)) * rotationAxis); - Quaternion_t exitFaceRotation(cos(0.5 * (bendAngle - entranceAngle - exitAngle)), - sin(0.5 * (bendAngle - entranceAngle - exitAngle)) * rotationAxis); + Quaternion_t exitFaceRotation(cos(0.5 * (bendAngle - entranceAngle - exitAngle_m)), + sin(0.5 * (bendAngle - entranceAngle - exitAngle_m)) * rotationAxis); Vector_t chord = getChordLength() * halfRotationAboutAxis.rotate(Vector_t(0, 0, 1)); beginToEnd_lcs_m = CoordinateSystemTrafo(chord, exitFaceRotation.conjugate()); beginToEnd_m = beginToEnd_lcs_m * CoordinateSystemTrafo(Vector_t(0.0), rotationAboutZ.conjugate()); @@ -1295,7 +1226,6 @@ bool Bend2D::setupBendGeometry(Inform &msg, double &startField, double &endField if (getType() == RBEND) { setEntranceAngle(getEntranceAngle()); } - findReferenceExitOrigin(xExit_m, zExit_m); /* * Set field map geometry. @@ -1797,7 +1727,5 @@ std::array<double,2> Bend2D::getExitFringeFieldLength() const { extFFL[0] = ( exitParameter3_m-exitParameter2_m ); //after edge extFFL[1] = ( exitParameter2_m-exitParameter1_m ); //before edge - - return extFFL; } \ No newline at end of file diff --git a/src/Classic/AbsBeamline/Bend2D.h b/src/Classic/AbsBeamline/Bend2D.h index aa87c162806ad90b7e80a1f76e5e1bb4970b78db..822a361ed4203bd8dbd678b13dd25dad5408ad46 100644 --- a/src/Classic/AbsBeamline/Bend2D.h +++ b/src/Classic/AbsBeamline/Bend2D.h @@ -33,6 +33,7 @@ #endif #include <array> +#include <string> #include <vector> class Fieldmap; @@ -110,9 +111,6 @@ public: /// Set quadrupole field component. void setK1(double k1); - void resetReinitializeFlag(); - void resetRecalcRefTrajFlag(); - void setExitAngle(double exitAngle); virtual double getExitAngle() const; @@ -185,7 +183,6 @@ private: virtual bool findChordLength(Inform &msg, double &chordLength) = 0; bool findIdealBendParameters(double chordLength); - void findReferenceExitOrigin(double &x, double &z); bool initializeFieldMap(Inform &msg); bool inMagnetCentralRegion(const Vector_t &R) const; bool inMagnetEntranceRegion(const Vector_t &R) const; @@ -194,7 +191,6 @@ private: bool isPositionInExitField(const Vector_t &R) const; void print(Inform &msg, double bendAngleX, double bendAngle); void readFieldMap(Inform &msg); - bool reinitialize(); void setBendEffectiveLength(double startField, double endField); void setBendStrength(); void setEngeOriginDelta(double delta); @@ -217,7 +213,7 @@ private: /// through the bend. Fieldmap *fieldmap_m; /// Magnet field map. - bool fast_m; /// Flag to turn on fast field calculation. + const bool fast_m = false; /// Flag to turn on fast field calculation. double designRadius_m; /// Bend design radius (m). @@ -239,8 +235,6 @@ private: */ bool reinitialize_m; - bool recalcRefTraj_m; /// Re-calculate reference trajectory. - /* * Enge function field map members. */ @@ -278,12 +272,6 @@ private: /// function origin that Enge function ends. int polyOrderEntry_m; /// Enge function order for entry region. - /* - * The ideal reference trajectory passes through (xExit_m, zExit_m). - */ - double xExit_m; - double zExit_m; - double deltaBeginExit_m; /// Perpendicular distance from exit Enge /// function origin that Enge function starts. double deltaEndExit_m; /// Perpendicular distance from exit Enge @@ -347,16 +335,6 @@ void Bend2D::setK1(double k1) { fieldIndex_m = k1; } -inline -void Bend2D::resetReinitializeFlag() { - reinitialize_m = true; -} - -inline -void Bend2D::resetRecalcRefTrajFlag() { - recalcRefTraj_m = true; -} - inline void Bend2D::setMessageHeader(const std::string & header) { @@ -385,6 +363,7 @@ inline void Bend2D::setExitAngle(double angle) { exitAngle_m = angle; + tanExitAngle_m = tan(exitAngle_m); } inline diff --git a/src/Classic/AbsBeamline/Cyclotron.cpp b/src/Classic/AbsBeamline/Cyclotron.cpp index 7d018d14502bf95cda640830bc13fac1deb283dc..df86975637e58890c5e4182058c940408dea20a8 100644 --- a/src/Classic/AbsBeamline/Cyclotron.cpp +++ b/src/Classic/AbsBeamline/Cyclotron.cpp @@ -1037,10 +1037,10 @@ void Cyclotron::initialise(PartBunchBase<double, 3> *bunch, const int &fieldflag void Cyclotron::getFieldFromFile_FFA(const double &scaleFactor) { /* - Field is read in from ascci file (COSY output) in the oder: + Field is read in from ascii file (COSY output) in the order: R(m) theta(Deg) x(m) y(m) Bz(T). - Theta is the fast varing variable + Theta is the fast varying variable 2.0000 0.0 2.0000 0.0000 0.0000000000000000 2.0000 1.0 1.9997 0.0349 0.0000000000000000 diff --git a/src/Classic/AbsBeamline/Ring.cpp b/src/Classic/AbsBeamline/Ring.cpp index d2478854037af880ed1262a42d1334276713d7b0..eece8fb5fd7d8e633249179b356d775520e648a3 100644 --- a/src/Classic/AbsBeamline/Ring.cpp +++ b/src/Classic/AbsBeamline/Ring.cpp @@ -275,7 +275,7 @@ void Ring::appendElement(const Component &element) { << section->getStartNormal()(0) << ", " << section->getStartNormal()(1) << ", " << section->getStartNormal()(2) << "), phi " << dphi << endl; - msg << "* End position (" + msg << "* End position (" << section->getEndPosition()(0) << ", " << section->getEndPosition()(1) << ", " << section->getEndPosition()(2) << ") normal (" diff --git a/src/Classic/AbsBeamline/ScalingFFAMagnet.cpp b/src/Classic/AbsBeamline/ScalingFFAMagnet.cpp index dbdaeb3e1f11aea50307f3617270254225576aa0..9bf76ef91ad077ad6fda6f62ba889dda49cdd858 100644 --- a/src/Classic/AbsBeamline/ScalingFFAMagnet.cpp +++ b/src/Classic/AbsBeamline/ScalingFFAMagnet.cpp @@ -46,9 +46,7 @@ ScalingFFAMagnet::ScalingFFAMagnet(const ScalingFFAMagnet &right) phiEnd_m(right.phiEnd_m), azimuthalExtent_m(right.azimuthalExtent_m), verticalExtent_m(right.verticalExtent_m), centre_m(right.centre_m), dfCoefficients_m(right.dfCoefficients_m) { - if (endField_m != NULL) { - delete endField_m; - } + delete endField_m; endField_m = right.endField_m->clone(); RefPartBunch_m = right.RefPartBunch_m; setElType(isDrift); @@ -57,9 +55,7 @@ ScalingFFAMagnet::ScalingFFAMagnet(const ScalingFFAMagnet &right) } ScalingFFAMagnet::~ScalingFFAMagnet() { - if (endField_m != NULL) { - delete endField_m; - } + delete endField_m; } ElementBase* ScalingFFAMagnet::clone() const { diff --git a/src/Classic/Algorithms/PartBunch.cpp b/src/Classic/Algorithms/PartBunch.cpp index aff6a5bf258b6306989998bddaf7cdda4272935d..9612695d365c1fc0e9325a961d46e859a42f0431 100644 --- a/src/Classic/Algorithms/PartBunch.cpp +++ b/src/Classic/Algorithms/PartBunch.cpp @@ -33,7 +33,6 @@ #include "Distribution/Distribution.h" // OPAL file #include "Structure/FieldSolver.h" // OPAL file #include "Utilities/GeneralClassicException.h" -#include "Structure/LossDataSink.h" #include "Algorithms/ListElem.h" diff --git a/src/Classic/Algorithms/PartBunchBase.h b/src/Classic/Algorithms/PartBunchBase.h index d60377d1d990d45b1df15b023ef4d4c9464c9f7c..fb1f57be2e40483a8e828f9eacabbec69c0c74d2 100644 --- a/src/Classic/Algorithms/PartBunchBase.h +++ b/src/Classic/Algorithms/PartBunchBase.h @@ -149,7 +149,7 @@ public: /** delete particles which are too far away from the center of beam*/ void boundp_destroy(); - /** This is only temporary in order to get the collimator and pepperpot workinh */ + /** This is only temporary in order to get the collimator and pepperpot working */ size_t boundp_destroyT(); size_t destroyT(); @@ -574,7 +574,7 @@ protected: double t_m; /// mean energy of the bunch (MeV) double eKin_m; - /// energy spread of the beam in keV + /// energy spread of the beam in MeV double dE_m; /// the position along design trajectory double spos_m; diff --git a/src/Classic/Algorithms/PartBunchBase.hpp b/src/Classic/Algorithms/PartBunchBase.hpp index 1078afed21cb48affe09e90168282e1e8995bbfa..3df9b5d0bc5beecc0234897ff3928ada7a0cbece 100644 --- a/src/Classic/Algorithms/PartBunchBase.hpp +++ b/src/Classic/Algorithms/PartBunchBase.hpp @@ -1278,7 +1278,6 @@ void PartBunchBase<T, Dim>::get_PBounds(Vector_t &min, Vector_t &max) const { template <class T, unsigned Dim> void PartBunchBase<T, Dim>::calcBeamParameters() { - using Physics::c; Vector_t eps2, fac, rsqsum, psqsum, rpsum; @@ -1330,10 +1329,10 @@ void PartBunchBase<T, Dim>::calcBeamParameters() { rprms_m = rpsum * fac; - Dx_m = moments_m(0, 5) / N; + Dx_m = moments_m(0, 5) / N; DDx_m = moments_m(1, 5) / N; - Dy_m = moments_m(2, 5) / N; + Dy_m = moments_m(2, 5) / N; DDy_m = moments_m(3, 5) / N; // Find unnormalized emittance. @@ -1717,7 +1716,6 @@ double PartBunchBase<T, Dim>::calcMeanPhi() { return meanPhi; } - // this function reset the BinID for each particles according to its current beta*gamma // it is for multi-turn extraction cyclotron with small energy gain // the bin number can be different with the bunch number @@ -1901,7 +1899,7 @@ template <class T, unsigned Dim> void PartBunchBase<T, Dim>::calcEMean() { const double totalNp = static_cast<double>(getTotalNum()); - const double locNp = static_cast<double>(getLocalNum()); + const double locNp = static_cast<double>(getLocalNum()); eKin_m = 0.0; diff --git a/src/Classic/Algorithms/Vektor.h b/src/Classic/Algorithms/Vektor.h index 87fd5d29b8057a87c494b938fc7c29075f5a639f..fe755ff0f8987866a9916606361f9019c25106b0 100644 --- a/src/Classic/Algorithms/Vektor.h +++ b/src/Classic/Algorithms/Vektor.h @@ -1,8 +1,6 @@ #ifndef OPAL_VEKTOR_HH #define OPAL_VEKTOR_HH -#include <limits> - #include "AppTypes/Vektor.h" typedef Vektor<double, 3> Vector_t; diff --git a/src/Classic/Fields/Interpolation/Interpolator3dGridTo3d.cpp b/src/Classic/Fields/Interpolation/Interpolator3dGridTo3d.cpp index 3e4ee2d863498c8dbea3ceb9bf2e0f49044d41c5..1553606ffeb672e1e49955b241429da46d195dca 100644 --- a/src/Classic/Fields/Interpolation/Interpolator3dGridTo3d.cpp +++ b/src/Classic/Fields/Interpolation/Interpolator3dGridTo3d.cpp @@ -73,7 +73,7 @@ void Interpolator3dGridTo3d::setAll(ThreeDGrid* grid, break; default: throw(LogicalError( - "Interpolator3dGridTo3d::Set", + "Interpolator3dGridTo3d::setAll", "Did not recognise interpolation algorithm" ) ); diff --git a/src/Classic/Fields/Interpolation/TriLinearInterpolator.cpp b/src/Classic/Fields/Interpolation/TriLinearInterpolator.cpp index f89764e88fc7e5f46f8995021caef6152038b7d8..95d4febfa6096dc6342aa7be77cebd876336f741 100644 --- a/src/Classic/Fields/Interpolation/TriLinearInterpolator.cpp +++ b/src/Classic/Fields/Interpolation/TriLinearInterpolator.cpp @@ -64,9 +64,9 @@ void TriLinearInterpolator::function double dx = (Point[0]-coordinates_m->x(i+1))/ (coordinates_m->x(i+2)-coordinates_m->x(i+1)); double f_x[2][2]; - f_x[0][0] = (f_m[i+1][j][k] - f_m[i][j][k]) * dx + f_m[i][j][k]; - f_x[1][0] = (f_m[i+1][j+1][k] - f_m[i][j+1][k]) * dx + f_m[i][j+1][k]; - f_x[0][1] = (f_m[i+1][j][k+1] - f_m[i][j][k+1]) * dx + f_m[i][j][k+1]; + f_x[0][0] = (f_m[i+1][j] [k] - f_m[i][j] [k]) * dx + f_m[i][j] [k]; + f_x[1][0] = (f_m[i+1][j+1][k] - f_m[i][j+1][k]) * dx + f_m[i][j+1][k]; + f_x[0][1] = (f_m[i+1][j] [k+1] - f_m[i][j] [k+1]) * dx + f_m[i][j] [k+1]; f_x[1][1] = (f_m[i+1][j+1][k+1] - f_m[i][j+1][k+1]) * dx + f_m[i][j+1][k+1]; // interpolation in y double f_xy[2]; diff --git a/src/Classic/Fields/SectorMagneticFieldMap.cpp b/src/Classic/Fields/SectorMagneticFieldMap.cpp index c518a3a65720b67b1dfb77ba5447fc6bf6054dce..28c6fafb4d0f8d89db919685a1be403b139822a3 100644 --- a/src/Classic/Fields/SectorMagneticFieldMap.cpp +++ b/src/Classic/Fields/SectorMagneticFieldMap.cpp @@ -337,7 +337,7 @@ VectorMap* SectorMagneticFieldMap::IO::getInterpolatorPolyPatch( // too lazy to write code to handle this case - not available to user anyway if (sym != SectorMagneticFieldMap::dipole) { throw(LogicalError( - "SectorMagneticFieldMap::IO::ReadMap", + "SectorMagneticFieldMap::IO::getInterpolatorPolyPatch", "Failed to recognise symmetry type" )); } diff --git a/src/Classic/Structure/MeshGenerator.cpp b/src/Classic/Structure/MeshGenerator.cpp index bf70f5eea17414b794c371886137745e7dca0ce8..5bc1f4573a1f00a5e83730381aa4907493b02a91 100644 --- a/src/Classic/Structure/MeshGenerator.cpp +++ b/src/Classic/Structure/MeshGenerator.cpp @@ -15,7 +15,6 @@ MeshGenerator::MeshGenerator(): { } void MeshGenerator::add(const ElementBase &element) { - auto apert = element.getAperture(); double start = 0.0; MeshData mesh; @@ -35,6 +34,7 @@ void MeshGenerator::add(const ElementBase &element) { double end, length; element.getElementDimensions(start, end); length = end - start; + auto apert = element.getAperture(); switch (apert.first) { case ElementBase::RECTANGULAR: diff --git a/src/Classic/Structure/PeakFinder.cpp b/src/Classic/Structure/PeakFinder.cpp index e7be6f2f70515d9a27d801c4dd5ece4a2ff05d01..5c6150269376d0a9103ad44440b32b5ea0f0d755 100644 --- a/src/Classic/Structure/PeakFinder.cpp +++ b/src/Classic/Structure/PeakFinder.cpp @@ -29,39 +29,39 @@ PeakFinder::PeakFinder(std::string elem, double min, void PeakFinder::addParticle(const Vector_t& R) { - + double radius = std::hypot(R(0),R(1)); radius_m.push_back(radius); - + peakRadius_m += radius; ++registered_m; } void PeakFinder::evaluate(const int& turn) { - + if ( first_m ) { turn_m = turn; first_m = false; } - + if ( turn_m != turn ) { finished_m = true; } - + bool globFinished = false; - + if ( !singlemode_m ) allreduce(finished_m, globFinished, 1, std::logical_and<bool>()); else globFinished = finished_m; - + if ( globFinished ) { - + this->computeCentroid_m(); - + turn_m = turn; - + // reset peakRadius_m = 0.0; registered_m = 0; @@ -71,29 +71,29 @@ void PeakFinder::evaluate(const int& turn) { void PeakFinder::save() { - + createHistogram_m(); - + // last turn is not yet computed this->computeCentroid_m(); - + if ( !peaks_m.empty() ) { // only rank 0 will go in here - + fn_m = element_m + std::string(".peaks"); hist_m = element_m + std::string(".hist"); - + INFOMSG("Save " << fn_m << " and " << hist_m << endl); - + if(OpalData::getInstance()->inRestartRun()) this->append_m(); else this->open_m(); - + this->saveASCII_m(); - + this->close_m(); - + } radius_m.clear(); @@ -104,16 +104,16 @@ void PeakFinder::save() { void PeakFinder::computeCentroid_m() { double globPeakRadius = 0.0; int globRegister = 0; - + //FIXME inefficient if ( !singlemode_m ) { reduce(peakRadius_m, globPeakRadius, 1, std::plus<double>()); - reduce(registered_m, globRegister, 1, std::plus<int>()); + reduce(registered_m, globRegister, 1, std::plus<int>()); } else { globPeakRadius = peakRadius_m; globRegister = registered_m; } - + if ( Ippl::myNode() == 0 ) { if ( globRegister > 0 ) peaks_m.push_back(globPeakRadius / double(globRegister)); @@ -124,7 +124,7 @@ void PeakFinder::createHistogram_m() { /* * create local histograms */ - + globHist_m.resize(nBins_m); container_t locHist(nBins_m,0.0); @@ -134,7 +134,7 @@ void PeakFinder::createHistogram_m() { if (bin < 0 || (unsigned int)bin >= nBins_m) continue; // Probe might save particles outside its boundary ++locHist[bin]; } - + /* * create global histograms */ @@ -144,7 +144,7 @@ void PeakFinder::createHistogram_m() { } else { globHist_m.swap(locHist); } - + // reduce(locHist.data(), globHist_m.data(), locHist.size(), std::plus<double>()); } diff --git a/src/Classic/Structure/PeakFinder.h b/src/Classic/Structure/PeakFinder.h index 2fbf8e5c342721f86ee224f56a0ef11ab3a57990..a6b0777e8c4318f3defb3ba1f00db714318ecec6 100644 --- a/src/Classic/Structure/PeakFinder.h +++ b/src/Classic/Structure/PeakFinder.h @@ -23,34 +23,34 @@ #include <list> class PeakFinder { - + public: using container_t = std::vector<double>; - + public: - + PeakFinder() = delete; PeakFinder(std::string elem, double min, double max, double binwidth, bool singlemode); - + /*! * Append the particle coordinates to the container * @param R is a particle position (x, y, z) */ void addParticle(const Vector_t& R); - + /*! * Evaluate the centroid of a turn. */ void evaluate(const int& turn); - + void save(); - + private: - + // compute global histogram, involves some inter-node communication void createHistogram_m(); - + /*************** * Output file * ***************/ @@ -62,29 +62,29 @@ private: void close_m(); /// Write to output file void saveASCII_m(); - + void computeCentroid_m(); private: container_t radius_m; /// global histogram values container_t globHist_m; - + /// filename with extension (.peaks) std::string fn_m; - + /// histogram filename with extension (.hist) std::string hist_m; /// used to write out the data std::ofstream os_m; - + /// used to write out the histrogram std::ofstream hos_m; - + /// Element/probe name, for name output file std::string element_m; - + // Histogram details /// Number of bins unsigned int nBins_m; @@ -92,7 +92,7 @@ private: double binWidth_m; ///@{ histogram size double min_m, max_m; - + int turn_m; double peakRadius_m; int registered_m; diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp index 4393c8246ba264ca76cc15bfa472a45f827533bb..efee18c1018fff4889f734db0e61ea2e3d5f2ac2 100644 --- a/src/Distribution/Distribution.cpp +++ b/src/Distribution/Distribution.cpp @@ -1377,7 +1377,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub std::array<double,3> Emit = siggen->getEmittances(); if (rguess > 0) - *gmsg << "* RGUESS " << rguess/1000.0 << " (m) " << endl; + *gmsg << "* RGUESS " << rguess << " (m) " << endl; *gmsg << "* Converged (Ex, Ey, Ez) = (" << Emit[0] << ", " << Emit[1] << ", " << Emit[2] << ") pi mm mrad for E= " << E_m*1E-6 << " (MeV)" << endl; @@ -1552,7 +1552,7 @@ void Distribution::createOpalCycl(PartBunchBase<double, 3> *beam, // Set distribution type. setDistType(); - // Emitting particles in not supported in OpalCyclT. + // Emitting particles is not supported in OPAL-cycl. emitting_m = false; // Create distribution. @@ -1911,12 +1911,12 @@ size_t Distribution::emitParticles(PartBunchBase<double, 3> *beam, double eZ) { * is much slower than doing a swap and popping off the end * of the vector. */ - std::swap(xDist_m.at(*ptbErasedIt), xDist_m.back()); - std::swap(pxDist_m.at(*ptbErasedIt), pxDist_m.back()); - std::swap(yDist_m.at(*ptbErasedIt), yDist_m.back()); - std::swap(pyDist_m.at(*ptbErasedIt), pyDist_m.back()); + std::swap( xDist_m.at(*ptbErasedIt), xDist_m.back()); + std::swap(pxDist_m.at(*ptbErasedIt), pxDist_m.back()); + std::swap( yDist_m.at(*ptbErasedIt), yDist_m.back()); + std::swap(pyDist_m.at(*ptbErasedIt), pyDist_m.back()); std::swap(tOrZDist_m.at(*ptbErasedIt), tOrZDist_m.back()); - std::swap(pzDist_m.at(*ptbErasedIt), pzDist_m.back()); + std::swap(pzDist_m.at(*ptbErasedIt), pzDist_m.back()); if (additionalRNs_m.size() == xDist_m.size()) { std::swap(additionalRNs_m.at(*ptbErasedIt), additionalRNs_m.back()); @@ -2488,7 +2488,7 @@ void Distribution::generateGaussZ(size_t numberOfParticles) { } #endif /* - //Sets the GSL error handler off, exception will be handled internaly with a renormalization method + //Sets the GSL error handler off, exception will be handled internally with a renormalization method gsl_set_error_handler_off(); */ int errcode = gsl_linalg_cholesky_decomp(corMat); @@ -4497,34 +4497,27 @@ void Distribution::writeOutFileHeader() { } else { outputFile.setf(std::ios::left); outputFile << "# "; + outputFile.width(17); + outputFile << "x [m]"; + outputFile.width(17); + outputFile << "px [betax gamma]"; + outputFile.width(17); + outputFile << "y [m]"; + outputFile.width(17); + outputFile << "py [betay gamma]"; if (emitting_m) { - outputFile.width(17); - outputFile << "x [m]"; - outputFile.width(17); - outputFile << "px [betax gamma]"; - outputFile.width(17); - outputFile << "y [m]"; - outputFile.width(17); - outputFile << "py [betay gamma]"; outputFile.width(17); outputFile << "t [s]"; - outputFile.width(17); - outputFile << "pz [betaz gamma]" ; - outputFile.width(17); - outputFile << "Bin Number" << std::endl; } else { - outputFile.width(17); - outputFile << "x [m]"; - outputFile.width(17); - outputFile << "px [betax gamma]"; - outputFile.width(17); - outputFile << "y [m]"; - outputFile.width(17); - outputFile << "py [betay gamma]"; outputFile.width(17); outputFile << "z [m]"; + } + outputFile.width(17); + outputFile << "pz [betaz gamma]" ; + if (emitting_m) { outputFile.width(17); - outputFile << "pz [betaz gamma]"; + outputFile << "Bin Number" << std::endl; + } else { if (numberOfEnergyBins_m > 0) { outputFile.width(17); outputFile << "Bin Number"; @@ -4737,10 +4730,10 @@ void Distribution::adjustPhaseSpace(double massIneV) { } double avrg[6]; - avrg[0] = std::accumulate(xDist_m.begin(), xDist_m.end(), 0.0) / totalNumberParticles_m; - avrg[1] = std::accumulate(pxDist_m.begin(), pxDist_m.end(), 0.0) / totalNumberParticles_m; - avrg[2] = std::accumulate(yDist_m.begin(), yDist_m.end(), 0.0) / totalNumberParticles_m; - avrg[3] = std::accumulate(pyDist_m.begin(), pyDist_m.end(), 0.0) / totalNumberParticles_m; + avrg[0] = std::accumulate( xDist_m.begin(), xDist_m.end(), 0.0) / totalNumberParticles_m; + avrg[1] = std::accumulate(pxDist_m.begin(), pxDist_m.end(), 0.0) / totalNumberParticles_m; + avrg[2] = std::accumulate( yDist_m.begin(), yDist_m.end(), 0.0) / totalNumberParticles_m; + avrg[3] = std::accumulate(pyDist_m.begin(), pyDist_m.end(), 0.0) / totalNumberParticles_m; avrg[4] = std::accumulate(tOrZDist_m.begin(), tOrZDist_m.end(), 0.0) / totalNumberParticles_m; avrg[5] = 0.0; for (unsigned int i = 0; i < pzDist_m.size(); ++ i) { @@ -4756,10 +4749,10 @@ void Distribution::adjustPhaseSpace(double massIneV) { // \sum_{i = 0}^{N-1} \sqrt{(pz_i + \eps)^2 + px_i^2 + py_i^2} = N p double eps = avrgpz_m - avrg[5]; for (unsigned int i = 0; i < pzDist_m.size(); ++ i) { - xDist_m[i] -= avrg[0]; - pxDist_m[i] -= avrg[1]; - yDist_m[i] -= avrg[2]; - pyDist_m[i] -= avrg[3]; + xDist_m[i] -= avrg[0]; + pxDist_m[i] -= avrg[1]; + yDist_m[i] -= avrg[2]; + pyDist_m[i] -= avrg[3]; tOrZDist_m[i] -= avrg[4]; pzDist_m[i] += eps; } diff --git a/src/Elements/OpalBeamline.cpp b/src/Elements/OpalBeamline.cpp index 785db4170a0308bae99c69839f5c7ea832d2884d..24f3c4c4206d033850740c87aefddb7a33207e29 100644 --- a/src/Elements/OpalBeamline.cpp +++ b/src/Elements/OpalBeamline.cpp @@ -1,5 +1,5 @@ #include "Elements/OpalBeamline.h" -#include "Utilities/OpalException.h" + #include "Utilities/Options.h" #include "Utilities/Util.h" #include "AbstractObjects/OpalData.h" @@ -9,7 +9,6 @@ #include <boost/filesystem.hpp> #include <boost/regex.hpp> #include <fstream> -using namespace std; OpalBeamline::OpalBeamline(): elements_m(), @@ -47,10 +46,6 @@ std::set<std::shared_ptr<Component>> OpalBeamline::getElements(const Vector_t &x return elementSet; } -void OpalBeamline::getKFactors(const unsigned int &index, const Vector_t &pos, const long &sindex, const double &t, Vector_t &KR, Vector_t &KT) { - -} - unsigned long OpalBeamline::getFieldAt(const unsigned int &index, const Vector_t &pos, const long &sindex, const double &t, Vector_t &E, Vector_t &B) { unsigned long rtv = 0x00; @@ -131,16 +126,6 @@ void OpalBeamline::switchElements(const double &min, const double &max, const do } } - -void OpalBeamline::switchAllElements() { - for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) { - if(!(*flit).isOn()) { - (*flit).setOn(-1.0); - } - } - -} - void OpalBeamline::switchElementsOff(const double &min, ElementBase::ElementType eltype) { if(eltype == ElementBase::ANY) { for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) { @@ -174,11 +159,6 @@ void OpalBeamline::prepareSections() { void OpalBeamline::print(Inform &msg) const { } - -double OpalBeamline::calcBeamlineLength() { - return 0.0; -} - void OpalBeamline::swap(OpalBeamline & rhs) { std::swap(elements_m, rhs.elements_m); std::swap(prepared_m, rhs.prepared_m); @@ -210,7 +190,6 @@ FieldList OpalBeamline::getElementByType(ElementBase::ElementType type) { void OpalBeamline::compute3DLattice() { static unsigned int order = 0; - FieldList::iterator it = elements_m.begin(); const FieldList::iterator end = elements_m.end(); unsigned int minOrder = order; @@ -228,7 +207,7 @@ void OpalBeamline::compute3DLattice() { return edgeA < edgeB; }); - it = elements_m.begin(); + FieldList::iterator it = elements_m.begin(); for (; it != end; ++ it) { if ((*it).isPositioned()) { continue; @@ -293,7 +272,7 @@ void OpalBeamline::compute3DLattice() { double endPriorPathLength = 0.0; CoordinateSystemTrafo currentCoordTrafo = coordTransformationTo_m; - it = elements_m.begin(); + FieldList::iterator it = elements_m.begin(); for (; it != end; ++ it) { if ((*it).isPositioned()) continue; @@ -366,102 +345,6 @@ void OpalBeamline::compute3DLattice() { elements_m.sort(ClassicField::SortAsc); } -void OpalBeamline::plot3DLattice() { - auto opal = OpalData::getInstance(); - if (opal->isOptimizerRun() || Ippl::myNode() != 0) return; - - elements_m.sort([](const ClassicField& a, const ClassicField& b) { - double edgeA = 0.0, edgeB = 0.0; - if (a.getElement()->isElementPositionSet()) - edgeA = a.getElement()->getElementPosition(); - - if (b.getElement()->isElementPositionSet()) - edgeB = b.getElement()->getElementPosition(); - - return edgeA < edgeB; - }); - - FieldList::iterator it = elements_m.begin(); - FieldList::iterator end = elements_m.end(); - - Quaternion rotDiagonal(0.5, 0.5 * Vector_t(-1, 1, -1)); - - Vector_t origin = rotDiagonal.rotate(coordTransformationTo_m.getOrigin()); - Vector_t direction = rotDiagonal.rotate(coordTransformationTo_m.rotateFrom(Vector_t(0, 0, 1))); - Vector_t minX = Vector_t(999999.9), maxX(-999999.9); - std::map<std::string, std::vector<Vector_t > > elementCorners; - - for (; it != end; ++ it) { - std::shared_ptr<Component> element = (*it).getElement(); - CoordinateSystemTrafo toBegin = (*it).getCoordTransformationTo(); - std::vector<Vector_t> corners; - - if (element->getType() == ElementBase::RBEND || element->getType() == ElementBase::SBEND) { - std::vector<Vector_t> outline = static_cast<const Bend2D*>(element.get())->getOutline(); - - for (auto point: outline) { - corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(point))); - } - } else { - CoordinateSystemTrafo toEnd = element->getEdgeToEnd() * toBegin; - auto aperture = element->getAperture(); - double elementHeightFront = aperture.second[0]; - double elementHeightBack = aperture.second[0] * aperture.second[2]; - - corners.push_back(rotDiagonal.rotate(toEnd.transformFrom(Vector_t(0.0)))); - corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(Vector_t(0)))); - corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(-elementHeightFront * Vector_t(1, 0, 0)))); - corners.push_back(rotDiagonal.rotate(toEnd.transformFrom(-elementHeightBack * Vector_t(1, 0, 0)))); - corners.push_back(rotDiagonal.rotate(toEnd.transformFrom(Vector_t(0.0)))); - corners.push_back(rotDiagonal.rotate(toEnd.transformFrom(elementHeightBack * Vector_t(1, 0, 0)))); - corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(elementHeightFront * Vector_t(1, 0, 0)))); - corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(Vector_t(0)))); - } - - elementCorners.insert(std::make_pair(element->getName(), corners)); - const unsigned int numCorners = corners.size(); - for (unsigned int i = 0 ; i < numCorners; ++ i) { - const Vector_t & X = corners[i]; - - if (X(0) < minX(0)) minX(0) = X(0); - else if (X(0) > maxX(0)) maxX(0) = X(0); - - if (X(1) < minX(1)) minX(1) = X(1); - else if (X(1) > maxX(1)) maxX(1) = X(1); - } - } - - std::ofstream gpl; - std::string fileName = "data/" + OpalData::getInstance()->getInputBasename() + "_ElementPositions.gpl"; - if (OpalData::getInstance()->getOpenMode() == OpalData::OPENMODE::APPEND && - boost::filesystem::exists(fileName)) { - gpl.open(fileName, std::ios_base::app); - } else { - gpl.open(fileName); - } - gpl.precision(8); - - it = elements_m.begin(); - for (; it != end; ++ it) { - std::shared_ptr<Component> element = (*it).getElement(); - - if (element->getType() != ElementBase::DRIFT) { - const std::vector<Vector_t> &corners = elementCorners[element->getName()]; - const unsigned int numCorners = corners.size(); - - gpl << "# " << element->getName() << "\n"; - for (unsigned int i = 0; i < numCorners; ++ i) { - gpl << std::setw(18) << corners[i](0) - << std::setw(18) << -corners[i](1) << "\n"; - } - gpl << std::setw(18) << corners.front()(0) - << std::setw(18) << -corners.front()(1) << "\n\n"; - } - } - - elements_m.sort(ClassicField::SortAsc); -} - void OpalBeamline::save3DLattice() { if (Ippl::myNode() != 0 || OpalData::getInstance()->isOptimizerRun()) return; @@ -496,8 +379,8 @@ void OpalBeamline::save3DLattice() { 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(); + toEnd = bendElement->getBeginToEnd_local() * (*it).getCoordTransformationTo(); + exit3D = toEnd.getOrigin(); unsigned int size = designPath.size(); unsigned int minNumSteps = std::max(20.0, @@ -519,12 +402,12 @@ void OpalBeamline::save3DLattice() { for (unsigned int i = frequency; i + 1 < size; i += frequency) { - Vector_t position = (*it).getCoordTransformationTo().transformFrom(designPath[i]); + position = (*it).getCoordTransformationTo().transformFrom(designPath[i]); pos << std::setw(30) << std::left << std::string("\"MID: ") + element->getName() + std::string("\"") << std::setw(18) << std::setprecision(10) << position(2) << std::setw(18) << std::setprecision(10) << position(0) << std::setw(18) << std::setprecision(10) << position(1) - << endl; + << std::endl; } position = (*it).getCoordTransformationTo().transformFrom(designPath.back()); @@ -599,7 +482,7 @@ namespace { unsigned int getMinimalSignificantDigits(double num, const unsigned int maxDigits) { char buf[32]; snprintf(buf, 32, "%.*f", maxDigits + 1, num); - string numStr(buf); + std::string numStr(buf); unsigned int length = numStr.length(); unsigned int numDigits = maxDigits; diff --git a/src/Elements/OpalBeamline.h b/src/Elements/OpalBeamline.h index ba24ab19678798e14069ac36d3d0412727203eab..6ef4455d236edb5b56c76f5480af820b39cd47e9 100644 --- a/src/Elements/OpalBeamline.h +++ b/src/Elements/OpalBeamline.h @@ -1,8 +1,6 @@ #ifndef OPAL_BEAMLINE_H #define OPAL_BEAMLINE_H -#include <list> -#include <limits> #include <set> #include <string> @@ -19,21 +17,17 @@ #include "AbsBeamline/Septum.h" #include "AbsBeamline/Source.h" -#include "BasicActions/Option.h" -#include "Utilities/Options.h" #include "Utilities/ClassicField.h" #include "Algorithms/CoordinateSystemTrafo.h" -class Tracker; template <class T, unsigned Dim> class PartBunchBase; class ParticleMatterInteractionHandler; class BoundaryGeometry; -class WakeFunction; -#define BEAMLINE_EOL 0x80000000 // end of line -#define BEAMLINE_PARTICLEMATTERINTERACTION 0x08000000 // has particle matter interaction +//#define BEAMLINE_EOL 0x80000000 // end of line +//#define BEAMLINE_PARTICLEMATTERINTERACTION 0x08000000 // has particle matter interaction class OpalBeamline { @@ -66,19 +60,14 @@ public: double getEnd(const Vector_t &) const; void switchElements(const double &, const double &, const double &kineticEnergy, const bool &nomonitors = false); - void switchAllElements(); void switchElementsOff(const double &, ElementBase::ElementType eltype = ElementBase::ANY); void switchElementsOff(); - WakeFunction *getWakeFunction(const unsigned int &); - std::shared_ptr<const ElementBase> getWakeFunctionOwner(const unsigned int &); - ParticleMatterInteractionHandler *getParticleMatterInteractionHandler(const unsigned int &); BoundaryGeometry *getBoundaryGeometry(const unsigned int &); - void getKFactors(const unsigned int &index, const Vector_t &pos, const long &sindex, const double &t, Vector_t &KR, Vector_t &KT); unsigned long getFieldAt(const unsigned int &, const Vector_t &, const long &, const double &, Vector_t &, Vector_t &); unsigned long getFieldAt(const Vector_t &, const Vector_t &, const double &, Vector_t &, Vector_t &); @@ -87,18 +76,12 @@ public: void prepareSections(); void compute3DLattice(); - void plot3DLattice(); void save3DLattice(); void save3DInput(); void print(Inform &) const; FieldList getElementByType(ElementBase::ElementType); - // need this for autophasing in case we have multiple tracks - double calcBeamlineLength(); - - void removeElement(const std::string &ElName); - void swap(OpalBeamline & rhs); void merge(OpalBeamline &rhs); @@ -169,16 +152,6 @@ void OpalBeamline::visit<Septum>(const Septum &element, BeamlineVisitor &, PartB WARNMSG(element.getTypeString() << " not implemented yet!" << endl); } -inline -void OpalBeamline::removeElement(const std::string &ElName) { - for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) { - if(flit->getElement()->getName() == ElName) { - flit->setStart(-flit->getEnd()); - flit->setEnd(-flit->getEnd()); - } - } -} - inline Vector_t OpalBeamline::transformTo(const Vector_t &r) const { return coordTransformationTo_m.transformTo(r); diff --git a/src/Lines/Sequence.cpp b/src/Lines/Sequence.cpp index eb807065d5bd788dcfdc7932791f184ef89b0ca5..591a351a7d84d2106df6ad0e9b87cd2403e752bf 100644 --- a/src/Lines/Sequence.cpp +++ b/src/Lines/Sequence.cpp @@ -23,7 +23,6 @@ #include "BeamlineCore/DriftRep.h" #include "Beamlines/Beamline.h" #include "Elements/OpalDrift.h" -#include "Lines/FlatWriter.h" #include "Lines/Replacer.h" #include "Lines/SequenceParser.h" diff --git a/src/Optimize/OpalSimulation.cpp b/src/Optimize/OpalSimulation.cpp index c8d9f9ada774bd5e5edeaa8914cbb4f8f0af9d44..a614394170890811da0c72af5c67bf50afe78677 100644 --- a/src/Optimize/OpalSimulation.cpp +++ b/src/Optimize/OpalSimulation.cpp @@ -350,17 +350,8 @@ void OpalSimulation::run() { // setup OPAL command line options std::ostringstream inputFileName; inputFileName << simulationName_ << ".in"; - - char exe_name[] = "opal"; char *inputfile = new char[inputFileName.str().size()+1] ; strcpy(inputfile, inputFileName.str().c_str()); - char nocomm[] = "--nocomminit"; - char info[] = "--info"; - char info0[] = "0"; - char warn[] = "--warn"; - char warn0[] = "0"; - char *arg[] = { exe_name, inputfile, nocomm, info, info0, warn, warn0 }; - int seed = Options::seed; CmdArguments_t args = getArgs(); @@ -374,6 +365,14 @@ void OpalSimulation::run() { "Restart specified but no restart H5 file available."); } + char exe_name[] = "opal"; + char nocomm[] = "--nocomminit"; + char info[] = "--info"; + char info0[] = "0"; + char warn[] = "--warn"; + char warn0[] = "0"; + char *arg[] = { exe_name, inputfile, nocomm, info, info0, warn, warn0 }; + //FIXME: this seems to crash OPAL in some cases //redirectOutToFile(); #ifdef SUPRESS_OUTPUT diff --git a/src/Sample/SampleCmd.cpp b/src/Sample/SampleCmd.cpp index 4fa32f07a710885dc1a9d7b628f198f08bf668c6..01970fed664d9dce5ffe89f68b1109e5b43a971e 100644 --- a/src/Sample/SampleCmd.cpp +++ b/src/Sample/SampleCmd.cpp @@ -446,7 +446,7 @@ void SampleCmd::execute() { boost::shared_ptr<Comm_t> comm(new Comm_t(args, MPI_COMM_WORLD)); if (comm->isWorker()) stashEnvironment(); - + if ( comm->isOptimizer() ) { for (sampleMethods_t::iterator it = sampleMethods.begin(); it != sampleMethods.end(); ++it) diff --git a/src/Structure/Beam.cpp b/src/Structure/Beam.cpp index 8f9f54040a9dd0cafc6533cd8477f9df963c1549..856da90f664f80fc2a812956a0ca287145d0a9b7 100644 --- a/src/Structure/Beam.cpp +++ b/src/Structure/Beam.cpp @@ -268,7 +268,7 @@ void Beam::update() { if(gamma > 1.0) { reference.setGamma(gamma); } else { - throw OpalException("Beam::execute()", + throw OpalException("Beam::update()", "\"GAMMA\" should be greater than 1."); } } else if(itsAttr[ENERGY]) { @@ -276,7 +276,7 @@ void Beam::update() { if(energy > reference.getM()) { reference.setE(energy); } else { - throw OpalException("Beam::execute()", + throw OpalException("Beam::update()", "\"ENERGY\" should be greater than \"MASS\"."); } } else if(itsAttr[PC]) { @@ -284,7 +284,7 @@ void Beam::update() { if(pc > 0.0) { reference.setP(pc); } else { - throw OpalException("Beam::execute()", + throw OpalException("Beam::update()", "\"PC\" should be greater than 0."); } }; diff --git a/src/Structure/ElementPositionWriter.cpp b/src/Structure/ElementPositionWriter.cpp index 62d83c9bbc8b90618fca2664f5b2ef98d48b6c99..4cd4b9bfae44c786651a57c94e4a9bf68e4a4367 100644 --- a/src/Structure/ElementPositionWriter.cpp +++ b/src/Structure/ElementPositionWriter.cpp @@ -1,5 +1,6 @@ #include "ElementPositionWriter.h" -#include "Structure/LossDataSink.h" + +#include "Ippl.h" #include "AbstractObjects/OpalData.h" ElementPositionWriter::ElementPositionWriter(const std::string& fname) diff --git a/src/Structure/H5PartWrapper.cpp b/src/Structure/H5PartWrapper.cpp index a4bc79becc499f299deb17b78f1bc8cade8222d7..af3005f07187bbb35a6b183b5c24e8b0e5e489e5 100644 --- a/src/Structure/H5PartWrapper.cpp +++ b/src/Structure/H5PartWrapper.cpp @@ -13,8 +13,6 @@ #include <fstream> -extern Inform *gmsg; - namespace { const h5_int64_t H5TypesCHAR = H5_STRING_T; const h5_int64_t H5TypesFLOAT = H5_FLOAT32_T; diff --git a/src/Structure/H5PartWrapperForPC.cpp b/src/Structure/H5PartWrapperForPC.cpp index 14e5237250f09e5785d4c180abc38c1dc82bf975..abff19b4770ca52e41250d768d2bb83c7e2d0c1d 100644 --- a/src/Structure/H5PartWrapperForPC.cpp +++ b/src/Structure/H5PartWrapperForPC.cpp @@ -205,7 +205,7 @@ void H5PartWrapperForPC::readStepData(PartBunchBase<double, 3>* bunch, for(long int n = 0; n < numParticles; ++ n) { bunch->M[n] = f64buffer[n]; } - + if ( bunch->getNumBunch() > 1 ) { std::vector<char> ibuffer(numParticles * sizeof(h5_int64_t)); h5_int64_t *i64buffer = reinterpret_cast<h5_int64_t*>(&ibuffer[0]); @@ -459,7 +459,7 @@ void H5PartWrapperForPC::writeStepHeader(PartBunchBase<double, 3>* bunch, } catch (std::out_of_range & m) { ERRORMSG(m.what() << endl); - throw OpalException("H5PartWrapperForPC::wirteStepHeader", + throw OpalException("H5PartWrapperForPC::writeStepHeader", "some additional step attribute not found"); } } @@ -560,10 +560,10 @@ void H5PartWrapperForPC::writeStepData(PartBunchBase<double, 3>* bunch) { i64buffer[i] = bunch->Bin[i]; for (size_t i = skipID; i < numLocalParticles; ++ i) i64buffer[i] = bunch->Bin[i + 1]; - + WRITEDATA(Int64, file_m, "bin", i64buffer); } - + for (size_t i = 0; i < skipID; ++ i) i64buffer[i] = bunch->bunchNum[i]; for (size_t i = skipID; i < numLocalParticles; ++ i) diff --git a/src/Structure/H5PartWrapperForPC.h b/src/Structure/H5PartWrapperForPC.h index e9eb295824a85297689af390a25669f71c192c1d..0b40133ddd12dfad93b8174ec97cde282b0eaf71 100644 --- a/src/Structure/H5PartWrapperForPC.h +++ b/src/Structure/H5PartWrapperForPC.h @@ -2,7 +2,7 @@ #define OPAL_H5PARTWRAPPERFORPC_H /*! - H5PartWrapperForPT: a class that manages calls to H5Part for the cyclotron tracker + H5PartWrapperForPC: a class that manages calls to H5Part for the cyclotron tracker */ #include "Structure/H5PartWrapper.h" diff --git a/src/Structure/H5PartWrapperForPS.cpp b/src/Structure/H5PartWrapperForPS.cpp index 7823bf3c5376069a9757fc0389c0ee91a52ed37b..0b2e8f7db324738dab4c321e94b63f283fa6da72 100644 --- a/src/Structure/H5PartWrapperForPS.cpp +++ b/src/Structure/H5PartWrapperForPS.cpp @@ -16,8 +16,6 @@ #include <sstream> #include <set> -extern Inform *gmsg; - H5PartWrapperForPS::H5PartWrapperForPS(const std::string &fileName, h5_int32_t flags): H5PartWrapper(fileName, flags) { } @@ -334,7 +332,7 @@ void H5PartWrapperForPS::writeStepHeader(PartBunchBase<double, 3>* bunch, } catch (std::out_of_range & m) { ERRORMSG(m.what() << endl); - throw OpalException("H5PartWrapperForPC::wirteStepHeader", + throw OpalException("H5PartWrapperForPC::writeStepHeader", "some additional step attribute not found"); } ++ numSteps_m; diff --git a/src/Structure/H5PartWrapperForPT.cpp b/src/Structure/H5PartWrapperForPT.cpp index 47e6cd468ed843e925feeeef8a904f9617d50111..8fe3529b4308a97f2fd90ae006f1aac7c8dc3dfe 100644 --- a/src/Structure/H5PartWrapperForPT.cpp +++ b/src/Structure/H5PartWrapperForPT.cpp @@ -16,8 +16,6 @@ #include <sstream> #include <set> -extern Inform *gmsg; - H5PartWrapperForPT::H5PartWrapperForPT(const std::string &fileName, h5_int32_t flags): H5PartWrapper(fileName, flags) { } @@ -194,34 +192,34 @@ void H5PartWrapperForPT::writeHeader() { std::stringstream OPAL_version; OPAL_version << OPAL_PROJECT_NAME << " " << OPAL_PROJECT_VERSION << " # git rev. " << Util::getGitRevision(); WRITESTRINGFILEATTRIB(file_m, "OPAL_version", OPAL_version.str().c_str()); - - + + WRITESTRINGFILEATTRIB(file_m, "idUnit", "1"); WRITESTRINGFILEATTRIB(file_m, "xUnit", "m"); WRITESTRINGFILEATTRIB(file_m, "yUnit", "m"); WRITESTRINGFILEATTRIB(file_m, "zUnit", "m"); - + WRITESTRINGFILEATTRIB(file_m, "pxUnit", "#beta#gamma"); WRITESTRINGFILEATTRIB(file_m, "pyUnit", "#beta#gamma"); WRITESTRINGFILEATTRIB(file_m, "pzUnit", "#beta#gamma"); - + WRITESTRINGFILEATTRIB(file_m, "ptypeUnit", "1"); WRITESTRINGFILEATTRIB(file_m, "qUnit", "C"); - + if (Options::ebDump) { WRITESTRINGFILEATTRIB(file_m, "ExUnit", "MV/m"); WRITESTRINGFILEATTRIB(file_m, "EyUnit", "MV/m"); WRITESTRINGFILEATTRIB(file_m, "EzUnit", "MV/m"); - + WRITESTRINGFILEATTRIB(file_m, "BxUnit", "T"); WRITESTRINGFILEATTRIB(file_m, "ByUnit", "T"); WRITESTRINGFILEATTRIB(file_m, "BzUnit", "T"); } - + if (Options::rhoDump) { WRITESTRINGFILEATTRIB(file_m, "rhoUnit", "C/m^3"); } - + WRITESTRINGFILEATTRIB(file_m, "TaitBryantAnglesUnit", "rad"); WRITESTRINGFILEATTRIB(file_m, "SPOSUnit", "m"); @@ -236,12 +234,12 @@ void H5PartWrapperForPT::writeHeader() { WRITESTRINGFILEATTRIB(file_m, "centroidUnit", "m"); WRITESTRINGFILEATTRIB(file_m, "RMSPUnit", "#beta#gamma"); WRITESTRINGFILEATTRIB(file_m, "MEANPUnit", "#beta#gamma"); - + WRITESTRINGFILEATTRIB(file_m, "minXUnit", "m"); WRITESTRINGFILEATTRIB(file_m, "maxXUnit", "m"); WRITESTRINGFILEATTRIB(file_m, "minPUnit", "#beta#gamma"); WRITESTRINGFILEATTRIB(file_m, "maxPUnit", "#beta#gamma"); - + WRITESTRINGFILEATTRIB(file_m, "dEUnit", "MeV"); WRITESTRINGFILEATTRIB(file_m, "MASSUnit", "GeV"); @@ -259,7 +257,7 @@ void H5PartWrapperForPT::writeHeader() { WRITESTRINGFILEATTRIB(file_m, "RefPartRUnit", "m"); WRITESTRINGFILEATTRIB(file_m, "RefPartPUnit", "#beta#gamma"); WRITESTRINGFILEATTRIB(file_m, "SteptoLastInjUnit", "1"); - + WRITESTRINGFILEATTRIB(file_m, "dump frequencyUnit", "1") WRITESTRINGFILEATTRIB(file_m, "dPhiGlobalUnit", "rad") @@ -378,7 +376,7 @@ void H5PartWrapperForPT::writeStepHeader(PartBunchBase<double, 3>* bunch, const } catch (std::out_of_range & m) { ERRORMSG(m.what() << endl); - throw OpalException("H5PartWrapperForPC::wirteStepHeader", + throw OpalException("H5PartWrapperForPC::writeStepHeader", "some additional step attribute not found"); } diff --git a/src/Structure/MonitorStatisticsWriter.cpp b/src/Structure/MonitorStatisticsWriter.cpp index ee67a48a5516715b03931599b58b2be66e19fdea..3aaf4c5514b95a41e1c6941cba43ab074f526ab1 100644 --- a/src/Structure/MonitorStatisticsWriter.cpp +++ b/src/Structure/MonitorStatisticsWriter.cpp @@ -1,6 +1,7 @@ #include "MonitorStatisticsWriter.h" -#include "Structure/LossDataSink.h" +#include "Structure/LossDataSink.h" +#include "Ippl.h" MonitorStatisticsWriter::MonitorStatisticsWriter(const std::string& fname, bool restart) : SDDSWriter(fname, restart) diff --git a/src/ValueDefinitions/RealVariable.cpp b/src/ValueDefinitions/RealVariable.cpp index e2dda88adfac675ed51c707419cc66e977771af8..402c50b33500884b4ca55f18624314646587abe7 100644 --- a/src/ValueDefinitions/RealVariable.cpp +++ b/src/ValueDefinitions/RealVariable.cpp @@ -19,7 +19,6 @@ #include "ValueDefinitions/RealVariable.h" #include "AbstractObjects/OpalData.h" #include "Attributes/Attributes.h" -#include "Utilities/Options.h" #include <iostream> diff --git a/src/ValueDefinitions/StringConstant.cpp b/src/ValueDefinitions/StringConstant.cpp index 200fc2ff757aadb13b2352424bbd38d6d75ff2a7..0b14913967943db0bacaa83dd43a6acbf79a11bf 100644 --- a/src/ValueDefinitions/StringConstant.cpp +++ b/src/ValueDefinitions/StringConstant.cpp @@ -19,7 +19,6 @@ #include "ValueDefinitions/StringConstant.h" #include "AbstractObjects/OpalData.h" #include "Attributes/Attributes.h" -#include "Utilities/Options.h" #include "Utilities/Util.h" #include <iostream>