Commit d2573049 authored by kraus's avatar kraus
Browse files

current status

parent 1f67683a
...@@ -86,6 +86,7 @@ OrbitThreader::OrbitThreader(const PartData &ref, ...@@ -86,6 +86,7 @@ OrbitThreader::OrbitThreader(const PartData &ref,
void OrbitThreader::execute() { void OrbitThreader::execute() {
double initialPathLength = pathLength_m; double initialPathLength = pathLength_m;
computeMaximalImplicitDrift2();
double maxDistance = computeMaximalImplicitDrift(); double maxDistance = computeMaximalImplicitDrift();
auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY); auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
...@@ -388,6 +389,22 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std:: ...@@ -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() { double OrbitThreader::computeMaximalImplicitDrift() {
FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY); FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
double maxDrift = 0.0; double maxDrift = 0.0;
...@@ -401,7 +418,6 @@ double OrbitThreader::computeMaximalImplicitDrift() { ...@@ -401,7 +418,6 @@ double OrbitThreader::computeMaximalImplicitDrift() {
FieldList::iterator it = allElements.begin(); FieldList::iterator it = allElements.begin();
const FieldList::iterator end = allElements.end(); const FieldList::iterator end = allElements.end();
for (; it != end; ++ it) { for (; it != end; ++ it) {
auto element1 = it->getElement(); auto element1 = it->getElement();
const auto &toEdge = element1->getCSTrafoGlobal2Local(); const auto &toEdge = element1->getCSTrafoGlobal2Local();
...@@ -459,4 +475,4 @@ double OrbitThreader::computeMaximalImplicitDrift() { ...@@ -459,4 +475,4 @@ double OrbitThreader::computeMaximalImplicitDrift() {
maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift); maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift);
return maxDrift; return maxDrift;
} }
\ No newline at end of file
...@@ -64,6 +64,8 @@ private: ...@@ -64,6 +64,8 @@ private:
std::ofstream logger_m; std::ofstream logger_m;
size_t loggingFrequency_m; size_t loggingFrequency_m;
ElementBase::BoundingBox globalBoundingBox_m;
struct elementPosition { struct elementPosition {
double startField_m; double startField_m;
double endField_m; double endField_m;
...@@ -87,6 +89,7 @@ private: ...@@ -87,6 +89,7 @@ private:
void registerElement(const IndexMap::value_t &elementSet, double, const Vector_t &r, const Vector_t &p); void registerElement(const IndexMap::value_t &elementSet, double, const Vector_t &r, const Vector_t &p);
void processElementRegister(); void processElementRegister();
void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements); void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements);
void computeMaximalImplicitDrift2();
double computeMaximalImplicitDrift(); double computeMaximalImplicitDrift();
}; };
......
...@@ -1682,12 +1682,12 @@ std::array<double,2> Bend2D::getExitFringeFieldLength() const { ...@@ -1682,12 +1682,12 @@ std::array<double,2> Bend2D::getExitFringeFieldLength() const {
return extFFL; return extFFL;
} }
BoundaryBox Bend2D::getBoundaryBoxInLabCoords() const { ElementBase::BoundingBox Bend2D::getBoundingBoxInLabCoords() const {
CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m; CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m; CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
double &x = aperture_m.second[0]; const double &x = aperture_m.second[0];
double &y = aperture_m.second[1]; const double &y = aperture_m.second[1];
std::vector<Vector_t> corners(8); std::vector<Vector_t> corners(8);
for (int i = -1; i < 2; i += 2) { for (int i = -1; i < 2; i += 2) {
...@@ -1706,11 +1706,11 @@ BoundaryBox Bend2D::getBoundaryBoxInLabCoords() const { ...@@ -1706,11 +1706,11 @@ BoundaryBox Bend2D::getBoundaryBoxInLabCoords() const {
if (bb.lowerLeftCorner(d) > corners[i](d)) { if (bb.lowerLeftCorner(d) > corners[i](d)) {
bb.lowerLeftCorner(d) = corners[i](d); bb.lowerLeftCorner(d) = corners[i](d);
} }
if (bb.upperLeftCorner(d) < corners[i](d)) { if (bb.upperRightCorner(d) < corners[i](d)) {
bb.upperLeftCorner(d) = corners[i](d); bb.upperRightCorner(d) = corners[i](d);
} }
} }
} }
return bb; return bb;
} }
\ No newline at end of file
...@@ -129,7 +129,7 @@ public: ...@@ -129,7 +129,7 @@ public:
// Used to create fringe fields in ThickTracker, (before edge[m], after edge[m]) // Used to create fringe fields in ThickTracker, (before edge[m], after edge[m])
std::array<double,2> getExitFringeFieldLength() const; std::array<double,2> getExitFringeFieldLength() const;
BoundaryBox getBoundaryBoxInLabCoords() const override; BoundingBox getBoundingBoxInLabCoords() const override;
protected: protected:
void setMessageHeader(const std::string & header); void setMessageHeader(const std::string & header);
double getStartField() const; double getStartField() const;
......
...@@ -308,37 +308,87 @@ void ElementBase::setCurrentSCoordinate(double s) { ...@@ -308,37 +308,87 @@ void ElementBase::setCurrentSCoordinate(double s) {
bool ElementBase::isInsideTransverse(const Vector_t &r) const bool ElementBase::isInsideTransverse(const Vector_t &r) const
{ {
double x = aperture_m.second[0]; const double &xLimit = aperture_m.second[0];
double y = aperture_m.second[1]; const double &yLimit = aperture_m.second[1];
double f = 1.0; double factor = 1.0;
if (aperture_m.first == CONIC_RECTANGULAR || if (aperture_m.first == CONIC_RECTANGULAR ||
aperture_m.first == CONIC_ELLIPTIC) { aperture_m.first == CONIC_ELLIPTICAL) {
Vector_t rRelativeToBegin = getEdgeToBegin().transformTo(r); Vector_t rRelativeToBegin = getEdgeToBegin().transformTo(r);
double fractionLength = rRelativeToBegin(2) / getElementLength(); double fractionLength = rRelativeToBegin(2) / getElementLength();
f = fractionLength * aperture_m.second[2]; factor = fractionLength * aperture_m.second[2];
} }
switch(aperture_m.first) { switch(aperture_m.first) {
case RECTANGULAR: 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: 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: 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: 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: default:
return false; 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 toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m; CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
double &x = aperture_m.second[0]; const double &x = aperture_m.second[0];
double &y = aperture_m.second[1]; const double &y = aperture_m.second[1];
double &f = aperture_m.second[2]; const double &f = aperture_m.second[2];
std::vector<Vector_t> corners(8); std::vector<Vector_t> corners(8);
for (int i = -1; i < 2; i += 2) { for (int i = -1; i < 2; i += 2) {
...@@ -353,12 +403,12 @@ BoundaryBox ElementBase::getBoundaryBoxInLabCoords() const { ...@@ -353,12 +403,12 @@ BoundaryBox ElementBase::getBoundaryBoxInLabCoords() const {
bb.lowerLeftCorner = bb.upperRightCorner = corners[0]; bb.lowerLeftCorner = bb.upperRightCorner = corners[0];
for (unsigned int i = 1; i < 8u; ++ i) { 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)) { if (bb.lowerLeftCorner(d) > corners[i](d)) {
bb.lowerLeftCorner(d) = corners[i](d); bb.lowerLeftCorner(d) = corners[i](d);
} }
if (bb.upperLeftCorner(d) < corners[i](d)) { if (bb.upperRightCorner(d) < corners[i](d)) {
bb.upperLeftCorner(d) = corners[i](d); bb.upperRightCorner(d) = corners[i](d);
} }
} }
} }
......
...@@ -71,6 +71,8 @@ ...@@ -71,6 +71,8 @@
#include "Algorithms/CoordinateSystemTrafo.h" #include "Algorithms/CoordinateSystemTrafo.h"
#include "Utilities/GeneralClassicException.h" #include "Utilities/GeneralClassicException.h"
#include <boost/optional.hpp>
#include <string> #include <string>
#include <queue> #include <queue>
...@@ -359,10 +361,26 @@ public: ...@@ -359,10 +361,26 @@ public:
void setRotationAboutZ(double rotation); void setRotationAboutZ(double rotation);
double getRotationAboutZ() const; 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: protected:
bool isInsideTransverse(const Vector_t &r, double f = 1) const; bool isInsideTransverse(const Vector_t &r) const;
// Sharable flag. // Sharable flag.
// If this flag is true, the element is always shared. // If this flag is true, the element is always shared.
......
...@@ -31,6 +31,7 @@ public: ...@@ -31,6 +31,7 @@ public:
return (fle.getLength() < 1.e-6); return (fle.getLength() < 1.e-6);
} }
ElementBase::BoundingBox getBoundingBoxInLabCoords() const;
CoordinateSystemTrafo getCoordTransformationTo() const ; CoordinateSystemTrafo getCoordTransformationTo() const ;
void setCoordTransformationTo(const CoordinateSystemTrafo &trafo); void setCoordTransformationTo(const CoordinateSystemTrafo &trafo);
...@@ -98,4 +99,8 @@ void ClassicField::fixPosition() { ...@@ -98,4 +99,8 @@ void ClassicField::fixPosition() {
element_m->fixPosition(); element_m->fixPosition();
} }
inline
ElementBase::BoundingBox ClassicField::getBoundingBoxInLabCoords() const {
return element_m->getBoundingBoxInLabCoords();
}
#endif // CLASSIC_FIELD_H #endif // CLASSIC_FIELD_H
...@@ -70,4 +70,4 @@ void OpalMarker::update() { ...@@ -70,4 +70,4 @@ void OpalMarker::update() {
// Transmit "unknown" attributes. // Transmit "unknown" attributes.
MarkerRep *mark = dynamic_cast<MarkerRep *>(getElement()); MarkerRep *mark = dynamic_cast<MarkerRep *>(getElement());
OpalElement::updateUnknown(mark); OpalElement::updateUnknown(mark);
} }
\ No newline at end of file
...@@ -193,7 +193,7 @@ void OpalRBend::update() { ...@@ -193,7 +193,7 @@ void OpalRBend::update() {
bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({hapert, gap, 1.0})); bend->setAperture(ElementBase::RECTANGULAR, std::vector<double>({hapert, gap, 1.0}));
} else { } else {
double gap = Attributes::getReal(itsAttr[GAP]); 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]) if(itsAttr[LENGTH])
......
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