From d25730490b1466d4faa9b2732ca8521bfb5ed4fb Mon Sep 17 00:00:00 2001
From: Christof Kraus <christof.j.kraus@gmail.com>
Date: Thu, 9 Jul 2020 22:43:15 +0200
Subject: [PATCH] current status
---
src/Algorithms/OrbitThreader.cpp | 20 +++++-
src/Algorithms/OrbitThreader.h | 3 +
src/Classic/AbsBeamline/Bend2D.cpp | 12 ++--
src/Classic/AbsBeamline/Bend2D.h | 2 +-
src/Classic/AbsBeamline/ElementBase.cpp | 82 ++++++++++++++++++++-----
src/Classic/AbsBeamline/ElementBase.h | 22 ++++++-
src/Classic/Utilities/ClassicField.h | 5 ++
src/Elements/OpalMarker.cpp | 2 +-
src/Elements/OpalRBend.cpp | 2 +-
9 files changed, 121 insertions(+), 29 deletions(-)
diff --git a/src/Algorithms/OrbitThreader.cpp b/src/Algorithms/OrbitThreader.cpp
index 3559d039d..f24dfc100 100644
--- a/src/Algorithms/OrbitThreader.cpp
+++ b/src/Algorithms/OrbitThreader.cpp
@@ -86,6 +86,7 @@ OrbitThreader::OrbitThreader(const PartData &ref,
void OrbitThreader::execute() {
double initialPathLength = pathLength_m;
+ computeMaximalImplicitDrift2();
double maxDistance = computeMaximalImplicitDrift();
auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
@@ -388,6 +389,22 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
}
}
+void OrbitThreader::computeMaximalImplicitDrift2() {
+ FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
+ FieldList::iterator it = allElements.begin();
+ const FieldList::iterator end = allElements.end();
+ globalBoundingBox_m.lowerLeftCorner = Vector_t(std::numeric_limits<double>::max());
+ globalBoundingBox_m.upperRightCorner = Vector_t(-std::numeric_limits<double>::max());
+
+ for (; it != end; ++ it) {
+ if (it->getElement()->getType() == ElementBase::MARKER) {
+ continue;
+ }
+ ElementBase::BoundingBox other = it->getBoundingBoxInLabCoords();
+ globalBoundingBox_m.getCombinedBoundingBox(other);
+ }
+}
+
double OrbitThreader::computeMaximalImplicitDrift() {
FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
double maxDrift = 0.0;
@@ -401,7 +418,6 @@ double OrbitThreader::computeMaximalImplicitDrift() {
FieldList::iterator it = allElements.begin();
const FieldList::iterator end = allElements.end();
-
for (; it != end; ++ it) {
auto element1 = it->getElement();
const auto &toEdge = element1->getCSTrafoGlobal2Local();
@@ -459,4 +475,4 @@ double OrbitThreader::computeMaximalImplicitDrift() {
maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift);
return maxDrift;
-}
+}
\ No newline at end of file
diff --git a/src/Algorithms/OrbitThreader.h b/src/Algorithms/OrbitThreader.h
index 010454157..09fccd121 100644
--- a/src/Algorithms/OrbitThreader.h
+++ b/src/Algorithms/OrbitThreader.h
@@ -64,6 +64,8 @@ private:
std::ofstream logger_m;
size_t loggingFrequency_m;
+ ElementBase::BoundingBox globalBoundingBox_m;
+
struct elementPosition {
double startField_m;
double endField_m;
@@ -87,6 +89,7 @@ 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();
};
diff --git a/src/Classic/AbsBeamline/Bend2D.cpp b/src/Classic/AbsBeamline/Bend2D.cpp
index 4f64a29f0..56da740a8 100644
--- a/src/Classic/AbsBeamline/Bend2D.cpp
+++ b/src/Classic/AbsBeamline/Bend2D.cpp
@@ -1682,12 +1682,12 @@ std::array<double,2> Bend2D::getExitFringeFieldLength() const {
return extFFL;
}
-BoundaryBox Bend2D::getBoundaryBoxInLabCoords() const {
+ElementBase::BoundingBox Bend2D::getBoundingBoxInLabCoords() const {
CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
- double &x = aperture_m.second[0];
- double &y = aperture_m.second[1];
+ const double &x = aperture_m.second[0];
+ const double &y = aperture_m.second[1];
std::vector<Vector_t> corners(8);
for (int i = -1; i < 2; i += 2) {
@@ -1706,11 +1706,11 @@ BoundaryBox Bend2D::getBoundaryBoxInLabCoords() const {
if (bb.lowerLeftCorner(d) > corners[i](d)) {
bb.lowerLeftCorner(d) = corners[i](d);
}
- if (bb.upperLeftCorner(d) < corners[i](d)) {
- bb.upperLeftCorner(d) = corners[i](d);
+ if (bb.upperRightCorner(d) < corners[i](d)) {
+ bb.upperRightCorner(d) = corners[i](d);
}
}
}
return bb;
-}
+}
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/Bend2D.h b/src/Classic/AbsBeamline/Bend2D.h
index d8b7638c2..502bf2c09 100644
--- a/src/Classic/AbsBeamline/Bend2D.h
+++ b/src/Classic/AbsBeamline/Bend2D.h
@@ -129,7 +129,7 @@ public:
// Used to create fringe fields in ThickTracker, (before edge[m], after edge[m])
std::array<double,2> getExitFringeFieldLength() const;
- BoundaryBox getBoundaryBoxInLabCoords() const override;
+ BoundingBox getBoundingBoxInLabCoords() const override;
protected:
void setMessageHeader(const std::string & header);
double getStartField() const;
diff --git a/src/Classic/AbsBeamline/ElementBase.cpp b/src/Classic/AbsBeamline/ElementBase.cpp
index 180a10e7e..540d0731c 100644
--- a/src/Classic/AbsBeamline/ElementBase.cpp
+++ b/src/Classic/AbsBeamline/ElementBase.cpp
@@ -308,37 +308,87 @@ void ElementBase::setCurrentSCoordinate(double s) {
bool ElementBase::isInsideTransverse(const Vector_t &r) const
{
- double x = aperture_m.second[0];
- double y = aperture_m.second[1];
- double f = 1.0;
+ const double &xLimit = aperture_m.second[0];
+ const double &yLimit = aperture_m.second[1];
+ double factor = 1.0;
if (aperture_m.first == CONIC_RECTANGULAR ||
- aperture_m.first == CONIC_ELLIPTIC) {
+ aperture_m.first == CONIC_ELLIPTICAL) {
Vector_t rRelativeToBegin = getEdgeToBegin().transformTo(r);
double fractionLength = rRelativeToBegin(2) / getElementLength();
- f = fractionLength * aperture_m.second[2];
+ factor = fractionLength * aperture_m.second[2];
}
switch(aperture_m.first) {
case RECTANGULAR:
- return (std::abs(r[0]) < aperture_m.second[0] && std::abs(r[1]) < aperture_m.second[1]);
+ return (std::abs(r[0]) < xLimit && std::abs(r[1]) < yLimit);
case ELLIPTICAL:
- return (std::pow(r[0] / aperture_m.second[0], 2) + std::pow(r[1] / aperture_m.second[1], 2) < 1.0);
+ return (std::pow(r[0] / xLimit, 2) + std::pow(r[1] / yLimit, 2) < 1.0);
case CONIC_RECTANGULAR:
- return (std::abs(r[0]) < f * aperture_m.second[0] && std::abs(r[1]) < f * aperture_m.second[1]);
+ return (std::abs(r[0]) < factor * xLimit && std::abs(r[1]) < factor * yLimit);
case CONIC_ELLIPTICAL:
- return (std::pow(r[0] / (f * aperture_m.second[0]), 2) + std::pow(r[1] / (f * aperture_m.second[1]), 2) < 1.0);
+ return (std::pow(r[0] / (factor * xLimit), 2) + std::pow(r[1] / (factor * yLimit), 2) < 1.0);
default:
return false;
}
}
-BoundaryBox ElementBase::getBoundaryBoxInLabCoords() const {
+bool ElementBase::BoundingBox::isInside(const Vector_t & position) const {
+ Vector_t relativePosition = position - lowerLeftCorner;
+ Vector_t diagonal = upperRightCorner - lowerLeftCorner;
+
+ for (unsigned int d = 0; d < 3; ++ d) {
+ if (relativePosition[d] < 0.0 ||
+ relativePosition[d] > diagonal[d]) {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+boost::optional<Vector_t>
+ElementBase::BoundingBox::getPointOfIntersection(const Vector_t & position,
+ const Vector_t & direction) const {
+ Vector_t relativePosition = position - lowerLeftCorner;
+ Vector_t diagonal = upperRightCorner - lowerLeftCorner;
+
+ for (unsigned int i = 0; i < 2; ++ i) {
+ for (unsigned int d = 0; d < 3; ++ d) {
+ double projectedDirection = direction[d];
+ if (std::abs(projectedDirection) < 1e-10) {
+ continue;
+ }
+
+ double distanceNearestPoint = relativePosition[d];
+ double tau = distanceNearestPoint / projectedDirection;
+ if (tau < 0) {
+ continue;
+ }
+ Vector_t intersectionPoint = relativePosition + tau * direction;
+
+ 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]) {
+ continue;
+ }
+
+ return position + tau * direction;
+ }
+ relativePosition = upperRightCorner - position;
+ }
+
+ return boost::none;
+}
+
+ElementBase::BoundingBox ElementBase::getBoundingBoxInLabCoords() const {
CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
- double &x = aperture_m.second[0];
- double &y = aperture_m.second[1];
- double &f = aperture_m.second[2];
+ const double &x = aperture_m.second[0];
+ const double &y = aperture_m.second[1];
+ const double &f = aperture_m.second[2];
std::vector<Vector_t> corners(8);
for (int i = -1; i < 2; i += 2) {
@@ -353,12 +403,12 @@ BoundaryBox ElementBase::getBoundaryBoxInLabCoords() const {
bb.lowerLeftCorner = bb.upperRightCorner = corners[0];
for (unsigned int i = 1; i < 8u; ++ i) {
- for (unsigned int d; d < 3u; ++ d) {
+ for (unsigned int d = 0; d < 3u; ++ d) {
if (bb.lowerLeftCorner(d) > corners[i](d)) {
bb.lowerLeftCorner(d) = corners[i](d);
}
- if (bb.upperLeftCorner(d) < corners[i](d)) {
- bb.upperLeftCorner(d) = corners[i](d);
+ if (bb.upperRightCorner(d) < corners[i](d)) {
+ bb.upperRightCorner(d) = corners[i](d);
}
}
}
diff --git a/src/Classic/AbsBeamline/ElementBase.h b/src/Classic/AbsBeamline/ElementBase.h
index aaa2fd5a4..4c3f6d50c 100644
--- a/src/Classic/AbsBeamline/ElementBase.h
+++ b/src/Classic/AbsBeamline/ElementBase.h
@@ -71,6 +71,8 @@
#include "Algorithms/CoordinateSystemTrafo.h"
#include "Utilities/GeneralClassicException.h"
+#include <boost/optional.hpp>
+
#include <string>
#include <queue>
@@ -359,10 +361,26 @@ public:
void setRotationAboutZ(double rotation);
double getRotationAboutZ() const;
- virtual BoundaryBox getBoundaryBoxInLabCoords() const;
+ struct BoundingBox {
+ Vector_t lowerLeftCorner;
+ Vector_t upperRightCorner;
+
+ void getCombinedBoundingBox(const BoundingBox & other) {
+ for (unsigned int d = 0; d < 3; ++ d) {
+ lowerLeftCorner[d] = std::min(lowerLeftCorner[d], other.lowerLeftCorner[d]);
+ upperRightCorner[d] = std::max(upperRightCorner[d], other.upperRightCorner[d]);
+ }
+ }
+
+ bool isInside(const Vector_t &) const;
+ boost::optional<Vector_t> getPointOfIntersection(const Vector_t & position,
+ const Vector_t & direction) const;
+ };
+
+ virtual BoundingBox getBoundingBoxInLabCoords() const;
protected:
- bool isInsideTransverse(const Vector_t &r, double f = 1) const;
+ bool isInsideTransverse(const Vector_t &r) const;
// Sharable flag.
// If this flag is true, the element is always shared.
diff --git a/src/Classic/Utilities/ClassicField.h b/src/Classic/Utilities/ClassicField.h
index 180dd0256..60028adfc 100644
--- a/src/Classic/Utilities/ClassicField.h
+++ b/src/Classic/Utilities/ClassicField.h
@@ -31,6 +31,7 @@ public:
return (fle.getLength() < 1.e-6);
}
+ ElementBase::BoundingBox getBoundingBoxInLabCoords() const;
CoordinateSystemTrafo getCoordTransformationTo() const ;
void setCoordTransformationTo(const CoordinateSystemTrafo &trafo);
@@ -98,4 +99,8 @@ void ClassicField::fixPosition() {
element_m->fixPosition();
}
+inline
+ElementBase::BoundingBox ClassicField::getBoundingBoxInLabCoords() const {
+ return element_m->getBoundingBoxInLabCoords();
+}
#endif // CLASSIC_FIELD_H
diff --git a/src/Elements/OpalMarker.cpp b/src/Elements/OpalMarker.cpp
index 1bd5eb58b..b7fc5cc0f 100644
--- a/src/Elements/OpalMarker.cpp
+++ b/src/Elements/OpalMarker.cpp
@@ -70,4 +70,4 @@ void OpalMarker::update() {
// Transmit "unknown" attributes.
MarkerRep *mark = dynamic_cast<MarkerRep *>(getElement());
OpalElement::updateUnknown(mark);
-}
+}
\ No newline at end of file
diff --git a/src/Elements/OpalRBend.cpp b/src/Elements/OpalRBend.cpp
index af1f3327e..90fecb7d5 100644
--- a/src/Elements/OpalRBend.cpp
+++ b/src/Elements/OpalRBend.cpp
@@ -193,7 +193,7 @@ void OpalRBend::update() {
bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({hapert, gap, 1.0}));
} else {
double gap = Attributes::getReal(itsAttr[GAP]);
- bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({0.5, gap, 1.0});
+ bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({0.5, gap, 1.0}));
}
if(itsAttr[LENGTH])
--
GitLab