Skip to content

GitLab

Commit 22bf204f authored by kraus's avatar kraus
Browse files

Bounding Box and computation of point of intersection seem to work

parent d2573049
Hide whitespace changes
Inline Side-by-side
Showing with 90 additions and 108 deletions
src/Algorithms/OrbitThreader.cpp
View file @ 22bf204f
......@@ -31,7 +31,6 @@ OrbitThreader::OrbitThreader(const PartData &ref,
double maxDiffZBunch,
double t,
double dt,
size_t maxIntegSteps,
double zstop,
OpalBeamline &bl) :
r_m(r),
......@@ -39,7 +38,6 @@ OrbitThreader::OrbitThreader(const PartData &ref,
pathLength_m(s),
time_m(t),
dt_m(dt),
maxIntegSteps_m(maxIntegSteps),
zstop_m(zstop),
itsOpalBeamline_m(bl),
errorFlag_m(0),
......@@ -82,17 +80,16 @@ OrbitThreader::OrbitThreader(const PartData &ref,
}
distTrackBack_m = std::min(pathLength_m, std::max(0.0, maxDiffZBunch));
computeBoundingBox();
}
void OrbitThreader::execute() {
double initialPathLength = pathLength_m;
computeMaximalImplicitDrift2();
double maxDistance = computeMaximalImplicitDrift();
auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
std::set<std::string> visitedElements;
trackBack(2 * maxDistance);
trackBack();
Vector_t nextR = r_m / (Physics::c * dt_m);
integrator_m.push(nextR, p_m, dt_m);
......@@ -110,7 +107,8 @@ void OrbitThreader::execute() {
double initialS = pathLength_m;
Vector_t initialR = r_m;
Vector_t initialP = p_m;
integrate(elementSet, maxIntegSteps_m, 2 * maxDistance);
double maxDistance = computeDriftLengthToBoundingBox(elementSet, r_m, p_m);
integrate(elementSet, maxDistance);
registerElement(elementSet, initialS, initialR, initialP);
......@@ -147,7 +145,7 @@ void OrbitThreader::execute() {
processElementRegister();
}
void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t /*maxSteps*/, double maxDrift) {
void OrbitThreader::integrate(const IndexMap::value_t &activeSet, double maxDrift) {
static size_t step = 0;
CoordinateSystemTrafo labToBeamline = itsOpalBeamline_m.getCSTrafoLab2Local();
const double oldPathLength = pathLength_m;
......@@ -282,7 +280,7 @@ double OrbitThreader::getMaxDesignEnergy(const IndexMap::value_t &elementSet) co
return designEnergy;
}
void OrbitThreader::trackBack(double maxDrift) {
void OrbitThreader::trackBack() {
dt_m *= -1;
double initialPathLength = pathLength_m;
......@@ -290,11 +288,12 @@ void OrbitThreader::trackBack(double maxDrift) {
integrator_m.push(nextR, p_m, dt_m);
nextR *= Physics::c * dt_m;
maxDrift = std::min(maxDrift, distTrackBack_m);
while (std::abs(initialPathLength - pathLength_m) < distTrackBack_m) {
auto elementSet = itsOpalBeamline_m.getElements(nextR);
integrate(elementSet, 1000, maxDrift);
double maxDrift = computeDriftLengthToBoundingBox(elementSet, r_m, -p_m);
maxDrift = std::min(maxDrift, distTrackBack_m);
integrate(elementSet, maxDrift);
nextR = r_m / (Physics::c * dt_m);
integrator_m.push(nextR, p_m, dt_m);
......@@ -389,7 +388,7 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
}
}
void OrbitThreader::computeMaximalImplicitDrift2() {
void OrbitThreader::computeBoundingBox() {
FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
FieldList::iterator it = allElements.begin();
const FieldList::iterator end = allElements.end();
......@@ -403,76 +402,40 @@ void OrbitThreader::computeMaximalImplicitDrift2() {
ElementBase::BoundingBox other = it->getBoundingBoxInLabCoords();
globalBoundingBox_m.getCombinedBoundingBox(other);
}
}
double OrbitThreader::computeMaximalImplicitDrift() {
FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
double maxDrift = 0.0;
MarkerRep start("#S");
CoordinateSystemTrafo toEdge(r_m, getQuaternion(p_m, Vector_t(0, 0, 1)));
start.setElementLength(0.0);
start.setCSTrafoGlobal2Local(toEdge);
std::shared_ptr<Component> startPtr(static_cast<Marker*>(start.clone()));
allElements.push_front(ClassicField(startPtr, 0.0, 0.0));
FieldList::iterator it = allElements.begin();
const FieldList::iterator end = allElements.end();
for (; it != end; ++ it) {
auto element1 = it->getElement();
const auto &toEdge = element1->getCSTrafoGlobal2Local();
auto toEnd = element1->getEdgeToEnd() * toEdge;
Vector_t end1 = toEnd.transformFrom(Vector_t(0, 0, 0));
Vector_t directionEnd = toEnd.rotateFrom(Vector_t(0, 0, 1));
if (element1->getType() == ElementBase::RBEND ||
element1->getType() == ElementBase::SBEND ||
element1->getType() == ElementBase::RBEND3D ) {
auto toBegin = element1->getEdgeToBegin() * toEdge;
BendBase *bend = static_cast<BendBase*>(element1.get());
double angleRelativeToFace = bend->getEntranceAngle() - bend->getBendAngle();
directionEnd = toBegin.rotateFrom(Vector_t(sin(angleRelativeToFace), 0, cos(angleRelativeToFace)));
}
double minDistanceLocal = std::numeric_limits<double>::max();
FieldList::iterator it2 = allElements.begin();
for (; it2 != end; ++ it2) {
if (it == it2) continue;
auto element2 = it2->getElement();
const auto &toEdge = element2->getCSTrafoGlobal2Local();
auto toBegin = element2->getEdgeToBegin() * toEdge;
auto toEnd = element2->getEdgeToEnd() * toEdge;
Vector_t begin2 = toBegin.transformFrom(Vector_t(0, 0, 0));
Vector_t end2 = toEnd.transformFrom(Vector_t(0, 0, 0));
Vector_t directionBegin = toBegin.rotateFrom(Vector_t(0, 0, 1));
if (element2->getType() == ElementBase::RBEND ||
element2->getType() == ElementBase::SBEND ||
element2->getType() == ElementBase::RBEND3D ) {
BendBase *bend = static_cast<BendBase*>(element2.get());
double E1 = bend->getEntranceAngle();
directionBegin = toBegin.rotateFrom(Vector_t(sin(E1), 0, cos(E1)));
}
double distance = euclidean_norm(begin2 - end1);
double directionProjection = dot(directionEnd, directionBegin);
bool overlapping = dot(begin2 - end1, directionBegin) < 0.0? true: false;
if (!overlapping &&
directionProjection > 0.999 &&
minDistanceLocal > distance) {
minDistanceLocal = distance;
}
Vector_t const& X = globalBoundingBox_m.lowerLeftCorner;
Vector_t const& Y = globalBoundingBox_m.upperRightCorner;
Vector_t delta = Y - X;
Vector_t dX(delta(0), 0, 0), dY(0, delta(1), 0), dZ(0, 0, delta(2));
std::vector<Vector_t> corners{X, X + dX, X + dX + dY, X + dY, Y, Y - dX, Y - dX - dY, Y - dY};
std::vector<std::vector<unsigned int>> paths{{0, 1, 2, 3}, {4, 5, 6, 7}, {0, 1, 7, 6}, {1, 2, 4, 7}, {2, 3, 5, 4}, {3, 0, 6, 5}};
std::ofstream out("boundingBox.gpl");
std::ofstream traj("trajectories.gpl");
for (std::vector<unsigned int> const& path: paths) {
for (unsigned int i = 0; i < 5; ++ i) {
Vector_t const& point = corners[path[i % 4]];
out << point(0) << "\t" << point(1) << "\t" << point(2) << std::endl;
}
out << std::endl;
}
}
if (maxDrift < minDistanceLocal &&
minDistanceLocal != std::numeric_limits<double>::max()) {
maxDrift = minDistanceLocal;
double OrbitThreader::computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
const Vector_t & position,
const Vector_t & direction) const {
if (elements.empty() ||
(elements.size() == 1 && (*elements.begin())->getType() == ElementBase::ElementType::DRIFT)) {
boost::optional<Vector_t> intersectionPoint = globalBoundingBox_m.getPointOfIntersection(position, direction);
double maxDrift = intersectionPoint ? euclidean_norm(intersectionPoint.get() - position): 10.0;
if (intersectionPoint) {
std::ofstream traj("trajectories.gpl", std::ios::app);
traj << position(0) << "\t" << position(1) << "\t" << position(2) << std::endl;
traj << intersectionPoint.get()(0) << "\t" << intersectionPoint.get()(1) << "\t" << intersectionPoint.get()(2) << std::endl;
traj << std::endl;
}
return maxDrift;
}
maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift);
return maxDrift;
return std::numeric_limits<double>::max();
}
\ No newline at end of file
src/Algorithms/OrbitThreader.h
View file @ 22bf204f
......@@ -20,7 +20,6 @@ public:
double maxDiffZBunch,
double t,
double dT,
size_t maxIntegSteps,
double zstop,
OpalBeamline &bl);
......@@ -48,8 +47,6 @@ private:
/// the time step
double dt_m;
/// the number of time steps to track
const size_t maxIntegSteps_m;
/// final position in path length
const double zstop_m;
......@@ -80,8 +77,8 @@ private:
std::multimap<std::string, elementPosition> elementRegistry_m;
void trackBack(double maxDrift);
void integrate(const IndexMap::value_t &activeSet, size_t maxSteps, double maxDrift = 10.0);
void trackBack();
void integrate(const IndexMap::value_t &activeSet, double maxDrift = 10.0);
bool containsCavity(const IndexMap::value_t &activeSet);
void autophaseCavities(const IndexMap::value_t &activeSet, const std::set<std::string> &visitedElements);
double getMaxDesignEnergy(const IndexMap::value_t &elementSet) const;
......@@ -89,8 +86,11 @@ private:
void registerElement(const IndexMap::value_t &elementSet, double, const Vector_t &r, const Vector_t &p);
void processElementRegister();
void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements);
void computeMaximalImplicitDrift2();
double computeMaximalImplicitDrift();
void computeBoundingBox();
// double computeMaximalImplicitDrift();
double computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
const Vector_t & position,
const Vector_t & direction) const;
};
inline
......
src/Algorithms/ParallelTTracker.cpp
View file @ 22bf204f
......@@ -193,7 +193,6 @@ void ParallelTTracker::execute() {
prepareSections();
double minTimeStep = stepSizes_m.getMinTimeStep();
unsigned long long totalNumSteps = stepSizes_m.getNumStepsFinestResolution();
itsOpalBeamline_m.activateElements();
......@@ -248,7 +247,6 @@ void ParallelTTracker::execute() {
-rmin(2),
itsBunch_m->getT(),
(back_track? -minTimeStep: minTimeStep),
totalNumSteps,
stepSizes_m.getFinalZStop() + 2 * rmax(2),
itsOpalBeamline_m);
......@@ -1390,4 +1388,4 @@ void ParallelTTracker::evenlyDistributeParticles() {
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End:
// End:
\ No newline at end of file
src/Classic/AbsBeamline/ElementBase.cpp
View file @ 22bf204f
......@@ -348,37 +348,51 @@ bool ElementBase::BoundingBox::isInside(const Vector_t & position) const {
boost::optional<Vector_t>
ElementBase::BoundingBox::getPointOfIntersection(const Vector_t & position,
const Vector_t & direction) const {
Vector_t relativePosition = position - lowerLeftCorner;
const Vector_t & direction) const {
Vector_t relativePosition = lowerLeftCorner - position;
Vector_t diagonal = upperRightCorner - lowerLeftCorner;
Vector_t normalizedDirection = direction / euclidean_norm(direction);
// *gmsg << "position: " << position << "\n"
// << "normalizedDirection: " << normalizedDirection << "\n"
// << "lowerLeftCorner: " << lowerLeftCorner << "\n"
// << "upperRightCorner: " << upperRightCorner << "\n"
// << "diagonal: " << diagonal << endl;
for (unsigned int i = 0; i < 2; ++ i) {
for (int i = -1; i < 2; i += 2) {
for (unsigned int d = 0; d < 3; ++ d) {
double projectedDirection = direction[d];
double projectedDirection = normalizedDirection[d];
// *gmsg << "\n" << "projectedDirection: " << projectedDirection << "\n";
if (std::abs(projectedDirection) < 1e-10) {
continue;
}
double distanceNearestPoint = relativePosition[d];
double tau = distanceNearestPoint / projectedDirection;
// *gmsg << "distanceNearestPoint: " << distanceNearestPoint << "\n"
// << "relativePosition: " << relativePosition << "\n"
// << "tau: " << tau << "\n";
if (tau < 0) {
continue;
}
Vector_t intersectionPoint = relativePosition + tau * direction;
Vector_t delta = tau * normalizedDirection;
Vector_t relativeIntersectionPoint = i * (relativePosition - delta);
double sign = 1 - 2 * i;
if (sign * intersectionPoint[(d + 1) % 3] < 0.0 ||
sign * intersectionPoint[(d + 1) % 3] > diagonal[(d + 1) % 3] ||
sign * intersectionPoint[(d + 2) % 3] < 0.0 ||
sign * intersectionPoint[(d + 2) % 3] > diagonal[(d + 2) % 3]) {
// *gmsg << "delta: " << delta << "\n"
// << "relativeIntersectionPoint: " << relativeIntersectionPoint << endl;
// *gmsg << sign * intersectionPoint[(d + 1) % 3] << "\t" << sign * relativeIntersectionPoint[(d + 2) % 3] << "\t";
if (relativeIntersectionPoint[(d + 1) % 3] < 0.0 ||
relativeIntersectionPoint[(d + 1) % 3] > diagonal[(d + 1) % 3] ||
relativeIntersectionPoint[(d + 2) % 3] < 0.0 ||
relativeIntersectionPoint[(d + 2) % 3] > diagonal[(d + 2) % 3]) {
continue;
}
return position + tau * direction;
return position + tau * normalizedDirection;
}
relativePosition = upperRightCorner - position;
}
// *gmsg << endl;
return boost::none;
}
......
src/Classic/AbsBeamline/RFCavity.cpp
View file @ 22bf204f
......@@ -742,4 +742,4 @@ void RFCavity::getElementDimensions(double &begin,
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End:
// End:
\ No newline at end of file
src/Classic/Fields/Astra1DDynamic_fast.cpp
View file @ 22bf204f
#include "Fields/Astra1DDynamic_fast.h"
#include "Utilities/GeneralClassicException.h"
#include "Utilities/Util.h"
#include "Utilities/OpalException.h"
#include "Physics/Physics.h"
#include <fstream>
#include <ios>
extern Inform *gmsg;
Astra1DDynamic_fast::Astra1DDynamic_fast(std::string aFilename):
Astra1D_fast(aFilename)
{
......@@ -82,11 +84,16 @@ bool Astra1DDynamic_fast::getFieldstrength(const Vector_t &R, Vector_t &E, Vecto
// do fourier interpolation in z-direction
const double RR2 = R(0) * R(0) + R(1) * R(1);
double ez = gsl_spline_eval(onAxisInterpolants_m[0], R(2) - zbegin_m, onAxisAccel_m[0]);
double ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2) - zbegin_m, onAxisAccel_m[1]);
double ezpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2) - zbegin_m, onAxisAccel_m[2]);
double ezppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2) - zbegin_m, onAxisAccel_m[3]);
double ez, ezp, ezpp, ezppp;
try {
ez = gsl_spline_eval(onAxisInterpolants_m[0], R(2) - zbegin_m, onAxisAccel_m[0]);
ezp = gsl_spline_eval(onAxisInterpolants_m[1], R(2) - zbegin_m, onAxisAccel_m[1]);
ezpp = gsl_spline_eval(onAxisInterpolants_m[2], R(2) - zbegin_m, onAxisAccel_m[2]);
ezppp = gsl_spline_eval(onAxisInterpolants_m[3], R(2) - zbegin_m, onAxisAccel_m[3]);
} catch (OpalException const& e) {
throw OpalException("Astra1DDynamice_fast::getFieldstrength",
"The requested interpolation point, " + std::to_string(R(2)) + " is out of range");
}
// expand the field off-axis
const double f = -(ezpp + ez * xlrep_m * xlrep_m) / 16.;
const double fp = -(ezppp + ezp * xlrep_m * xlrep_m) / 16.;
......@@ -115,6 +122,7 @@ void Astra1DDynamic_fast::getFieldDimensions(double &zBegin, double &zEnd) const
zBegin = zbegin_m;
zEnd = zend_m;
}
void Astra1DDynamic_fast::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/, double &/*yIni*/, double &/*yFinal*/, double &/*zIni*/, double &/*zFinal*/) const {}
void Astra1DDynamic_fast::swap()
......@@ -205,5 +213,4 @@ int Astra1DDynamic_fast::stripFileHeader(std::ifstream &file) {
interpretLine<double>(file, tmpDouble);
return accuracy;
}
}
\ No newline at end of file
src/Classic/Fields/Fieldmap.cpp
View file @ 22bf204f
......@@ -513,7 +513,7 @@ void Fieldmap::checkMap(unsigned int accuracy,
if (std::sqrt(error / ezSquare) > 1e-1 || maxDiff > 1e-1 * ezMax) {
lowResolutionWarning(std::sqrt(error / ezSquare), maxDiff / ezMax);
throw GeneralClassicException("Astra2DDynamic_fast::readMap()",
throw GeneralClassicException("Fieldmap::checkMap",
"Field map can't be reproduced properly with the given number of fourier components");
}
if (std::sqrt(error / ezSquare) > 1e-2 || maxDiff > 1e-2 * ezMax) {
......
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