Commit 590f45fc authored by snuverink_j's avatar snuverink_j
Browse files

various cosmetic improvements

parent 29f9b866
......@@ -196,95 +196,95 @@ Canada, July 2006.
// returns the hypervolume of ps[0 ..] in 3D
// assumes that ps is sorted improving
{
avl_init_node(ps.points[ps.nPoints-1].tnode,ps.points[ps.nPoints-1].objectives);
avl_insert_top(tree,ps.points[ps.nPoints-1].tnode);
avl_init_node(ps.points[ps.nPoints-1].tnode,ps.points[ps.nPoints-1].objectives);
avl_insert_top(tree,ps.points[ps.nPoints-1].tnode);
double hypera = (ref.objectives[0] - ps.points[ps.nPoints-1].objectives[0]) *
(ref.objectives[1] - ps.points[ps.nPoints-1].objectives[1]);
double hypera = (ref.objectives[0] - ps.points[ps.nPoints-1].objectives[0]) *
(ref.objectives[1] - ps.points[ps.nPoints-1].objectives[1]);
double height;
if (ps.nPoints == 1)
height = ref.objectives[2] - ps.points[ps.nPoints-1].objectives[2];
else
height = ps.points[ps.nPoints-2].objectives[2] - ps.points[ps.nPoints-1].objectives[2];
double hyperv = hypera * height;
for (int i = ps.nPoints - 2; i >= 0; i--)
{
if (i == 0)
height = ref.objectives[2] - ps.points[i].objectives[2];
double height;
if (ps.nPoints == 1)
height = ref.objectives[2] - ps.points[ps.nPoints-1].objectives[2];
else
height = ps.points[i-1].objectives[2] - ps.points[i].objectives[2];
// search tree for point q to the right of current point
const double * prv_ip, * nxt_ip;
avl_node_t *tnode;
avl_init_node(ps.points[i].tnode, ps.points[i].objectives);
if (avl_search_closest(tree, ps.points[i].objectives, &tnode) <= 0) {
nxt_ip = (double *)(tnode->item);
tnode = tnode->prev;
} else {
nxt_ip = (tnode->next!=NULL)
? (double *)(tnode->next->item)
: ref.objectives;
}
// if p is not dominated
if (nxt_ip[0] > ps.points[i].objectives[0]) {
// insert p in tree
avl_insert_after(tree, tnode, ps.points[i].tnode);
if (tnode !=NULL) {
prv_ip = (double *)(tnode->item);
if (prv_ip[0] > ps.points[i].objectives[0]) {
const double * cur_ip;
tnode = ps.points[i].tnode->prev;
// cur_ip = point dominated by pp with highest [0]-coordinate
cur_ip = (double *)(tnode->item);
// for each point in s in tree dominated by p
while (tnode->prev) {
prv_ip = (double *)(tnode->prev->item);
// decrease area by contribution of s
hypera -= (prv_ip[1] - cur_ip[1])*(nxt_ip[0] - cur_ip[0]);
if (prv_ip[0] < ps.points[i].objectives[0])
break; // prv is not dominated by pp
cur_ip = prv_ip;
// remove s from tree
avl_unlink_node(tree,tnode);
tnode = tnode->prev;
}
// remove s from tree
avl_unlink_node(tree,tnode);
if (!tnode->prev) {
// decrease area by contribution of s
hypera -= (ref.objectives[1] - cur_ip[1])*(nxt_ip[0] - cur_ip[0]);
prv_ip = ref.objectives;
}
}
} else
height = ps.points[ps.nPoints-2].objectives[2] - ps.points[ps.nPoints-1].objectives[2];
double hyperv = hypera * height;
for (int i = ps.nPoints - 2; i >= 0; i--) {
if (i == 0)
height = ref.objectives[2] - ps.points[i].objectives[2];
else
height = ps.points[i-1].objectives[2] - ps.points[i].objectives[2];
// search tree for point q to the right of current point
const double * prv_ip, * nxt_ip;
avl_node_t *tnode;
avl_init_node(ps.points[i].tnode, ps.points[i].objectives);
if (avl_search_closest(tree, ps.points[i].objectives, &tnode) <= 0) {
nxt_ip = (double *)(tnode->item);
tnode = tnode->prev;
} else {
nxt_ip = (tnode->next!=NULL)
? (double *)(tnode->next->item)
: ref.objectives;
}
// if p is not dominated
if (nxt_ip[0] > ps.points[i].objectives[0]) {
// insert p in tree
avl_insert_after(tree, tnode, ps.points[i].tnode);
if (tnode !=NULL) {
prv_ip = (double *)(tnode->item);
if (prv_ip[0] > ps.points[i].objectives[0]) {
const double * cur_ip;
tnode = ps.points[i].tnode->prev;
// cur_ip = point dominated by pp with highest [0]-coordinate
cur_ip = (double *)(tnode->item);
// for each point in s in tree dominated by p
while (tnode->prev) {
prv_ip = (double *)(tnode->prev->item);
// decrease area by contribution of s
hypera -= (prv_ip[1] - cur_ip[1])*(nxt_ip[0] - cur_ip[0]);
if (prv_ip[0] < ps.points[i].objectives[0])
break; // prv is not dominated by pp
cur_ip = prv_ip;
// remove s from tree
avl_unlink_node(tree,tnode);
tnode = tnode->prev;
}
// remove s from tree
avl_unlink_node(tree,tnode);
if (!tnode->prev) {
// decrease area by contribution of s
hypera -= (ref.objectives[1] - cur_ip[1])*(nxt_ip[0] - cur_ip[0]);
prv_ip = ref.objectives;
// increase area by contribution of p
hypera += (prv_ip[1] -
ps.points[i].objectives[1])*(nxt_ip[0] -
ps.points[i].objectives[0]);
}
}
} else
prv_ip = ref.objectives;
// increase area by contribution of p
hypera += (prv_ip[1] -
ps.points[i].objectives[1])*(nxt_ip[0] -
ps.points[i].objectives[0]);
if (height > 0)
hyperv += hypera * height;
}
avl_clear_tree(tree);
return hyperv;
}
if (height > 0)
hyperv += hypera * height;
}
avl_clear_tree(tree);
return hyperv;
}
double inclhv(POINT p)
......
......@@ -167,9 +167,9 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
if(j == 0) {
phi = initialPhase;
E = Emax;
j = -1;
// j = -1;
do {
j ++;
// j ++;
Emax = E;
initialPhase = phi;
phi += dphi;
......
......@@ -104,7 +104,7 @@ Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);
#define PSdim 6
//#define PSdim 6
extern Inform *gmsg;
......@@ -507,7 +507,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << endl;
double sym = elptr->getSymmetry();
*gmsg << "* " << sym << "-fold field symmerty " << endl;
*gmsg << "* " << sym << "-fold field symmetry " << endl;
// ckr: this just returned the default value as defined in Component.h
// double rff = elptr->getRfFrequ();
......@@ -2725,7 +2725,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
//itsBunch_m->R *= Vector_t(0.001); // mm --> m
// -------------- If flag DumpLocalFrame is not set, dump bunch in global frame ------------- //
// -------------- If flag psDumpFrame is global, dump bunch in global frame ------------- //
if (Options::psDumpFreq < std::numeric_limits<int>::max() ){
if (Options::psDumpFrame == Options::GLOBAL) {
......@@ -2757,7 +2757,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
}
}
// ---------------- If flag DumpLocalFrame is set, dump bunch in local frame ---------------- //
// ---------------- If flag psDumpFrame is local, dump bunch in local frame ---------------- //
else {
// The bunch is transformed into a local coordinate system with meanP in direction of y-axis
......
#include "Algorithms/ParallelSliceTracker.h"
#include <cfloat>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <memory>
#include <string>
#include "AbstractObjects/OpalData.h"
#include "Beamlines/Beamline.h"
......@@ -259,21 +257,21 @@ void ParallelSliceTracker::computeSpaceChargeFields() {
void ParallelSliceTracker::dumpStats(long long step) {
double sposRef = itsBunch_m->get_sPos();
if (step != 0 && (step % Options::psDumpFreq == 0 || step % Options::statDumpFreq == 0))
writePhaseSpace(step, sposRef);
double sposRef = itsBunch_m->get_sPos();
if (step != 0 && (step % Options::psDumpFreq == 0 || step % Options::statDumpFreq == 0))
writePhaseSpace(step, sposRef);
}
void ParallelSliceTracker::switchElements(double scaleMargin) {
double margin = 1.0;
double margin = 1.0;
itsOpalBeamline_m->resetStatus();
currentSimulationTime_m = itsBunch_m->getT();
itsOpalBeamline_m->switchElements(itsBunch_m->zTail() - margin,
itsBunch_m->zHead() + margin,
itsBunch_m->Eavg() * 1e-6);
itsOpalBeamline_m->resetStatus();
currentSimulationTime_m = itsBunch_m->getT();
itsOpalBeamline_m->switchElements(itsBunch_m->zTail() - margin,
itsBunch_m->zHead() + margin,
itsBunch_m->Eavg() * 1e-6);
}
......
#ifndef OPAL_ParallelSliceTracker_HH
#define OPAL_ParallelSliceTracker_HH
#include <vector>
#include <list>
#include <memory>
#include "Algorithms/bet/EnvelopeBunch.h"
......
......@@ -740,8 +740,6 @@ void ParallelTTracker::computeWakefield(IndexMap::value_t &elements) {
void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth) {
Inform msg("ParallelTTracker ", *gmsg);
IndexMap::value_t::const_iterator it = elements.begin();
const IndexMap::value_t::const_iterator end = elements.end();
std::set<IndexMap::value_t::value_type> elementsWithParticleMatterInteraction;
std::set<ParticleMatterInteractionHandler*> particleMatterinteractionHandlers;
std::pair<double, double> currentRange(0.0, 0.0);
......
......@@ -872,7 +872,7 @@ std::vector<double> StatisticalErrors::interpolateSDDSData(const std::vector<dou
std::vector<double> interpolated(newSize);
for (unsigned int i = 0; i < newSize; ++ i) {
while (oldPositions[j] <= newPositions[i] && j < oldSize) {
while (j < oldSize && oldPositions[j] <= newPositions[i]) {
++ j;
}
if (j == 0 || j == oldSize) {
......
......@@ -40,7 +40,10 @@ extern Inform *gmsg;
beta = 0.85 -> 460 keV
*/
#define BETA_MIN1 0.30 // minimum beta-value for space-charge calculations: start
#ifndef USE_HOMDYN_SC_MODEL
#define BETA_MIN2 0.45 // minimum beta-value for space-charge calculations: full impact
#endif
// Hack allows odeint in rk.C to be called with a class member function
static EnvelopeBunch *thisBunch = NULL; // pointer to access calling bunch
......
......@@ -740,14 +740,14 @@ void Bend::findBendEffectiveLength(double startField, double endField) {
bendAngle1 = newBendAngle;
delta1 = delta;
} else {
bendAngle2 = newBendAngle;
// bendAngle2 = newBendAngle;
delta2 = delta;
}
} else {
if(newBendAngle - angle_m < 0.0) {
bendAngle2 = newBendAngle;
// bendAngle2 = newBendAngle;
delta2 = delta;
} else {
bendAngle1 = newBendAngle;
......@@ -822,14 +822,14 @@ void Bend::findBendStrength(double mass,
bendAngle1 = newBendAngle;
amplitude1 = fieldAmplitude_m;
} else {
bendAngle2 = newBendAngle;
// bendAngle2 = newBendAngle;
amplitude2 = fieldAmplitude_m;
}
} else {
if(newBendAngle - angle_m < 0.0) {
bendAngle2 = newBendAngle;
// bendAngle2 = newBendAngle;
amplitude2 = fieldAmplitude_m;
} else {
bendAngle1 = newBendAngle;
......
......@@ -32,8 +32,6 @@
extern Inform *gmsg;
using namespace std;
// Class Collimator
// ------------------------------------------------------------------------
......@@ -371,8 +369,8 @@ void Collimator::print() {
std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
ERRORMSG(errormsg << endl);
if (Ippl::myNode() == 0) {
ofstream omsg("errormsg.txt", ios_base::app);
omsg << errormsg << endl;
std::ofstream omsg("errormsg.txt", std::ios_base::app);
omsg << errormsg << std::endl;
omsg.close();
}
informed_m = true;
......@@ -416,7 +414,7 @@ void Collimator::setOutputFN(std::string fn) {
filename_m = fn;
}
string Collimator::getOutputFN() {
std::string Collimator::getOutputFN() {
if (filename_m == std::string(""))
return getName();
else
......@@ -432,7 +430,7 @@ ElementBase::ElementType Collimator::getType() const {
return COLLIMATOR;
}
string Collimator::getCollimatorShape() {
std::string Collimator::getCollimatorShape() {
if (isAPepperPot_m)
return "PeperPot";
else if (isASlit_m)
......
......@@ -27,8 +27,6 @@
#include "Fields/EMField.h"
#include "Algorithms/Vektor.h"
#define EPS_MISALIGNMENT 1e-8
class PartData;
template <class T, unsigned Dim>
......
......@@ -89,7 +89,10 @@ Cyclotron::~Cyclotron() {
void Cyclotron::applyTrimCoil(const double r, const double z, double *br, double *bz) {
/// update bz and br with trim coil contributions
// for some discussion on the formulas see https://gitlab.psi.ch/OPAL/src/issues/110
// for some discussion on the formulas see
// http://doi.org/10.1103/PhysRevSTAB.14.054402
// https://gitlab.psi.ch/OPAL/src/issues/157
// https://gitlab.psi.ch/OPAL/src/issues/110
// unitless constants
const double Amax1 = 1;
......@@ -750,7 +753,7 @@ double Cyclotron::gutdf5d(double *f, double dx, const int kor, const int krl, co
}
// evaulate other derivative of magnetic field.
// evaluate other derivative of magnetic field.
void Cyclotron::getdiffs() {
Bfield.dbr.resize(Bfield.ntot);
......
......@@ -160,7 +160,7 @@ bool VariableRFCavity::apply(const Vector_t &R, const Vector_t &P,
}
double E0 = _amplitude_td->getValue(t);
double f = _frequency_td->getValue(t) * 1.0E-3; // neet GHz on the element we have MHz
double f = _frequency_td->getValue(t) * 1.0E-3; // need GHz on the element we have MHz
double phi = _phase_td->getValue(t);
E = Vector_t(0., 0., E0*sin(Physics::two_pi * f * t + phi));
return false;
......
......@@ -4361,7 +4361,7 @@ void Distribution::setupParticleBins(double massIneV, PartBunchBase<double, 3> *
if (!itsAttr[Attrib::Legacy::Distribution::PT].defaultUsed())
throw OpalException("Distribution::setupParticleBins",
"PT is obsolet. The moments of the beam is defined with OFFSETPZ");
"PT is obsolete. The moments of the beam is defined with OFFSETPZ");
// we get gamma from PC of the beam
const double pz = beam->getP()/beam->getM();
......@@ -4460,7 +4460,7 @@ void Distribution::shiftDistCoordinates(double massIneV) {
double deltaPz = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETPZ]);
if (Attributes::getReal(currDist->itsAttr[Attrib::Legacy::Distribution::PT])!=0.0)
WARNMSG("PT & PZ are obsolet and will be ignored. The moments of the beam is defined with PC" << endl);
WARNMSG("PT & PZ are obsolete and will be ignored. The moments of the beam is defined with PC" << endl);
// Check input momentum units.
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
......@@ -4485,61 +4485,62 @@ void Distribution::shiftDistCoordinates(double massIneV) {
void Distribution::writeOutFileHeader() {
if (Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE])) {
if (Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE]) == false)
return;
unsigned int totalNum = tOrZDist_m.size();
reduce(totalNum, totalNum, OpAddAssign());
if (Ippl::myNode() == 0) {
std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
*gmsg << "* **********************************************************************************" << endl;
*gmsg << "* Write initial distribution to file \"" << fname << "\"" << endl;
*gmsg << "* **********************************************************************************" << endl;
std::ofstream outputFile(fname);
if (outputFile.bad()) {
*gmsg << "Unable to open output file \"" << fname << "\"" << endl;
} else {
outputFile.setf(std::ios::left);
outputFile << "# ";
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 (numberOfEnergyBins_m > 0) {
outputFile.width(17);
outputFile << "Bin Number";
}
outputFile << std::endl;
unsigned int totalNum = tOrZDist_m.size();
reduce(totalNum, totalNum, OpAddAssign());
if (Ippl::myNode() != 0)
return;
outputFile << "# " << totalNum << std::endl;
}
std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
*gmsg << "* **********************************************************************************" << endl;
*gmsg << "* Write initial distribution to file \"" << fname << "\"" << endl;
*gmsg << "* **********************************************************************************" << endl;
std::ofstream outputFile(fname);
if (outputFile.bad()) {
*gmsg << "Unable to open output file \"" << fname << "\"" << endl;
} else {
outputFile.setf(std::ios::left);
outputFile << "# ";
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 (numberOfEnergyBins_m > 0) {
outputFile.width(17);
outputFile << "Bin Number";
}
outputFile.close();
outputFile << std::endl;
outputFile << "# " << totalNum << std::endl;
}
}
outputFile.close();
}
void Distribution::writeOutFileEmission() {
......@@ -4669,55 +4670,55 @@ void Distribution::writeOutFileEmission() {
void Distribution::writeOutFileInjection() {
if (Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE])) {
if (Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE]) == false)
return;
std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
// Nodes take turn writing particles to file.
for (int nodeIndex = 0; nodeIndex < Ippl::getNodes(); nodeIndex++) {
// Write to file if its our turn.
size_t numberOfParticles = 0;
if (Ippl::myNode() == nodeIndex) {
std::ofstream outputFile(fname, std::fstream::app);
if (outputFile.bad()) {
*gmsg << "Node " << Ippl::myNode() << " unable to write"
<< "to file \"" << fname << "\"" << endl;
} else {
std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
// Nodes take turn writing particles to file.
for (int nodeIndex = 0; nodeIndex < Ippl::getNodes(); nodeIndex++) {
outputFile.precision(9);
outputFile.setf(std::ios::scientific);
outputFile.setf(std::ios::right);
// Write to file if its our turn.
size_t numberOfParticles = 0;
if (Ippl::myNode() == nodeIndex) {
std::ofstream outputFile(fname, std::fstream::app);
if (outputFile.bad()) {