diff --git a/src/Algorithms/CavityAutophaser.cpp b/src/Algorithms/CavityAutophaser.cpp
index 9fa6b29ecacc0e848295ab56a1fcd8a27e7dc419..11a9d4e9c1e38361cb8ef22f0bbd2c721d8dc592 100644
--- a/src/Algorithms/CavityAutophaser.cpp
+++ b/src/Algorithms/CavityAutophaser.cpp
@@ -19,8 +19,7 @@ CavityAutophaser::CavityAutophaser(const PartData &ref,
itsCavity_m(cavity)
{
double zbegin = 0.0, zend = 0.0;
- PartBunchBase<double, 3> *fakeBunch = NULL;
- cavity->initialise(fakeBunch, zbegin, zend);
+ cavity->getDimensions(zbegin, zend);
initialR_m = Vector_t(0, 0, zbegin);
}
diff --git a/src/Algorithms/OrbitThreader.cpp b/src/Algorithms/OrbitThreader.cpp
index f24dfc1004356e3af2e5cad51f441a9a1445aae8..c340a5374bc007e39c4696c3e2c4df8531019eba 100644
--- a/src/Algorithms/OrbitThreader.cpp
+++ b/src/Algorithms/OrbitThreader.cpp
@@ -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
diff --git a/src/Algorithms/OrbitThreader.h b/src/Algorithms/OrbitThreader.h
index 09fccd121538487fba4fb8d125be33c91c0e0eb5..c56b4a7a4b9feca0ffeaf993d68b9caab65c1842 100644
--- a/src/Algorithms/OrbitThreader.h
+++ b/src/Algorithms/OrbitThreader.h
@@ -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
diff --git a/src/Algorithms/ParallelTTracker.cpp b/src/Algorithms/ParallelTTracker.cpp
index f38da15e80bd5e36a201c9a9637a7aeaee16c12c..49f56f3e87fbec86534b53c7719095938fa338d9 100644
--- a/src/Algorithms/ParallelTTracker.cpp
+++ b/src/Algorithms/ParallelTTracker.cpp
@@ -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
diff --git a/src/Classic/AbsBeamline/ElementBase.cpp b/src/Classic/AbsBeamline/ElementBase.cpp
index 540d0731c813156f72a8bd7f68bbdbdc9d53ca28..2e7336b231c0a850b43e56f691ffefad2c36e6f4 100644
--- a/src/Classic/AbsBeamline/ElementBase.cpp
+++ b/src/Classic/AbsBeamline/ElementBase.cpp
@@ -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;
}
diff --git a/src/Classic/AbsBeamline/RFCavity.cpp b/src/Classic/AbsBeamline/RFCavity.cpp
index b8ebb79459043ae993e9c3b3b05712ad5183e541..5f08002c23f702ca573d735971b26e5c847e2fbe 100644
--- a/src/Classic/AbsBeamline/RFCavity.cpp
+++ b/src/Classic/AbsBeamline/RFCavity.cpp
@@ -740,4 +740,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
diff --git a/src/Classic/Fields/Astra1DDynamic_fast.cpp b/src/Classic/Fields/Astra1DDynamic_fast.cpp
index dcfab8e1e462f9c9353a9187764ee346d75b1281..dc57e81a44bfbbbee56d829cfd5cf808909a8ae8 100644
--- a/src/Classic/Fields/Astra1DDynamic_fast.cpp
+++ b/src/Classic/Fields/Astra1DDynamic_fast.cpp
@@ -1,11 +1,13 @@
#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
diff --git a/src/Classic/Fields/Fieldmap.cpp b/src/Classic/Fields/Fieldmap.cpp
index 3a420bfe468bb20c4c0f9b0bf624737ee55b4e45..b8eb3e6b4e75e6c8a58ca56ef90b2d4e931f2ad4 100644
--- a/src/Classic/Fields/Fieldmap.cpp
+++ b/src/Classic/Fields/Fieldmap.cpp
@@ -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) {