Skip to content

GitLab

Commit a6e65681 authored by kraus's avatar kraus
Browse files

cleaning up Bend, ParallelTTracker and OpalBeamline

parent cf4044fc
Hide whitespace changes
Inline Side-by-side
Showing with 13 additions and 239 deletions
src/Algorithms/ParallelTTracker.cpp
View file @ a6e65681
......@@ -84,7 +84,6 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
particleMaterStatus_m(false),
totalParticlesInSimulation_m(0)
// , logger_m("designPath_" + std::to_string(Ippl::myNode()) + ".dat")
{
CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
......@@ -130,7 +129,6 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
particleMaterStatus_m(false),
totalParticlesInSimulation_m(0)
// , logger_m("designPath_" + std::to_string(Ippl::myNode()) + ".dat")
{
CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
......@@ -359,7 +357,7 @@ void ParallelTTracker::execute() {
itsBunch_m->setT(t);
if (t > 0.0) {
updateRefToLabCSTrafo(pusher);
updateReference(pusher);
}
if (deletedParticles_m) {
......@@ -1123,12 +1121,15 @@ void ParallelTTracker::writePhaseSpace(const long long step, bool psDump, bool s
}
}
void ParallelTTracker::updateReference(const BorisPusher &pusher) {
updateReferenceParticle(pusher);
updateRefToLabCSTrafo();
}
void ParallelTTracker::updateReferenceParticle(const BorisPusher &pusher) {
//static size_t step = 0;
const double dt = std::min(itsBunch_m->getT(), itsBunch_m->getdT());
const double scaleFactor = Physics::c * dt;
Vector_t Ef(0.0), Bf(0.0);
// Vector_t oldR = RefPartR_m;
RefPartR_m /= scaleFactor;
pusher.push(RefPartR_m, RefPartP_m, dt);
......@@ -1158,35 +1159,11 @@ void ParallelTTracker::updateReferenceParticle(const BorisPusher &pusher) {
Bf += refToLocalCSTrafo.rotateFrom(localB);
}
// if (step % loggingFrequency_m == 0) {
// Vector_t R = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
// Vector_t P = referenceToLabCSTrafo_m.rotateTo(RefPartP_m);
// Vector_t lEf = referenceToLabCSTrafo_m.rotateTo(Ef);
// Vector_t lBf = referenceToLabCSTrafo_m.rotateTo(Bf);
// logger_m << std::setw(18) << std::setprecision(8) << pathLength_m + euclidean_norm(RefPartR_m - oldR)
// << std::setw(18) << std::setprecision(8) << R(0)
// << std::setw(18) << std::setprecision(8) << R(1)
// << std::setw(18) << std::setprecision(8) << R(2)
// << std::setw(18) << std::setprecision(8) << P(0)
// << std::setw(18) << std::setprecision(8) << P(1)
// << std::setw(18) << std::setprecision(8) << P(2)
// << std::setw(18) << std::setprecision(8) << lEf(0)
// << std::setw(18) << std::setprecision(8) << lEf(1)
// << std::setw(18) << std::setprecision(8) << lEf(2)
// << std::setw(18) << std::setprecision(8) << lBf(0)
// << std::setw(18) << std::setprecision(8) << lBf(1)
// << std::setw(18) << std::setprecision(8) << lBf(2)
// << std::setw(18) << std::setprecision(8) << Util::getEnergy(RefPartP_m, itsBunch_m->getM() * 1e-6)
// << std::setw(18) << std::setprecision(8) << (itsBunch_m->getT() - 0.5 * itsBunch_m->getdT()) * 1e9
// << std::endl;
// }
pusher.kick(RefPartR_m, RefPartP_m, Ef, Bf, dt);
RefPartR_m /= scaleFactor;
pusher.push(RefPartR_m, RefPartP_m, dt);
RefPartR_m *= scaleFactor;
//++ step;
}
void ParallelTTracker::transformBunch(const CoordinateSystemTrafo &trafo) {
......@@ -1199,9 +1176,7 @@ void ParallelTTracker::transformBunch(const CoordinateSystemTrafo &trafo) {
}
}
void ParallelTTracker::updateRefToLabCSTrafo(const BorisPusher &pusher) {
updateReferenceParticle(pusher);
void ParallelTTracker::updateRefToLabCSTrafo() {
pathLength_m += euclidean_norm(RefPartR_m);
CoordinateSystemTrafo update(RefPartR_m,
......
src/Algorithms/ParallelTTracker.h
View file @ a6e65681
......@@ -315,16 +315,14 @@ private:
void transformBunch(const CoordinateSystemTrafo &trafo);
void updateRefToLabCSTrafo(const BorisPusher &pusher);
void updateReference(const BorisPusher &pusher);
void updateRefToLabCSTrafo();
void findStartPosition(const BorisPusher &pusher);
void autophaseCavities(const BorisPusher &pusher);
void evenlyDistributeParticles();
static unsigned long long getMaxSteps(std::queue<unsigned long long> numSteps);
// std::ofstream logger_m;
// size_t loggingFrequency_m;
};
inline void ParallelTTracker::visitAlignWrapper(const AlignWrapper &wrap) {
......@@ -570,4 +568,4 @@ void ParallelTTracker::kickParticlesDKS() {
}
#endif
#endif // OPAL_ParallelTTracker_HH
#endif // OPAL_ParallelTTracker_HH
\ No newline at end of file
src/Classic/AbsBeamline/Bend.cpp
View file @ a6e65681
......@@ -1436,76 +1436,6 @@ bool Bend::treatAsDrift(Inform &msg, double chordLength) {
return false;
}
void Bend::retrieveDesignEnergy(double startField) {
// energyEvolution_t::iterator it = OpalData::getInstance()->getFirstEnergyData();
// energyEvolution_t::iterator end = OpalData::getInstance()->getLastEnergyData();
// for (; it != end; ++ it) {
// if ((*it).first > startField) break;
// designEnergy_m = (*it).second;
// }
}
std::pair<Vector_t, Vector_t> Bend::getDesignPathSecant(double startsAtDistFromEdge, double length) const {
length = std::abs(length);
Vector_t startPosition(0.0);
Vector_t tangent(0, 0, 1);
double pathLength = refTrajMap_m[0](2);
unsigned int size = refTrajMap_m.size();
for (unsigned int i = 1; i < size; ++ i) {
Vector_t step = refTrajMap_m[i] - refTrajMap_m[i-1];
double stepSize = euclidean_norm(step);
if (pathLength + stepSize > startsAtDistFromEdge) {
double diff = startsAtDistFromEdge - pathLength;
tangent = step / stepSize;
startPosition = refTrajMap_m[i-1] + diff * tangent;
if (length <= 1e-12) {
return std::make_pair(startPosition, tangent);
}
for (unsigned j = i; j < size; ++ j) {
Vector_t position = refTrajMap_m[j];
if (euclidean_norm(position - startPosition) >= length) {
step = refTrajMap_m[j-1] - refTrajMap_m[j];
double tau = (dot(startPosition - position, step) - sqrt(std::pow(dot(startPosition - position, step), 2) - dot(step, step) * (dot(startPosition - position, startPosition - position) - std::pow(length, 2)))) / dot(step, step);
tangent = position + tau * step - startPosition;
tangent /= euclidean_norm(tangent);
return std::make_pair(startPosition, tangent);
}
}
Vector_t position = refTrajMap_m[size - 1];
Vector_t step = position - refTrajMap_m[size - 2];
step /= euclidean_norm(step);
double tau = (dot(startPosition - position, step) + sqrt(std::pow(dot(startPosition - position, step), 2) - dot(step, step) * (dot(startPosition - position, startPosition - position) - std::pow(length, 2)))) / dot(step, step);
tangent = position + tau * step - startPosition;
tangent /= euclidean_norm(tangent);
return std::make_pair(startPosition, tangent);
}
pathLength += stepSize;
}
double diff = startsAtDistFromEdge - pathLength;
unsigned int i = size - 1;
Vector_t step = refTrajMap_m[i] - refTrajMap_m[i-1];
double stepSize = euclidean_norm(step);
tangent = step / stepSize;
startPosition = refTrajMap_m[i] + diff * tangent;
//tangent = tangent;
return std::make_pair(startPosition, tangent);
}
std::vector<Vector_t> Bend::getOutline() const {
std::vector<Vector_t> outline;
Vector_t rotationCenter = Vector_t(-designRadius_m * cosEntranceAngle_m, 0.0, designRadius_m * sinEntranceAngle_m);
......@@ -1907,10 +1837,4 @@ std::array<double,2> Bend::getExitFringeFieldLength() const {
return extFFL;
}
}
\ No newline at end of file
src/Classic/AbsBeamline/Bend.h
View file @ a6e65681
......@@ -124,8 +124,6 @@ public:
double getExitAngle() const;
double getMapLength() const;
std::pair<Vector_t, Vector_t> getDesignPathSecant(double startsAtDistFromEdge, double length) const;
std::vector<Vector_t> getOutline() const;
MeshData getSurfaceMesh() const;
......@@ -213,7 +211,6 @@ private:
void setFieldBoundaries(double startField, double endField);
void setupPusher(PartBunchBase<double, 3> *bunch);
bool treatAsDrift(Inform &msg, double chordlength);
void retrieveDesignEnergy(double startField);
void setCSTrafoToEntranceRegion(const CoordinateSystemTrafo &trafo);
void setCSTrafoToExitRegion(const CoordinateSystemTrafo &trafo);
......
src/Elements/OpalBeamline.cpp
View file @ a6e65681
......@@ -354,33 +354,11 @@ void OpalBeamline::compute3DLattice() {
endPriorPathLength = beginThisPathLength + arcLength;
} else {
// Quaternion rotation(1, 0, 0, 0);
// FieldList::iterator priorDipole = partiallyInsideDipole(it, elements_m.begin(), elements_m.end(), minOrder);
// if (priorDipole != it) {
// Bend * bendElement = static_cast<Bend*>((*priorDipole).getElement().get());
// double pathDifference = beginThisPathLength - bendElement->getElementPosition();
// auto secant = bendElement->getDesignPathSecant(pathDifference, thisLength);
// Vector_t position = (*priorDipole).getCoordTransformationTo().transformFrom(secant.first);
// Vector_t orientation = (*priorDipole).getCoordTransformationTo().rotateFrom(secant.second);
// beginThis3D = currentCoordTrafo.transformTo(position);
// rotation = getQuaternion(orientation,
// currentCoordTrafo.rotateFrom(Vector_t(0, 0, 1)));
// CoordinateSystemTrafo fromLastToThis(beginThis3D, rotation);
// fromLastToThis *= currentCoordTrafo;
// Vector_t origin = fromLastToThis.getOrigin();
// Vector_t end = origin + thisLength * fromLastToThis.rotateFrom(Vector_t(0,0,1));
// }
double rotationAngleAboutZ = (*it).getElement()->getRotationAboutZ();
Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ),
sin(-0.5 * rotationAngleAboutZ) * Vector_t(0, 0, 1));
CoordinateSystemTrafo fromLastToThis(beginThis3D, rotationAboutZ);// * rotation);
CoordinateSystemTrafo fromLastToThis(beginThis3D, rotationAboutZ);
(*it).setCoordTransformationTo(fromLastToThis * currentCoordTrafo);
}
......@@ -712,55 +690,6 @@ void OpalBeamline::save3DInput() {
pos << input << std::endl;
}
FieldList::iterator OpalBeamline::partiallyInsideDipole(const FieldList::iterator &it,
const FieldList::iterator &begin,
const FieldList::iterator &end,
const unsigned int &minOrder) {
if (it == begin) return it;
FieldList::iterator prior = it;
-- prior;
while (true) {
std::shared_ptr<Component> element = (*prior).getElement();
if ((*prior).getEnd() > /*(*it).getStart()*/ (*it).getElement()->getElementPosition() &&
(element->getType() == ElementBase::SBEND ||
element->getType() == ElementBase::RBEND)) {
if ((*prior).order_m >= minOrder)
return prior;
}
if (prior == begin) break;
-- prior;
}
if (it == end) return it;
FieldList::iterator next = it;
++ next;
if (next == end) return it;
while (true) {
std::shared_ptr<Component> element = (*next).getElement();
if ((element->getType() == ElementBase::SBEND ||
element->getType() == ElementBase::RBEND) &&
(*it).getElement()->getElementPosition() > (*next).getStart() &&
(*it).getElement()->getElementPosition() < (*next).getEnd()) {
if ((*next).order_m >= minOrder)
return next;
}
++ next;
if (next == end) break;
}
return it;
}
void OpalBeamline::activateElements() {
auto it = elements_m.begin();
const auto end = elements_m.end();
......
src/Elements/OpalBeamline.h
View file @ a6e65681
......@@ -65,25 +65,10 @@ public:
CoordinateSystemTrafo getCSTrafoLab2Local(const std::shared_ptr<Component> &comp) const;
CoordinateSystemTrafo getCSTrafoLab2Local() const;
CoordinateSystemTrafo getMisalignment(const std::shared_ptr<Component> &comp) const;
// void getSectionIndexAt(const Vector_t &, long &) const;
// double getSectionStart(const long &) const;
// double getSectionEnd(const unsigned int &) const;
// double getSectionEnd(const Vector_t &, long);
double getStart(const Vector_t &) const;
double getEnd(const Vector_t &) const;
// void setOrientation(const Vector_t &, const Vector_t &);
// void setOrientation(const Vector_t &, const unsigned int &);
// void updateOrientation(const Vector_t &, const Vector_t &, const double &, const double &);
// const Vector_t &getOrientation(const Vector_t &) const;
// const Vector_t &getOrientation(const long &) const;
// void resetStatus();
// void setStatus(const unsigned int &, const bool &);
// const bool &getStatus(const unsigned int &) const;
void switchElements(const double &, const double &, const double &kineticEnergy, const bool &nomonitors = false);
void switchAllElements();
......@@ -123,11 +108,6 @@ public:
bool containsSource();
private:
FieldList::iterator partiallyInsideDipole(const FieldList::iterator &it,
const FieldList::iterator &begin,
const FieldList::iterator &end,
const unsigned int &minOrder);
FieldList elements_m;
bool prepared_m;
bool containsSource_m;
......@@ -138,35 +118,6 @@ private:
static OpalSection dummy_section_m;
};
// inline WakeFunction *OpalBeamline::getWakeFunction(const unsigned int &index) {
// if(index < sections_m.size()) {
// return sections_m[index].getWakeFunction();
// }
// return NULL;
// }
// inline std::shared_ptr<const ElementBase> OpalBeamline::getWakeFunctionOwner(const unsigned int &index) {
// if(index < sections_m.size()) {
// return sections_m[index].getWakeFunctionOwner();
// }
// return NULL;
// }
// inline BoundaryGeometry *OpalBeamline::getBoundaryGeometry(const unsigned int &index) {
// if(index < sections_m.size()) {
// return sections_m[index].getBoundaryGeometry();
// }
// return NULL;
// }
// inline ParticleMatterInteractionHandler *OpalBeamline::getParticleMatterInteractionHandler(const unsigned int &index) {
// if(index < sections_m.size()) {
// return sections_m[index].getParticleMatterInteractionHandler();
// }
// return 0;
// }
template<class T> inline
void OpalBeamline::visit(const T &element, BeamlineVisitor &, PartBunchBase<double, 3> *bunch) {
Inform msg("OPAL ");
......@@ -298,4 +249,4 @@ inline
bool OpalBeamline::containsSource() {
return containsSource_m;
}
#endif // OPAL_BEAMLINE_H
#endif // OPAL_BEAMLINE_H
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment