diff --git a/src/Algorithms/IndexMap.cpp b/src/Algorithms/IndexMap.cpp
index 330b424c93a2c0054880268a2c8c5e56d9ca8e1a..a825cb9b7244aea4f53f1cc8155fc63b6c737cf8 100644
--- a/src/Algorithms/IndexMap.cpp
+++ b/src/Algorithms/IndexMap.cpp
@@ -117,11 +117,19 @@ void IndexMap::add(key_t::first_type initialS, key_t::second_type finalS, const
}
}
-void IndexMap::tidyUp() {
+void IndexMap::tidyUp(double zstop) {
map_t::reverse_iterator rit = mapRange2Element_m.rbegin();
- if (rit != mapRange2Element_m.rend() && (*rit).second.size() == 0) {
+ if (rit != mapRange2Element_m.rend() &&
+ (*rit).second.size() == 0 &&
+ zstop > (*rit).first.first &&
+ zstop < (*rit).first.second) {
+
+ key_t key((*rit).first.first, zstop);
+ value_t val;
+
mapRange2Element_m.erase(std::next(rit).base());
+ mapRange2Element_m.insert(std::pair<key_t, value_t>(key, val));
}
}
@@ -135,17 +143,18 @@ enum elements {
SOLENOID,
RFCAVITY,
MONITOR,
+ OTHER,
SIZE
};
-void IndexMap::saveSDDS(double startS) const {
+void IndexMap::saveSDDS(double initialPathLength) const {
auto opal = OpalData::getInstance();
if (opal->isOptimizerRun()) return;
std::string fileName("data/" + OpalData::getInstance()->getInputBasename() + "_ElementPositions.sdds");
std::ofstream sdds;
if (OpalData::getInstance()->hasPriorTrack() && boost::filesystem::exists(fileName)) {
- Util::rewindLinesSDDS(fileName, startS, false);
+ Util::rewindLinesSDDS(fileName, initialPathLength, false);
sdds.open(fileName, std::ios::app);
} else {
sdds.open(fileName);
@@ -222,6 +231,12 @@ void IndexMap::saveSDDS(double startS) const {
<< indent << "units=1, \n"
<< indent << "description=\"10 monitor present\" \n"
<< "&end\n";
+ sdds << "&column \n"
+ << indent << "name=other, \n"
+ << indent << "type=float, \n"
+ << indent << "units=1, \n"
+ << indent << "description=\"11 other element present\" \n"
+ << "&end\n";
sdds << "&column \n"
<< indent << "name=element_names, \n"
<< indent << "type=string, \n"
@@ -238,7 +253,7 @@ void IndexMap::saveSDDS(double startS) const {
std::vector<std::vector<int> > allItems(SIZE);
std::vector<double> allPositions;
std::vector<std::string> allNames;
- std::vector<double> typeMultipliers = {3.3333e-1, 1.0, 0.5, 0.25, 1.0, 1.0, 1.0, 1.0, 1.0};
+ std::vector<double> typeMultipliers = {3.3333e-1, 1.0, 0.5, 0.25, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
unsigned int counter = 0;
auto mapIti = mapRange2Element_m.begin();
@@ -298,7 +313,8 @@ void IndexMap::saveSDDS(double startS) const {
items[MONITOR] = 1;
break;
default:
- continue;
+ items[OTHER] = 1;
+ break;
}
names += element->getName() + ", ";
@@ -328,6 +344,26 @@ void IndexMap::saveSDDS(double startS) const {
if (counter == 0) return;
+ if (allPositions.front() > initialPathLength) {
+ {
+ auto tmp = allPositions;
+ allPositions = std::vector<double>(4, initialPathLength);
+ allPositions.insert(allPositions.end(), tmp.begin(), tmp.end());
+ }
+
+ {
+ auto tmp = allNames;
+ allNames = std::vector<std::string>(4, "");
+ allNames.insert(allNames.end(), tmp.begin(), tmp.end());
+ }
+
+ for (unsigned int i = 0; i < SIZE; ++ i) {
+ auto tmp = allItems[i];
+ allItems[i] = std::vector<int>(4, 0);
+ allItems[i].insert(allItems[i].end(), tmp.begin(), tmp.end());
+ }
+ }
+
const unsigned int totalSize = counter;
for (unsigned int i = 0; i < SIZE; ++ i) {
for (unsigned int j = 0; j < totalSize - 1; ++ j) {
diff --git a/src/Algorithms/IndexMap.h b/src/Algorithms/IndexMap.h
index 78ff89e9a8b4bc9a3a759a396525839f156a5327..2b24769c4cb05bbe9368ec3d2f95fa76bcfe6627 100644
--- a/src/Algorithms/IndexMap.h
+++ b/src/Algorithms/IndexMap.h
@@ -23,7 +23,7 @@ public:
value_t query(key_t::first_type s, key_t::second_type ds);
- void tidyUp();
+ void tidyUp(double zstop);
void print(std::ostream&) const;
void saveSDDS(double startS) const;
@@ -92,4 +92,4 @@ Inform& operator<< (Inform &out, const IndexMap &im) {
return out;
}
-#endif
+#endif
\ No newline at end of file
diff --git a/src/Algorithms/OrbitThreader.cpp b/src/Algorithms/OrbitThreader.cpp
index a66bbb56b4a358df0446d1bc5ee17c9cb73097ab..d60b75b63d975219c6ff49e78b08c71a467acfe6 100644
--- a/src/Algorithms/OrbitThreader.cpp
+++ b/src/Algorithms/OrbitThreader.cpp
@@ -4,6 +4,7 @@
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/BendBase.h"
#include "AbsBeamline/TravelingWave.h"
+#include "BeamlineCore/MarkerRep.h"
#include "AbstractObjects/OpalData.h"
#include "BasicActions/Option.h"
@@ -132,7 +133,7 @@ void OrbitThreader::execute() {
} while (errorFlag_m != HITMATERIAL &&
errorFlag_m != EOL);
- imap_m.tidyUp();
+ imap_m.tidyUp(zstop_m);
*gmsg << level1 << "\n" << imap_m << endl;
imap_m.saveSDDS(initialPathLength);
@@ -387,6 +388,14 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
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();
@@ -414,7 +423,9 @@ double OrbitThreader::computeMaximalImplicitDrift() {
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 ||
@@ -427,15 +438,22 @@ double OrbitThreader::computeMaximalImplicitDrift() {
double distance = euclidean_norm(begin2 - end1);
double directionProjection = dot(directionEnd, directionBegin);
- if (directionProjection > 0.999 && minDistanceLocal > distance) {
+ bool overlapping = dot(begin2 - end1, directionBegin) < 0.0? true: false;
+
+ if (!overlapping &&
+ directionProjection > 0.999 &&
+ minDistanceLocal > distance) {
minDistanceLocal = distance;
}
}
- if (maxDrift < minDistanceLocal) {
+ if (maxDrift < minDistanceLocal &&
+ minDistanceLocal != std::numeric_limits<double>::max()) {
maxDrift = minDistanceLocal;
}
}
+ maxDrift = std::min(maxIntegSteps_m * dt_m * Physics::c, maxDrift);
+
return maxDrift;
}
\ No newline at end of file