Commit a7f37041 authored by snuverink_j's avatar snuverink_j
Browse files

cleaning, refactoring, spelling, spacing

parent 7dff022f
......@@ -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);
......
......@@ -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
......
......@@ -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) {
......
......@@ -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();
......
......@@ -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
......@@ -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
......
......@@ -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
......
......@@ -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 ("
......
......@@ -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 {
......
......@@ -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"
......
......@@ -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;
......
......@@ -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;
......
#ifndef OPAL_VEKTOR_HH
#define OPAL_VEKTOR_HH
#include <limits>
#include "AppTypes/Vektor.h"
typedef Vektor<double, 3> Vector_t;
......
......@@ -73,7 +73,7 @@ void Interpolator3dGridTo3d::setAll(ThreeDGrid* grid,
break;
default:
throw(LogicalError(
"Interpolator3dGridTo3d::Set",
"Interpolator3dGridTo3d::setAll",
"Did not recognise interpolation algorithm"
)
);
......
......@@ -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];
......
......@@ -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"
));
}
......
......@@ -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:
......
......@@ -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) {