diff --git a/src/Algorithms/CMakeLists.txt b/src/Algorithms/CMakeLists.txt
index ebe21b7f4c87d1121ce38d73c3b214b5de43c5d2..25315a638f0f4f9c25a646fbc75a83f304368436 100644
--- a/src/Algorithms/CMakeLists.txt
+++ b/src/Algorithms/CMakeLists.txt
@@ -10,9 +10,10 @@ set (_SRCS
MPSplitIntegrator.cpp
MultiBunchHandler.cpp
OrbitThreader.cpp
- ParallelTTracker.cpp
ParallelCyclotronTracker.cpp
ParallelSliceTracker.cpp
+ ParallelTTracker.cpp
+ StepSizeConfig.cpp
ThickMapper.cpp
ThickTracker.cpp
TransportMapper.cpp
@@ -51,6 +52,7 @@ set (HDRS
ParallelCyclotronTracker.h
ParallelSliceTracker.h
ParallelTTracker.h
+ StepSizeConfig.h
ThickMapper.h
ThickTracker.h
TransportMapper.h
@@ -68,4 +70,4 @@ set (HDRS
bet/math/svdfit.h
)
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Algorithms")
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Algorithms")
\ No newline at end of file
diff --git a/src/Algorithms/IndexMap.cpp b/src/Algorithms/IndexMap.cpp
index 317b30e3c1058c4064e312523ea9b64f40d7fae8..4b2e9fc8b35a84fd5cc482f1bc337f323b9fe5e4 100644
--- a/src/Algorithms/IndexMap.cpp
+++ b/src/Algorithms/IndexMap.cpp
@@ -51,7 +51,7 @@ void IndexMap::print(std::ostream &out) const {
}
IndexMap::value_t IndexMap::query(key_t::first_type s, key_t::second_type ds) {
- const double lowerLimit = (ds < s? s - ds: 0);
+ const double lowerLimit = s - ds;//(ds < s? s - ds: 0);
const double upperLimit = std::min(totalPathLength_m, s + ds);
value_t elementSet;
@@ -88,7 +88,11 @@ IndexMap::value_t IndexMap::query(key_t::first_type s, key_t::second_type ds) {
}
void IndexMap::add(key_t::first_type initialS, key_t::second_type finalS, const value_t &val) {
+ if (initialS > finalS) {
+ std::swap(initialS, finalS);
+ }
key_t key(initialS, finalS * oneMinusEpsilon_m);
+
mapRange2Element_m.insert(std::pair<key_t, value_t>(key, val));
totalPathLength_m = (*mapRange2Element_m.rbegin()).first.second;
diff --git a/src/Algorithms/OrbitThreader.cpp b/src/Algorithms/OrbitThreader.cpp
index 99fb846c5e40ba5a1e3e2c3023e7a5a950f6622f..fcc672913610d412a5015f6aebe57bceb5368dc9 100644
--- a/src/Algorithms/OrbitThreader.cpp
+++ b/src/Algorithms/OrbitThreader.cpp
@@ -72,12 +72,12 @@ OrbitThreader::OrbitThreader(const PartData &ref,
logger_m.open(fileName, std::ios_base::app);
}
- loggingFrequency_m = std::max(1.0, floor(1e-11 / dt_m + 0.5));
+ loggingFrequency_m = std::max(1.0, floor(1e-11 / std::abs(dt_m) + 0.5));
} else {
loggingFrequency_m = std::numeric_limits<size_t>::max();
}
- dZ_m = std::min(pathLength_m, std::max(0.0, maxDiffZBunch));
+ distTrackBack_m = std::min(pathLength_m, std::max(0.0, maxDiffZBunch));
}
void OrbitThreader::execute() {
@@ -110,7 +110,9 @@ void OrbitThreader::execute() {
registerElement(elementSet, initialS, initialR, initialP);
if (errorFlag_m == HITMATERIAL) {
- pathLength_m += 1.0;
+ // Shouldn't be reached because reference particle
+ // isn't stopped by collimators
+ pathLength_m += std::copysign(1.0, dt_m);
}
imap_m.add(initialS, pathLength_m, elementSet);
@@ -144,7 +146,6 @@ void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxStep
static size_t step = 0;
CoordinateSystemTrafo labToBeamline = itsOpalBeamline_m.getCSTrafoLab2Local();
const double oldPathLength = pathLength_m;
- double stepLength;
Vector_t nextR;
do {
errorFlag_m = EVERYTHINGFINE;
@@ -211,15 +212,13 @@ void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxStep
nextR *= Physics::c * dt_m;
if ((activeSet.size() == 0 && std::abs(pathLength_m - oldPathLength) > maxDrift) ||
- (activeSet.size() > 0 && pathLength_m > zstop_m)) {
+ (activeSet.size() > 0 && dt_m * pathLength_m > dt_m * zstop_m) ||
+ (pathLength_m < 0.0 && dt_m < 0.0)) {
errorFlag_m = EOL;
return;
}
- stepLength = std::abs(dt_m) * euclidean_norm(p_m) / sqrt(dot(p_m, p_m) + 1) * Physics::c;
- } while ((dt_m > 0.0 ||
- euclidean_norm(labToBeamline.transformTo(r_m)) > stepLength) &&
- (activeSet == itsOpalBeamline_m.getElements(nextR)));
+ } while (activeSet == itsOpalBeamline_m.getElements(nextR));
}
bool OrbitThreader::containsCavity(const IndexMap::value_t &activeSet) {
@@ -286,8 +285,8 @@ void OrbitThreader::trackBack(double maxDrift) {
integrator_m.push(nextR, p_m, dt_m);
nextR *= Physics::c * dt_m;
- maxDrift = std::min(maxDrift, dZ_m);
- while (initialPathLength - pathLength_m < dZ_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);
@@ -453,7 +452,7 @@ double OrbitThreader::computeMaximalImplicitDrift() {
}
}
- maxDrift = std::min(maxIntegSteps_m * dt_m * Physics::c, maxDrift);
+ 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 af8247c1eb52b50fd57416ac1a73775afbdc5d5d..010454157a710711b71dace66c0f60fbe9ea34ab 100644
--- a/src/Algorithms/OrbitThreader.h
+++ b/src/Algorithms/OrbitThreader.h
@@ -34,13 +34,23 @@ public:
IndexMap::value_t getTouchingElements(const std::pair<double, double> &range);
private:
+ /// position of reference particle in lab coordinates
Vector_t r_m;
+ /// momentum of reference particle
Vector_t p_m;
+ /// position of reference particle in path length
double pathLength_m;
- double dZ_m;
+ /// distance to track back before tracking forward
+ /// (length of bunch but not beyond cathode)
+ double distTrackBack_m;
+ /// the simulated time
double time_m;
+ /// 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;
OpalBeamline &itsOpalBeamline_m;
diff --git a/src/Algorithms/ParallelTTracker.cpp b/src/Algorithms/ParallelTTracker.cpp
index 486bef3650742cb2d24a8b4f36acf8e381abd0be..380b1ac76c22a62c4f9b9f90e8cf9587fceeba59 100644
--- a/src/Algorithms/ParallelTTracker.cpp
+++ b/src/Algorithms/ParallelTTracker.cpp
@@ -64,10 +64,7 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
wakeFunction_m(NULL),
pathLength_m(0.0),
zstart_m(0.0),
- zStop_m(),
dtCurrentTrack_m(0.0),
- dtAllTracks_m(),
- localTrackSteps_m(),
minStepforReBin_m(-1),
minBinEmitted_m(std::numeric_limits<size_t>::max()),
repartFreq_m(-1),
@@ -120,16 +117,13 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
particleMatterStatus_m(false),
totalParticlesInSimulation_m(0)
{
- for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
- localTrackSteps_m.push(*it);
- }
- for (std::vector<double>::const_iterator it = dt.begin(); it != dt.end(); ++ it) {
- dtAllTracks_m.push(*it);
- }
- for (std::vector<double>::const_iterator it = zstop.begin(); it != zstop.end(); ++ it) {
- zStop_m.push(*it);
+ for (unsigned int i = 0; i < zstop.size(); ++ i) {
+ stepSizes_m.push_back(dt[i], zstop[i], maxSteps[i]);
}
+ stepSizes_m.sortAscendingZStop();
+ stepSizes_m.resetIterator();
+
#ifdef OPAL_DKS
if (IpplInfo::DKSEnabled)
setupDKS();
@@ -202,6 +196,8 @@ void ParallelTTracker::execute() {
const double globalTimeShift = itsBunch_m->weHaveEnergyBins()? OpalData::getInstance()->getGlobalPhaseShift(): 0.0;
OpalData::getInstance()->setGlobalPhaseShift(0.0);
+ // the time step needs to be positive in the setup
+ itsBunch_m->setdT(std::abs(itsBunch_m->getdT()));
dtCurrentTrack_m = itsBunch_m->getdT();
evenlyDistributeParticles();
@@ -212,20 +208,8 @@ void ParallelTTracker::execute() {
prepareSections();
- std::queue<double> timeStepSizes(dtAllTracks_m);
- std::queue<unsigned long long> numSteps(localTrackSteps_m);
- double minTimeStep = timeStepSizes.front();
- unsigned long long totalNumSteps = 0;
- while (timeStepSizes.size() > 0) {
- if (minTimeStep > timeStepSizes.front()) {
- totalNumSteps = std::ceil(totalNumSteps * minTimeStep / timeStepSizes.front());
- minTimeStep = timeStepSizes.front();
- }
- totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);
-
- numSteps.pop();
- timeStepSizes.pop();
- }
+ double minTimeStep = stepSizes_m.getMinTimeStep();
+ unsigned long long totalNumSteps = stepSizes_m.getNumStepsFinestResolution();
itsOpalBeamline_m.activateElements();
@@ -259,19 +243,29 @@ void ParallelTTracker::execute() {
}
}
+ stepSizes_m.advanceToPos(zstart_m);
+ if (back_track) {
+ itsBunch_m->setdT(-std::abs(itsBunch_m->getdT()));
+ stepSizes_m.reverseDirection();
+ if (pathLength_m < stepSizes_m.getZStop()) {
+ ++ stepSizes_m;
+ }
+ }
+
Vector_t rmin(0.0), rmax(0.0);
if (itsBunch_m->getTotalNum() > 0) {
itsBunch_m->get_bounds(rmin, rmax);
}
+
OrbitThreader oth(itsReference,
itsBunch_m->RefPartR_m,
itsBunch_m->RefPartP_m,
pathLength_m,
-rmin(2),
itsBunch_m->getT(),
- minTimeStep,
+ (back_track? -minTimeStep: minTimeStep),
totalNumSteps,
- zStop_m.back() + 2 * rmax(2),
+ stepSizes_m.getFinalZStop() + 2 * rmax(2),
itsOpalBeamline_m);
oth.execute();
@@ -298,7 +292,7 @@ void ParallelTTracker::execute() {
*gmsg << level1
<< "Executing ParallelTTracker, initial dt= " << Util::getTimeString(itsBunch_m->getdT()) << ";\n"
- << "max integration steps " << getMaxSteps(localTrackSteps_m) << ", next step= " << step << endl;
+ << "max integration steps " << stepSizes_m.getMaxSteps() << ", next step= " << step << endl;
setOptionalVariables();
@@ -307,17 +301,20 @@ void ParallelTTracker::execute() {
allocateDeviceMemory();
#endif
- // loggingFrequency_m = floor(1e-11/itsBunch_m->getdT() + 0.5);
globalEOL_m = false;
wakeStatus_m = false;
deletedParticles_m = false;
OpalData::getInstance()->setInPrepState(false);
- while (localTrackSteps_m.size() > 0) {
- localTrackSteps_m.front() += step;
- dtCurrentTrack_m = dtAllTracks_m.front();
- changeDT();
-
- for (; step < localTrackSteps_m.front(); ++step) {
+ while (!stepSizes_m.reachedEnd()) {
+ unsigned long long trackSteps = stepSizes_m.getNumSteps() + step;
+ dtCurrentTrack_m = stepSizes_m.getdT();
+ changeDT(back_track);
+
+ for (; step < trackSteps; ++ step) {
+ Vector_t rmin(0.0), rmax(0.0);
+ if (itsBunch_m->getTotalNum() > 0) {
+ itsBunch_m->get_bounds(rmin, rmax);
+ }
timeIntegration1(pusher);
@@ -326,9 +323,9 @@ void ParallelTTracker::execute() {
computeSpaceChargeFields(step);
- selectDT();
+ selectDT(back_track);
emitParticles(step);
- selectDT();
+ selectDT(back_track);
computeExternalFields(oth);
@@ -336,7 +333,8 @@ void ParallelTTracker::execute() {
itsBunch_m->incrementT();
- if (itsBunch_m->getT() > 0.0) {
+ if (itsBunch_m->getT() > 0.0 ||
+ itsBunch_m->getdT() < 0.0) {
updateReference(pusher);
}
@@ -344,7 +342,6 @@ void ParallelTTracker::execute() {
evenlyDistributeParticles();
deletedParticles_m = false;
}
-
itsBunch_m->set_sPos(pathLength_m);
if (hasEndOfLineReached()) break;
@@ -356,17 +353,16 @@ void ParallelTTracker::execute() {
itsBunch_m->incTrackSteps();
double beta = euclidean_norm(itsBunch_m->RefPartP_m / Util::getGamma(itsBunch_m->RefPartP_m));
- double driftPerTimeStep = itsBunch_m->getdT() * Physics::c * beta;
- if (std::abs(zStop_m.front() - pathLength_m) < 0.5 * driftPerTimeStep)
- localTrackSteps_m.front() = step;
+ double driftPerTimeStep = std::abs(itsBunch_m->getdT()) * Physics::c * beta;
+ if (std::abs(stepSizes_m.getZStop() - pathLength_m) < 0.5 * driftPerTimeStep) {
+ break;
+ }
}
if (globalEOL_m)
break;
- dtAllTracks_m.pop();
- localTrackSteps_m.pop();
- zStop_m.pop();
+ ++ stepSizes_m;
}
itsBunch_m->set_sPos(pathLength_m);
@@ -465,7 +461,7 @@ void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
IpplTimings::stopTimer(timeIntegrationTimer2_m);
}
-void ParallelTTracker::selectDT() {
+void ParallelTTracker::selectDT(bool backTrack) {
if (itsBunch_m->getIfBeamEmitting()) {
double dt = itsBunch_m->getEmissionDeltaT();
@@ -474,10 +470,13 @@ void ParallelTTracker::selectDT() {
double dt = dtCurrentTrack_m;
itsBunch_m->setdT(dt);
}
+ if (backTrack) {
+ itsBunch_m->setdT(-std::abs(itsBunch_m->getdT()));
+ }
}
-void ParallelTTracker::changeDT() {
- selectDT();
+void ParallelTTracker::changeDT(bool backTrack) {
+ selectDT(backTrack);
const unsigned int localNum = itsBunch_m->getLocalNum();
for (unsigned int i = 0; i < localNum; ++ i) {
itsBunch_m->dt[i] = itsBunch_m->getdT();
@@ -645,7 +644,7 @@ void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
if (ne > 0) {
msg << level1 << "* Deleted " << ne << " particles, "
- << "remaining " << itsBunch_m->getTotalNum() << " particles" << endl;
+ << "remaining " << totalParticlesInSimulation_m << " particles" << endl;
}
}
@@ -1054,7 +1053,7 @@ void ParallelTTracker::writePhaseSpace(const long long step, bool psDump, bool s
if (psDump && (itsBunch_m->getTotalNum() > 0)) {
// Write fields to .h5 file.
const size_t localNum = itsBunch_m->getLocalNum();
- double distToLastStop = zStop_m.back() - pathLength_m;
+ double distToLastStop = stepSizes_m.getFinalZStop() - pathLength_m;
Vector_t beta = itsBunch_m->RefPartP_m / Util::getGamma(itsBunch_m->RefPartP_m);
Vector_t driftPerTimeStep = itsBunch_m->getdT() * Physics::c * itsBunch_m->toLabTrafo_m.rotateFrom(beta);
bool driftToCorrectPosition = std::abs(distToLastStop) < 0.5 * euclidean_norm(driftPerTimeStep);
@@ -1080,7 +1079,7 @@ void ParallelTTracker::writePhaseSpace(const long long step, bool psDump, bool s
Quaternion(1.0, 0.0, 0.0, 0.0));
itsBunch_m->toLabTrafo_m = itsBunch_m->toLabTrafo_m * update.inverted();
- itsBunch_m->set_sPos(zStop_m.back());
+ itsBunch_m->set_sPos(stepSizes_m.getFinalZStop());
itsBunch_m->calcBeamParameters();
}
@@ -1111,7 +1110,9 @@ void ParallelTTracker::updateReference(const BorisPusher &pusher) {
}
void ParallelTTracker::updateReferenceParticle(const BorisPusher &pusher) {
- const double dt = std::min(itsBunch_m->getT(), itsBunch_m->getdT());
+ const double direction = back_track? -1: 1;
+ const double dt = direction * std::min(itsBunch_m->getT(),
+ direction * itsBunch_m->getdT());
const double scaleFactor = Physics::c * dt;
Vector_t Ef(0.0), Bf(0.0);
@@ -1162,10 +1163,9 @@ void ParallelTTracker::transformBunch(const CoordinateSystemTrafo &trafo) {
void ParallelTTracker::updateRefToLabCSTrafo() {
Vector_t R = itsBunch_m->toLabTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
-
Vector_t P = itsBunch_m->toLabTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
- pathLength_m += euclidean_norm(R);
+ pathLength_m += std::copysign(1, itsBunch_m->getdT()) * euclidean_norm(R);
CoordinateSystemTrafo update(R, getQuaternion(P, Vector_t(0, 0, 1)));
@@ -1174,15 +1174,40 @@ void ParallelTTracker::updateRefToLabCSTrafo() {
itsBunch_m->toLabTrafo_m = itsBunch_m->toLabTrafo_m * update.inverted();
}
+void ParallelTTracker::applyFractionalStep(const BorisPusher &pusher, double tau) {
+ double t = itsBunch_m->getT();
+ t += tau;
+ itsBunch_m->setT(t);
+
+ // the push method below pushes for half a time step. Hence the ref particle
+ // should be pushed for 2 * tau
+ itsBunch_m->RefPartR_m /= (Physics::c * 2 * tau);
+ pusher.push(itsBunch_m->RefPartR_m, itsBunch_m->RefPartP_m, tau);
+ itsBunch_m->RefPartR_m *= (Physics::c * 2 * tau);
+
+ pathLength_m = zstart_m;
+
+ Vector_t R = itsBunch_m->toLabTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
+ Vector_t P = itsBunch_m->toLabTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
+ CoordinateSystemTrafo update(R, getQuaternion(P, Vector_t(0, 0, 1)));
+ itsBunch_m->toLabTrafo_m = itsBunch_m->toLabTrafo_m * update.inverted();
+}
+
void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
+ StepSizeConfig stepSizesCopy(stepSizes_m);
+ if (back_track) {
+ stepSizesCopy.shiftZStopLeft(zstart_m);
+ }
+
double t = 0.0;
itsBunch_m->setT(t);
- dtCurrentTrack_m = dtAllTracks_m.front();
- changeDT();
+ dtCurrentTrack_m = stepSizesCopy.getdT();
+ selectDT();
- if (Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
+ if ((back_track && itsOpalBeamline_m.containsSource()) ||
+ Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
double gamma = 0.1 / itsBunch_m->getM() + 1.0;
double beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
itsBunch_m->RefPartP_m = itsBunch_m->toLabTrafo_m.rotateTo(beta * gamma * Vector_t(0, 0, 1));
@@ -1198,37 +1223,32 @@ void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
updateReferenceParticle(pusher);
pathLength_m += euclidean_norm(itsBunch_m->RefPartR_m - oldR);
- if (pathLength_m > zStop_m.front()) {
- if (localTrackSteps_m.size() == 0) return;
-
- dtAllTracks_m.pop();
- localTrackSteps_m.pop();
- zStop_m.pop();
-
- changeDT();
- }
-
double speed = euclidean_norm(itsBunch_m->RefPartP_m * Physics::c / Util::getGamma(itsBunch_m->RefPartP_m));
- if (std::abs(pathLength_m - zstart_m) <= 0.5 * itsBunch_m->getdT() * speed) {
- double tau = (pathLength_m - zstart_m) / speed;
- t += tau;
- itsBunch_m->setT(t);
+ if (pathLength_m > stepSizesCopy.getZStop()) {
+ ++ stepSizesCopy;
- itsBunch_m->RefPartR_m /= (Physics::c * tau);
- pusher.push(itsBunch_m->RefPartR_m, itsBunch_m->RefPartP_m, tau);
- itsBunch_m->RefPartR_m *= (Physics::c * tau);
+ if (stepSizesCopy.reachedEnd()) {
+ -- stepSizesCopy;
+ double tau = (stepSizesCopy.getZStop() - pathLength_m) / speed;
+ applyFractionalStep(pusher, tau);
- pathLength_m = zstart_m;
+ break;
+ }
- Vector_t R = itsBunch_m->toLabTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
- Vector_t P = itsBunch_m->toLabTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
- CoordinateSystemTrafo update(R, getQuaternion(P, Vector_t(0, 0, 1)));
- itsBunch_m->toLabTrafo_m = itsBunch_m->toLabTrafo_m * update.inverted();
+ dtCurrentTrack_m = stepSizesCopy.getdT();
+ selectDT();
+ }
- return;
+ if (std::abs(pathLength_m - zstart_m) <= 0.5 * itsBunch_m->getdT() * speed) {
+ double tau = (zstart_m - pathLength_m) / speed;
+ applyFractionalStep(pusher, tau);
+
+ break;
}
}
+
+ changeDT();
}
void ParallelTTracker::autophaseCavities(const BorisPusher &pusher) {
diff --git a/src/Algorithms/ParallelTTracker.h b/src/Algorithms/ParallelTTracker.h
index dc35e2c715accf70afbbc7d4d5e179738b47ff04..c495cf0b98b2749f8c884d39ab87c75aada0c2fd 100644
--- a/src/Algorithms/ParallelTTracker.h
+++ b/src/Algorithms/ParallelTTracker.h
@@ -21,6 +21,7 @@
#include "Algorithms/Tracker.h"
#include "Steppers/BorisPusher.h"
#include "Structure/DataSink.h"
+#include "Algorithms/StepSizeConfig.h"
#include "BasicActions/Option.h"
#include "Utilities/Options.h"
@@ -62,9 +63,7 @@
#include <list>
#include <vector>
-#include <queue>
-class BorisPusher;
class ParticleMatterInteractionHandler;
class ParallelTTracker: public Tracker {
@@ -214,14 +213,13 @@ private:
/// where to start
double zstart_m;
- /// where to stop
- std::queue<double> zStop_m;
+ /// stores informations where to change the time step and
+ /// where to stop the simulation,
+ /// the time step sizes and
+ /// the number of time steps with each configuration
+ StepSizeConfig stepSizes_m;
double dtCurrentTrack_m;
- std::queue<double> dtAllTracks_m;
-
- /// The maximal number of steps the system is integrated per TRACK
- std::queue<unsigned long long> localTrackSteps_m;
// This variable controls the minimal number of steps of emission (using bins)
// before we can merge the bins
@@ -292,8 +290,8 @@ private:
void timeIntegration1(BorisPusher & pusher);
void timeIntegration2(BorisPusher & pusher);
- void selectDT();
- void changeDT();
+ void selectDT(bool backTrack = false);
+ void changeDT(bool backTrack = false);
void emitParticles(long long step);
void computeExternalFields(OrbitThreader &oth);
void computeWakefield(IndexMap::value_t &elements);
@@ -312,12 +310,11 @@ private:
void updateReference(const BorisPusher &pusher);
void updateRefToLabCSTrafo();
+ void applyFractionalStep(const BorisPusher &pusher, double tau);
void findStartPosition(const BorisPusher &pusher);
void autophaseCavities(const BorisPusher &pusher);
void evenlyDistributeParticles();
-
- static unsigned long long getMaxSteps(std::queue<unsigned long long> numSteps);
};
inline void ParallelTTracker::visitAlignWrapper(const AlignWrapper &wrap) {
@@ -455,18 +452,6 @@ inline void ParallelTTracker::pushParticles(const BorisPusher &pusher) {
itsBunch_m->switchOffUnitlessPositions(true);
}
-inline
-unsigned long long ParallelTTracker::getMaxSteps(std::queue<unsigned long long> numSteps) {
- unsigned long long totalNumSteps = 0;
-
- while (numSteps.size() > 0) {
- totalNumSteps += numSteps.front();
- numSteps.pop();
- }
-
- return totalNumSteps;
-}
-
#ifdef OPAL_DKS
inline void ParallelTTracker::setupDKS() {
dksbase = new DKSOPAL;
diff --git a/src/Algorithms/StepSizeConfig.cpp b/src/Algorithms/StepSizeConfig.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f4808a8bf7eac9e7e6dd0892fe141165d2fc57a4
--- /dev/null
+++ b/src/Algorithms/StepSizeConfig.cpp
@@ -0,0 +1,181 @@
+#include "Algorithms/StepSizeConfig.h"
+#include "Utilities/OpalException.h"
+
+#include <algorithm>
+#include <numeric>
+#include <iterator>
+#include <cmath>
+
+void StepSizeConfig::sortAscendingZStop() {
+ configurations_m.sort([] (const entry_t &a,
+ const entry_t &b) -> bool
+ {
+ return std::get<1>(a) < std::get<1>(b);
+ });
+}
+
+void StepSizeConfig::reverseDirection() {
+ unsigned int posIterator = std::distance(it_m, configurations_m.end()) - 1;
+ configurations_m.reverse();
+ it_m = configurations_m.begin();
+ std::advance(it_m, posIterator);
+}
+
+StepSizeConfig& StepSizeConfig::advanceToPos(double spos) {
+ while (getZStop() < spos && std::next(it_m) != configurations_m.end()) {
+ ++ it_m;
+ }
+
+ return *this;
+}
+
+StepSizeConfig& StepSizeConfig::operator++() {
+ if (reachedEnd()) {
+ throw OpalException("StepSizeConfig::operator++",
+ "iterator is at end of list of configurations");
+ }
+
+ ++ it_m;
+
+ return *this;
+}
+
+StepSizeConfig& StepSizeConfig::operator--() {
+ if (reachedStart()) {
+ throw OpalException("StepSizeConfig::operator--",
+ "iterator is at begin of list of configurations");
+ }
+
+ -- it_m;
+
+ return *this;
+}
+
+void StepSizeConfig::shiftZStopRight(double front) {
+ auto it = configurations_m.begin();
+ while (std::get<1>(*it) < front &&
+ std::next(it) != configurations_m.end()) {
+ ++ it;
+ }
+
+ double zstop = std::get<1>(*it);
+ if (zstop < front) return;
+
+ std::get<1>(*it) = front;
+ for (++ it ;
+ it != configurations_m.end();
+ ++ it) {
+ std::swap(zstop, std::get<1>(*it));
+ }
+}
+
+void StepSizeConfig::shiftZStopLeft(double back) {
+ auto it = configurations_m.rbegin();
+ while (std::get<1>(*it) > back &&
+ std::next(it) != configurations_m.rend()) {
+ ++ it;
+ }
+
+ double zstop = std::get<1>(*it);
+ if (zstop > back) return;
+
+ std::get<1>(*it) = back;
+ for (++ it ;
+ it != configurations_m.rend();
+ ++ it) {
+ std::swap(zstop, std::get<1>(*it));
+ }
+}
+
+double StepSizeConfig::getdT() const {
+ if (reachedEnd()) {
+ throw OpalException("StepSizeConfig::getdT",
+ "iterator is at end of list of configurations");
+ }
+
+ return std::get<0>(*it_m);
+}
+
+double StepSizeConfig::getZStop() const {
+ if (reachedEnd()) {
+ throw OpalException("StepSizeConfig::getZStop",
+ "iterator is at end of list of configurations");
+ }
+
+ return std::get<1>(*it_m);
+}
+
+unsigned long StepSizeConfig::getNumSteps() const {
+ if (reachedEnd()) {
+ throw OpalException("StepSizeConfig::getNumSteps",
+ "iterator is at end of list of configurations");
+ }
+
+ return std::get<2>(*it_m);
+}
+
+unsigned long long StepSizeConfig::getMaxSteps() const {
+ unsigned long long maxSteps = 0;
+ for (const auto config: configurations_m) {
+ maxSteps += std::get<2>(config);
+ }
+
+ return maxSteps;
+}
+
+unsigned long long StepSizeConfig::getNumStepsFinestResolution() const {
+ double minTimeStep = std::get<0>(configurations_m.front());
+ unsigned long long totalNumSteps = 0;
+
+ for (const auto config: configurations_m) {
+ const double &dt = std::get<0>(config);
+ const unsigned long &numSteps = std::get<2>(config);
+
+ if (minTimeStep > dt) {
+ totalNumSteps = std::ceil(totalNumSteps * minTimeStep / dt);
+ minTimeStep = dt;
+ }
+
+ totalNumSteps += std::ceil(numSteps * dt / minTimeStep);
+ }
+
+ return totalNumSteps;
+}
+
+double StepSizeConfig::getMinTimeStep() const {
+ double minTimeStep = std::get<0>(configurations_m.front());
+ for (const auto config: configurations_m) {
+ if (minTimeStep > std::get<0>(config)) {
+ minTimeStep = std::get<0>(config);
+ }
+ }
+
+ return minTimeStep;
+}
+
+double StepSizeConfig::getFinalZStop() const {
+ return std::get<1>(configurations_m.back());
+}
+
+void StepSizeConfig::print(Inform &out) const {
+ out << std::scientific << " "
+ << std::setw(20) << "dt [ns] "
+ << std::setw(20) << "zStop [m] "
+ << std::setw(20) << "num Steps [1]"
+ << endl;
+
+ for (auto it = configurations_m.begin();
+ it != configurations_m.end();
+ ++ it) {
+ if (it_m == it) {
+ out << "-> ";
+ } else {
+ out << " ";
+ }
+
+ out << std::setw(20) << std::get<0>(*it)
+ << std::setw(20) << std::get<1>(*it)
+ << std::setw(20) << std::get<2>(*it)
+ << endl;
+ }
+}
\ No newline at end of file
diff --git a/src/Algorithms/StepSizeConfig.h b/src/Algorithms/StepSizeConfig.h
new file mode 100644
index 0000000000000000000000000000000000000000..2ece19c7a3e2aacc8aaf9d96ff491cbf69dabe49
--- /dev/null
+++ b/src/Algorithms/StepSizeConfig.h
@@ -0,0 +1,106 @@
+#ifndef STEPSIZECONFIG_H
+#define STEPSIZECONFIG_H
+
+#include "Utility/Inform.h"
+
+#include <list>
+#include <tuple>
+
+class StepSizeConfig {
+public:
+ StepSizeConfig();
+
+ StepSizeConfig(const StepSizeConfig &right);
+
+ void operator=(const StepSizeConfig &) = delete;
+
+ void push_back(double dt,
+ double zstop,
+ unsigned long numSteps);
+
+ void sortAscendingZStop();
+
+ void resetIterator();
+
+ bool reachedStart() const;
+
+ bool reachedEnd() const;
+
+ void clear();
+
+ void reverseDirection();
+
+ StepSizeConfig& advanceToPos(double spos);
+
+ StepSizeConfig& operator++();
+
+ StepSizeConfig& operator--();
+
+ void shiftZStopRight(double front);
+ void shiftZStopLeft(double back);
+
+ double getdT() const;
+
+ double getZStop() const;
+
+ unsigned long getNumSteps() const;
+
+ unsigned long long getMaxSteps() const;
+
+ unsigned long long getNumStepsFinestResolution() const;
+
+ double getMinTimeStep() const;
+
+ double getFinalZStop() const;
+
+ void print(Inform &out) const;
+
+private:
+ typedef std::tuple<double, double, unsigned long> entry_t;
+ typedef std::list<entry_t> container_t;
+
+ container_t configurations_m;
+ container_t::iterator it_m;
+};
+
+inline
+StepSizeConfig::StepSizeConfig():
+ configurations_m(),
+ it_m(configurations_m.begin())
+{ }
+
+inline
+StepSizeConfig::StepSizeConfig(const StepSizeConfig &right):
+ configurations_m(right.configurations_m),
+ it_m(configurations_m.begin())
+{ }
+
+inline
+void StepSizeConfig::push_back(double dt,
+ double zstop,
+ unsigned long numSteps) {
+ configurations_m.push_back(std::make_tuple(dt, zstop, numSteps));
+}
+
+inline
+void StepSizeConfig::resetIterator() {
+ it_m = configurations_m.begin();
+}
+
+inline
+bool StepSizeConfig::reachedStart() const {
+ return (it_m == configurations_m.begin());
+}
+
+inline
+bool StepSizeConfig::reachedEnd() const {
+ return (it_m == configurations_m.end());
+}
+
+inline
+void StepSizeConfig::clear() {
+ configurations_m.clear();
+ it_m = configurations_m.begin();
+}
+
+#endif
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/Monitor.cpp b/src/Classic/AbsBeamline/Monitor.cpp
index a460fc7c170c5a40ffb5cc3c453e6575a75e0303..78a8c454b71a2abf2356df6065241813335c5f3f 100644
--- a/src/Classic/AbsBeamline/Monitor.cpp
+++ b/src/Classic/AbsBeamline/Monitor.cpp
@@ -78,7 +78,8 @@ bool Monitor::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B)
const double recpgamma = Physics::c * dt / Util::getGamma(P);
const double middle = 0.5 * getElementLength();
if (online_m && type_m == SPATIAL) {
- if (R(2) < middle && R(2) + P(2) * recpgamma > middle) {
+ if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking
+ dt * (R(2) + P(2) * recpgamma) > dt * middle) {
double frac = (middle - R(2)) / (P(2) * recpgamma);
lossDs_m->addParticle(R + frac * recpgamma * P,
@@ -99,7 +100,8 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R,
const double recpgamma = Physics::c * dt / Util::getGamma(P);
const double middle = 0.5 * getElementLength();
- if (R(2) < middle && R(2) + P(2) * recpgamma > middle) {
+ if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking
+ dt * (R(2) + P(2) * recpgamma) > dt * middle) {
double frac = (middle - R(2)) / (P(2) * recpgamma);
double time = t + frac * dt;
Vector_t dR = (0.5 + frac) * P * recpgamma;
@@ -115,7 +117,7 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R,
for (unsigned int i = 0; i < localNum; ++ i) {
const double recpgamma = Physics::c * dt / Util::getGamma(RefPartBunch_m->P[i]);
- lossDs_m->addParticle(RefPartBunch_m->R[i] + frac * RefPartBunch_m->P[i] * recpgamma,
+ lossDs_m->addParticle(RefPartBunch_m->R[i] + frac * RefPartBunch_m->P[i] * recpgamma - halfLength_s,
RefPartBunch_m->P[i], RefPartBunch_m->ID[i],
time, 0);
}
diff --git a/src/Classic/AbsBeamline/Source.cpp b/src/Classic/AbsBeamline/Source.cpp
index 1a6cb856d26cec3497f4f0e583d2c9c8f7cce5cf..c96f309137ffe7cd5b5e14a9cb0248ade7bc49af 100644
--- a/src/Classic/AbsBeamline/Source.cpp
+++ b/src/Classic/AbsBeamline/Source.cpp
@@ -14,8 +14,6 @@ extern Inform *gmsg;
// Class Source
// ------------------------------------------------------------------------
-const std::string Source::defaultDistribution("DISTRIBUTION");
-
Source::Source():
Source("")
{}
@@ -65,26 +63,37 @@ void Source::addKR(int i, double t, Vector_t &K) {
void Source::addKT(int i, double t, Vector_t &K) {
}
-bool Source::apply(const double &t) {
- Vector_t externalE = Vector_t(0.0);
- Vector_t externalB = Vector_t(0.0);
- beamline_m->getFieldAt(csTrafoGlobal2Local_m.getOrigin(),
- Vector_t(0.0),
- RefPartBunch_m->getT() + 0.5 * RefPartBunch_m->getdT(),
- externalE,
- externalB);
- for (Distribution* dist: distrs_m) {
- dist->emitParticles(RefPartBunch_m, externalE(2));
+bool Source::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
+ const Vector_t &R = RefPartBunch_m->R[i];
+ const Vector_t &P = RefPartBunch_m->P[i];
+ const double &dt = RefPartBunch_m->dt[i];
+ const double recpgamma = Physics::c * dt / Util::getGamma(P);
+ const double end = getElementLength();
+ if (online_m && dt < 0.0) {
+ if (R(2) > end &&
+ (R(2) + P(2) * recpgamma) < end) {
+ double frac = (end - R(2)) / (P(2) * recpgamma);
+
+ lossDs_m->addParticle(R + frac * recpgamma * P,
+ P, RefPartBunch_m->ID[i], t + frac * dt, 0);
+
+ return true;
+ }
}
return false;
}
void Source::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
- // Inform msg("Source ", *gmsg);
- // double zBegin = 0.0, zEnd = 0.0, rBegin = 0.0, rEnd = 0.0;
-
RefPartBunch_m = bunch;
+ endField = startField;
+ startField -= getElementLength();
+
+ std::string filename = getName();
+ lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename,
+ !Options::asciidump,
+ getType()));
+
}
void Source::finalise()
@@ -101,6 +110,7 @@ void Source::goOnline(const double &) {
void Source::goOffline() {
online_m = false;
+ lossDs_m->save();
}
void Source::getDimensions(double &zBegin, double &zEnd) const {
@@ -111,16 +121,4 @@ void Source::getDimensions(double &zBegin, double &zEnd) const {
ElementBase::ElementType Source::getType() const {
return SOURCE;
-}
-
-void Source::setDistribution(std::vector<std::string> distNames) {
- const size_t numberOfDistributions = distNames.size();
- if (numberOfDistributions == 0) {
- distrs_m.push_back(Distribution::find(defaultDistribution));
- } else {
- for (std::string name: distNames) {
- distrs_m.push_back(Distribution::find(name));
- distrs_m.back()->setNumberOfDistributions(numberOfDistributions);
- }
- }
}
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/Source.h b/src/Classic/AbsBeamline/Source.h
index e62ea0c329a1caf798c1bb5d2a9ff69600fded52..830eb96546abc174bdab970b0823e7406ce93829 100644
--- a/src/Classic/AbsBeamline/Source.h
+++ b/src/Classic/AbsBeamline/Source.h
@@ -2,7 +2,8 @@
#define CLASSIC_SOURCE_HH
#include "AbsBeamline/Component.h"
-#include "Distribution/Distribution.h"
+#include "Structure/LossDataSink.h"
+
class OpalBeamline;
template <class T, unsigned Dim>
@@ -26,8 +27,7 @@ public:
virtual void addKT(int i, double t, Vector_t &K) override;
- using Component::apply;
- virtual bool apply(const double &t);
+ virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) override;
virtual void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) override;
@@ -43,26 +43,15 @@ public:
virtual void getDimensions(double &zBegin, double &zEnd) const override;
- void setDistribution(std::vector<std::string> distNames);
-
- void setBeamline(OpalBeamline *beamline);
private:
double ElementEdge_m;
double startField_m; /**< startingpoint of field, m*/
double endField_m;
- std::vector<Distribution*> distrs_m;
+ std::unique_ptr<LossDataSink> lossDs_m;
- OpalBeamline *beamline_m;
// Not implemented.
void operator=(const Source &);
-
- static const std::string defaultDistribution;
};
-
-inline
-void Source::setBeamline(OpalBeamline *beamline) {
- beamline_m = beamline;
-}
-#endif // CLASSIC_SOURCE_HH
+#endif // CLASSIC_SOURCE_HH
\ No newline at end of file
diff --git a/src/Classic/Algorithms/PartBunchBase.hpp b/src/Classic/Algorithms/PartBunchBase.hpp
index 3df9b5d0bc5beecc0234897ff3928ada7a0cbece..2c1d57cd69fdcd56772142ee96b5320f51de051b 100644
--- a/src/Classic/Algorithms/PartBunchBase.hpp
+++ b/src/Classic/Algorithms/PartBunchBase.hpp
@@ -1191,7 +1191,7 @@ Vector_t PartBunchBase<T, Dim>::get_pmean() const {
template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_pmean_Distribution() const {
- if (dist_m && dist_m->getType() != DistrTypeT::FROMFILE)
+ if (dist_m)// && dist_m->getType() != DistrTypeT::FROMFILE)
return dist_m->get_pmean();
double gamma = 0.1 / getM() + 1; // set default 0.1 eV
diff --git a/src/Classic/Structure/LossDataSink.cpp b/src/Classic/Structure/LossDataSink.cpp
index 4238a77a977620cfd4bf2944e91d7713dd2ef046..e549fb3ee14fe56bc70f15fa0dfc47aab0bc8d1f 100644
--- a/src/Classic/Structure/LossDataSink.cpp
+++ b/src/Classic/Structure/LossDataSink.cpp
@@ -103,8 +103,14 @@
SetStatistics::SetStatistics():
+ element_m(""),
+ spos_m(0.0),
+ refTime_m(0.0),
tmean_m(0.0),
trms_m(0.0),
+ nTotal_m(0),
+ RefPartR_m(0.0),
+ RefPartP_m(0.0),
rmin_m(0.0),
rmax_m(0.0),
rmean_m(0.0),
@@ -627,6 +633,7 @@ SetStatistics LossDataSink::computeSetStatistics(unsigned int setIdx) {
new_reduce(rminmax, 6, std::greater<double>());
if (Ippl::myNode() != 0) return stat;
+ if (plainData[0] == 0.0) return stat;
double *centroid = plainData + 1;
double *moments = plainData + 7;
diff --git a/src/Classic/Structure/LossDataSink.h b/src/Classic/Structure/LossDataSink.h
index 8f5113130c72a3b205083babec69ed3f86d57e9c..900cb6723e27d28d94679a1d7b98f8d6df966377 100644
--- a/src/Classic/Structure/LossDataSink.h
+++ b/src/Classic/Structure/LossDataSink.h
@@ -182,7 +182,10 @@ std::set<SetStatistics> LossDataSink::computeStatistics(unsigned int numStatisti
splitSets(numStatistics);
for (unsigned int i = 0; i < numStatistics; ++ i) {
- stats.insert(computeSetStatistics(i));
+ auto setStats = computeSetStatistics(i);
+ if (setStats.nTotal_m > 0) {
+ stats.insert(setStats);
+ }
}
return stats;
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index efee18c1018fff4889f734db0e61ea2e3d5f2ac2..322b8772bbc38a2287f19875fedcb5b5e5fa0bbe 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -1328,10 +1328,10 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
double rguess =
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]);
-
+
double denergy = 1000.0 *
Attributes::getReal(itsAttr[Attrib::Distribution::DENERGY]);
-
+
if ( denergy < 0.0 )
throw OpalException("Distribution:CreateMatchedGaussDistribution()",
"DENERGY < 0");
@@ -4401,7 +4401,7 @@ void Distribution::shiftBeam(double &maxTOrZ, double &minTOrZ) {
maxTOrZ -= maxTOrZ;
}
- } else {
+ } else if (distrTypeT_m != DistrTypeT::FROMFILE) {
double avgZ[] = {0.0, 1.0 * tOrZDist_m.size()};
for (double tOrZ : tOrZDist_m)
avgZ[0] += tOrZ;
diff --git a/src/Elements/OpalBeamline.cpp b/src/Elements/OpalBeamline.cpp
index 24f3c4c4206d033850740c87aefddb7a33207e29..a852ce7379771a406ecec9330aaaa07d39999a52 100644
--- a/src/Elements/OpalBeamline.cpp
+++ b/src/Elements/OpalBeamline.cpp
@@ -214,58 +214,60 @@ void OpalBeamline::compute3DLattice() {
}
(*it).order_m = minOrder;
std::shared_ptr<Component> element = (*it).getElement();
- if (element->getType() == ElementBase::SBEND ||
- element->getType() == ElementBase::RBEND ||
- element->getType() == ElementBase::RBEND3D) {
-
- double beginThisPathLength = element->getElementPosition();
- Vector_t beginThis3D(0, 0, beginThisPathLength - endPriorPathLength);
- BendBase * bendElement = static_cast<BendBase*>(element.get());
- double thisLength = bendElement->getChordLength();
- double bendAngle = bendElement->getBendAngle();
- double entranceAngle = bendElement->getEntranceAngle();
- double arcLength = (thisLength * std::abs(bendAngle) / (2 * sin(std::abs(bendAngle) / 2)));
-
- double rotationAngleAboutZ = bendElement->getRotationAboutZ();
- Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ),
- sin(-0.5 * rotationAngleAboutZ) * Vector_t(0, 0, 1));
-
- Vector_t effectiveRotationAxis = rotationAboutZ.rotate(Vector_t(0, -1, 0));
- effectiveRotationAxis /= euclidean_norm(effectiveRotationAxis);
-
- Quaternion_t rotationAboutAxis(cos(0.5 * bendAngle),
- sin(0.5 * bendAngle) * effectiveRotationAxis);
- Quaternion_t halfRotationAboutAxis(cos(0.25 * bendAngle),
- sin(0.25 * bendAngle) * effectiveRotationAxis);
- Quaternion_t entryFaceRotation(cos(0.5 * entranceAngle),
- sin(0.5 * entranceAngle) * effectiveRotationAxis);
-
- if (!Options::idealized) {
- std::vector<Vector_t> truePath = bendElement->getDesignPath();
- Quaternion_t directionExitHardEdge(cos(0.5 * (0.5 * bendAngle - entranceAngle)),
- sin(0.5 * (0.5 * bendAngle - entranceAngle)) * effectiveRotationAxis);
- Vector_t exitHardEdge = thisLength * directionExitHardEdge.rotate(Vector_t(0, 0, 1));
- double distanceEntryHETruePath = euclidean_norm(truePath.front());
- double distanceExitHETruePath = euclidean_norm(rotationAboutZ.rotate(truePath.back()) - exitHardEdge);
- double pathLengthTruePath = (*it).getEnd() - (*it).getStart();
- arcLength = pathLengthTruePath - distanceEntryHETruePath - distanceExitHETruePath;
- }
- Vector_t chord = thisLength * halfRotationAboutAxis.rotate(Vector_t(0, 0, 1));
- Vector_t endThis3D = beginThis3D + chord;
- double endThisPathLength = beginThisPathLength + arcLength;
+ if (element->getType() != ElementBase::SBEND &&
+ element->getType() != ElementBase::RBEND &&
+ element->getType() != ElementBase::RBEND3D) {
+ continue;
+ }
+
+ double beginThisPathLength = element->getElementPosition();
+ Vector_t beginThis3D(0, 0, beginThisPathLength - endPriorPathLength);
+ BendBase * bendElement = static_cast<BendBase*>(element.get());
+ double thisLength = bendElement->getChordLength();
+ double bendAngle = bendElement->getBendAngle();
+ double entranceAngle = bendElement->getEntranceAngle();
+ double arcLength = (thisLength * std::abs(bendAngle) / (2 * sin(std::abs(bendAngle) / 2)));
- CoordinateSystemTrafo fromEndLastToBeginThis(beginThis3D,
- (entryFaceRotation * rotationAboutZ).conjugate());
- CoordinateSystemTrafo fromEndLastToEndThis(endThis3D,
- rotationAboutAxis.conjugate());
+ double rotationAngleAboutZ = bendElement->getRotationAboutZ();
+ Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ),
+ sin(-0.5 * rotationAngleAboutZ) * Vector_t(0, 0, 1));
- (*it).setCoordTransformationTo(fromEndLastToBeginThis * currentCoordTrafo);
+ Vector_t effectiveRotationAxis = rotationAboutZ.rotate(Vector_t(0, -1, 0));
+ effectiveRotationAxis /= euclidean_norm(effectiveRotationAxis);
- currentCoordTrafo = (fromEndLastToEndThis * currentCoordTrafo);
+ Quaternion_t rotationAboutAxis(cos(0.5 * bendAngle),
+ sin(0.5 * bendAngle) * effectiveRotationAxis);
+ Quaternion_t halfRotationAboutAxis(cos(0.25 * bendAngle),
+ sin(0.25 * bendAngle) * effectiveRotationAxis);
+ Quaternion_t entryFaceRotation(cos(0.5 * entranceAngle),
+ sin(0.5 * entranceAngle) * effectiveRotationAxis);
- endPriorPathLength = endThisPathLength;
+ if (!Options::idealized) {
+ std::vector<Vector_t> truePath = bendElement->getDesignPath();
+ Quaternion_t directionExitHardEdge(cos(0.5 * (0.5 * bendAngle - entranceAngle)),
+ sin(0.5 * (0.5 * bendAngle - entranceAngle)) * effectiveRotationAxis);
+ Vector_t exitHardEdge = thisLength * directionExitHardEdge.rotate(Vector_t(0, 0, 1));
+ double distanceEntryHETruePath = euclidean_norm(truePath.front());
+ double distanceExitHETruePath = euclidean_norm(rotationAboutZ.rotate(truePath.back()) - exitHardEdge);
+ double pathLengthTruePath = (*it).getEnd() - (*it).getStart();
+ arcLength = pathLengthTruePath - distanceEntryHETruePath - distanceExitHETruePath;
}
+
+ Vector_t chord = thisLength * halfRotationAboutAxis.rotate(Vector_t(0, 0, 1));
+ Vector_t endThis3D = beginThis3D + chord;
+ double endThisPathLength = beginThisPathLength + arcLength;
+
+ CoordinateSystemTrafo fromEndLastToBeginThis(beginThis3D,
+ (entryFaceRotation * rotationAboutZ).conjugate());
+ CoordinateSystemTrafo fromEndLastToEndThis(endThis3D,
+ rotationAboutAxis.conjugate());
+
+ (*it).setCoordTransformationTo(fromEndLastToBeginThis * currentCoordTrafo);
+
+ currentCoordTrafo = (fromEndLastToEndThis * currentCoordTrafo);
+
+ endPriorPathLength = endThisPathLength;
}
}
@@ -283,6 +285,9 @@ void OpalBeamline::compute3DLattice() {
double thisLength = element->getElementLength();
Vector_t beginThis3D(0, 0, beginThisPathLength - endPriorPathLength);
+ if (element->getType() == ElementBase::SOURCE) {
+ beginThis3D(2) -= thisLength;
+ }
if (element->getType() == ElementBase::MONITOR) {
beginThis3D(2) -= 0.5 * thisLength;
}
@@ -587,4 +592,4 @@ void OpalBeamline::activateElements() {
(*it).setOn(designEnergy);
// element->goOnline(designEnergy);
}
-}
\ No newline at end of file
+}
diff --git a/src/Elements/OpalBeamline.h b/src/Elements/OpalBeamline.h
index 6ef4455d236edb5b56c76f5480af820b39cd47e9..9eced4f53b60d1823f470ddaed92da84f76e3f2e 100644
--- a/src/Elements/OpalBeamline.h
+++ b/src/Elements/OpalBeamline.h
@@ -111,6 +111,15 @@ void OpalBeamline::visit(const T &element, BeamlineVisitor &, PartBunchBase<doub
template<> inline
void OpalBeamline::visit<Source>(const Source &element, BeamlineVisitor &, PartBunchBase<double, 3> *bunch) {
containsSource_m = true;
+ double startField = 0.0;
+ double endField = 0.0;
+ std::shared_ptr<Source> elptr(dynamic_cast<Source *>(element.removeWrappers()->clone()));
+
+ if (elptr->isElementPositionSet())
+ startField = elptr->getElementPosition();
+
+ elptr->initialise(bunch, startField, endField);
+ elements_m.push_back(ClassicField(elptr, startField, endField));
}
template<> inline
diff --git a/src/Elements/OpalSource.cpp b/src/Elements/OpalSource.cpp
index 330f44d22d661535c8960af47718aafb94ecd6f5..0a912825f8d3d712bf3a3f418c0b01e195fba352 100644
--- a/src/Elements/OpalSource.cpp
+++ b/src/Elements/OpalSource.cpp
@@ -49,13 +49,10 @@ void OpalSource::update() {
SourceRep *sol =
dynamic_cast<SourceRep *>(getElement()->removeWrappers());
- double length = Attributes::getReal(itsAttr[LENGTH]);
- std::vector<std::string> distribution = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
+ double length = 0.05;
sol->setElementLength(length);
- sol->setDistribution(distribution);
-
// Transmit "unknown" attributes.
OpalElement::updateUnknown(sol);
}
\ No newline at end of file
diff --git a/src/Track/TrackRun.cpp b/src/Track/TrackRun.cpp
index 40822dd81f03949678dff62a6dd4e50e4c5657f1..4eafbebd935ec03c19aa38ee314d488f69885f56 100644
--- a/src/Track/TrackRun.cpp
+++ b/src/Track/TrackRun.cpp
@@ -81,6 +81,7 @@ namespace {
BOUNDARYGEOMETRY, // The boundary geometry
DISTRIBUTION, // The particle distribution
MULTIPACTING, // MULTIPACTING flag
+ TRACKBACK,
SIZE
};
}
@@ -125,7 +126,8 @@ TrackRun::TrackRun():
("DISTRIBUTION", "List of particle distributions to be used ");
itsAttr[MULTIPACTING] = Attributes::makeBool
("MULTIPACTING", "Multipacting flag, default: false. Set true to initialize primary particles according to BoundaryGeometry", false);
-
+ itsAttr[TRACKBACK] = Attributes::makeBool
+ ("TRACKBACK", "Track in reverse direction, default: false", false);
registerOwnership(AttributeHandler::SUB_COMMAND);
opal = OpalData::getInstance();
@@ -318,19 +320,15 @@ void TrackRun::setupSliceTracker() {
if (Track::block->bunch->getTotalNum() > 0) {
double spos = Track::block->zstart;
auto &zstop = Track::block->zstop;
- auto &timeStep = Track::block->localTimeSteps;
- auto &dT = Track::block->dT;
+ auto it = Track::block->dT.begin();
unsigned int i = 0;
while (i + 1 < zstop.size() && zstop[i + 1] < spos) {
++ i;
+ ++ it;
}
- zstop.erase(zstop.begin(), zstop.begin() + i);
- timeStep.erase(timeStep.begin(), timeStep.begin() + i);
- dT.erase(dT.begin(), dT.begin() + i);
-
- Track::block->bunch->setdT(dT.front());
+ Track::block->bunch->setdT(*it);
} else {
Track::block->zstart = 0.0;
}
@@ -550,19 +548,15 @@ void TrackRun::setupTTracker(){
if (Track::block->bunch->getTotalNum() > 0) {
double spos = Track::block->zstart;
auto &zstop = Track::block->zstop;
- auto &timeStep = Track::block->localTimeSteps;
- auto &dT = Track::block->dT;
+ auto it = Track::block->dT.begin();
unsigned int i = 0;
while (i + 1 < zstop.size() && zstop[i + 1] < spos) {
++ i;
+ ++ it;
}
- zstop.erase(zstop.begin(), zstop.begin() + i);
- timeStep.erase(timeStep.begin(), timeStep.begin() + i);
- dT.erase(dT.begin(), dT.begin() + i);
-
- Track::block->bunch->setdT(dT.front());
+ Track::block->bunch->setdT(*it);
} else {
Track::block->zstart = 0.0;
}
@@ -585,7 +579,7 @@ void TrackRun::setupTTracker(){
*ds,
Track::block->reference,
false,
- false,
+ Attributes::getBool(itsAttr[TRACKBACK]),
Track::block->localTimeSteps,
Track::block->zstart,
Track::block->zstop,
diff --git a/tools/emacs/opal.el b/tools/emacs/opal.el
index c2efbfaf467978c7d16ce6514e31da71bfb99466..9bbb23bea007d3d04c7585d03364f3cb382b4100 100644
--- a/tools/emacs/opal.el
+++ b/tools/emacs/opal.el
@@ -138,11 +138,11 @@
)
"Highlighting expressions for OPAL mode (seqediting).")
-;(concat "\\<" (regexp-opt '("A" "ADD" "AMPLITUDE_MODEL" "ANGLE" "APERTURE" "APVETO" "AT" "AUTOPHASE" "B" "BBOXINCR" "BCFFTZ" "BCFFTT" "BCFFTX" "BCFFTY" "BCURRENT" "BEAM_PHIINIT" "BEAM_PRINIT" "BEAM_RINIT" "BFREQ" "BIRTH_CONTROL" "BMAX" "BOUNDPDESTROYFQ" "BSCALE" "BY" "CALLS" "CATHTEMP" "CENTRE" "CHARGE" "CLASS" "CLEAR" "CMD" "COEFDENOM" "COEFDENOMPHI" "COEFNUM" "COEFNUMPHI" "COLUMN" "CONDUCT" "CONSTRAINTS" "CONST_LENGTH" "CONV_HVOL_PROG" "CORRX" "CORRY" "CORRZ" "CSRDUMP" "CROSSOVER" "CUTOFFLONG" "CUTOFFPX" "CUTOFFPY" "CUTOFFPZ" "CUTOFFR" "CUTOFFX" "CUTOFFY" "CYHARMON" "CZERO" "DESIGNENERGY" "DK1" "DK1S" "DK2" "DK2S" "DKN" "DKNR" "DKS" "DKSR" "DLAG" "DPHI" "DPSI" "DS" "DT" "DTHETA" "DUMP" "DUMP_DAT" "DUMP_FREQ" "DUMP_OFFSPRING" "DVARS" "DVOLT" "DX" "DY" "DZ" "E1" "E2" "EBDUMP" "ECHO" "EKIN" "ELASER" "ELEMEDGE" "ELEMENT" "EMISSIONMODEL" "EMISSIONSTEPS" "EMITTED" "ENABLEHDF5" "ENABLERUTHERFORD" "ENERGY" "ENDSEQUENCE" "END_NORMAL_X" "END_NORMAL_Y" "END_POSITION_X" "END_POSITION_Y" "EPSILON" "ESCALE" "ET" "EVERYSTEP" "EX" "EXPECTED_HYPERVOL" "EXPR" "EY" "FE" "FGEOM" "FIELDMAPDIR" "FIELDMAPDIR" "FILE" "FINT" "FMAPFN" "FNAME" "FORM" "FREQ" "FREQUENCY_MODEL" "FROM" "FSTYPE" "FTOSCAMPLITUDE" "FTOSCPERIODS" "FULL" "GAMMA" "GAP" "GAPWIDTH" "GENE_MUTATION_PROBABILITY" "GEOMETRY" "GREENSF" "H1" "H2" "HAPERT" "HARMON" "HARMONIC_NUMBER" "HGAP" "HKICK" "HYPERVOLREFERENCE" "IDEALIZED" "INITIAL_OPTIMIZATION" "IMAGENAME" "INFO" "INITIALPOPULATION" "INPUT" "INPUTMOUNITS" "INTENSITYCUT" "INTERPL" "IS_CLOSED" "ITSOLVER" "K0" "K0S" "K1" "K1S" "K2" "K2S" "K3" "K3S" "KEYWORD" "KICK" "KN" "KS" "L" "LAG" "LASERPROFFN" "LATTICE_PHIINIT" "LATTICE_RINIT" "LATTICE_THETAINIT" "LENGTH" "LEVEL" "LOGBENDTRAJECTORY" "LOWER" "LOWERBOUND" "MASS" "MATERIAL" "MAXGENERATIONS" "MAXITERS" "MAXR" "MAXSTEPS" "MAXZ" "METHOD" "MINR" "MINZ" "MODE" "MREX" "MREY" "MSCALX" "MSCALY" "MT" "MUTATION_PROBABILITY" "MX" "MY" "NBIN" "NFREQ" "NLEFT" "NO" "NPART" "NPOINTS" "NRIGHT" "NUMCELLS" "NUM_COWORKERS" "NUM_IND_GEN" "NUM_MASTERS" "OBJECTIVES" "OFFSETPX" "OFFSETPY" "OFFSETPZ" "OFFSETT" "OFFSETX" "OFFSETY" "OFFSETZ" "ONE_PILOT_CONVERGE" "OPCHARGE" "OPMASS" "OPYIELD" "ORDER" "ORIENTATION" "ORIGIN" "OUTDIR" "OUTFN" "OUTPUT" "P1" "P2" "P3" "P4" "PARFFTT" "PARFFTX" "PARFFTY" "PARTICLE" "PARTICLEMATTERINTERACTION" "PATTERN" "PC" "PDIS" "PHI" "PHI0" "PHIINIT" "PHIMAX" "PHIMIN" "PLANE" "POLYORDER" "PRECMODE" "PRINIT" "PSDUMPEACHTURN" "PSDUMPFRAME" "PSDUMPFREQ" "PSI" "PTC" "PYMULT" "PZINIT" "PZMULT" "R51" "R52" "R61" "R62" "RADIUS" "RANDOM" "RANGE" "RASTER" "RECOMBINATION_PROBABILITY" "REFER" "REFPOS" "REPARTFREQ" "RESET" "RFFREQ" "RFMAPFN" "RFPHI" "RINIT" "RMAX" "RMIN" "ROTATION" "ROW" "S" "SAMPLINGS" "SCALABLE" "SEED" "SELECTED" "SEQUENCE" "SIGMA" "SIGMAPX" "SIGMAPY" "SIGMAPZ" "SIGMAR" "SIGMAT" "SIGMAX" "SIGMAY" "SIGMAZ" "SIGX" "SIGY" "SIMBIN_CROSSOVER_NU" "SIMTMPDIR" "SLPTC" "SOL_SYNCH" "SPLIT" "SPTDUMPFREQ" "STARTPOPULATION" "STATDUMPFREQ" "STEP" "STEPSPERTURN" "STOP" "STRING" "SUPERPOSE" "SYMMETRY" "T0" "TABLE" "TAU" "TELL" "TEMPLATEDIR" "TFALL" "THETA" "THIN" "THRESHOLD" "TIME" "TIMEINTEGRATOR" "TMULT" "TO" "TOL" "TOLERANCE" "TPULSEFWHM" "TRACE" "TRIMCOILTHRESHOLD" "TRISE" "TURNS" "TYPE" "UPPER" "UPPERBOUND" "VARIABLE" "VERIFY" "VERSION" "VKICK" "VMAX" "VMIN" "VOLT" "W" "WAKEF" "WARN" "WARP" "WEIGHT" "WIDTH" "WRITETOFILE" "X" "XEND" "XMA" "XMULT" "XSIZE" "XSTART" "Y" "YEND" "YMA" "YMULT" "YSIZE" "YSTART" "Z" "Z0" "ZEND" "ZINIT" "ZSTART" "ZSTOP") t) "\\>")
+;(concat "\\<" (regexp-opt '("A" "ADD" "AMPLITUDE_MODEL" "ANGLE" "APERTURE" "APVETO" "AT" "AUTOPHASE" "B" "BBOXINCR" "BCFFTZ" "BCFFTT" "BCFFTX" "BCFFTY" "BCURRENT" "BEAM_PHIINIT" "BEAM_PRINIT" "BEAM_RINIT" "BFREQ" "BIRTH_CONTROL" "BMAX" "BOUNDPDESTROYFQ" "BSCALE" "BY" "CALLS" "CATHTEMP" "CENTRE" "CHARGE" "CLASS" "CLEAR" "CMD" "COEFDENOM" "COEFDENOMPHI" "COEFNUM" "COEFNUMPHI" "COLUMN" "CONDUCT" "CONSTRAINTS" "CONST_LENGTH" "CONV_HVOL_PROG" "CORRX" "CORRY" "CORRZ" "CSRDUMP" "CROSSOVER" "CUTOFFLONG" "CUTOFFPX" "CUTOFFPY" "CUTOFFPZ" "CUTOFFR" "CUTOFFX" "CUTOFFY" "CYHARMON" "CZERO" "DESIGNENERGY" "DK1" "DK1S" "DK2" "DK2S" "DKN" "DKNR" "DKS" "DKSR" "DLAG" "DPHI" "DPSI" "DS" "DT" "DTHETA" "DUMP" "DUMP_DAT" "DUMP_FREQ" "DUMP_OFFSPRING" "DVARS" "DVOLT" "DX" "DY" "DZ" "E1" "E2" "EBDUMP" "ECHO" "EKIN" "ELASER" "ELEMEDGE" "ELEMENT" "EMISSIONMODEL" "EMISSIONSTEPS" "EMITTED" "ENABLEHDF5" "ENABLERUTHERFORD" "ENERGY" "ENDSEQUENCE" "END_NORMAL_X" "END_NORMAL_Y" "END_POSITION_X" "END_POSITION_Y" "EPSILON" "ESCALE" "ET" "EVERYSTEP" "EX" "EXPECTED_HYPERVOL" "EXPR" "EY" "FE" "FGEOM" "FIELDMAPDIR" "FIELDMAPDIR" "FILE" "FINT" "FMAPFN" "FNAME" "FORM" "FREQ" "FREQUENCY_MODEL" "FROM" "FSTYPE" "FTOSCAMPLITUDE" "FTOSCPERIODS" "FULL" "GAMMA" "GAP" "GAPWIDTH" "GENE_MUTATION_PROBABILITY" "GEOMETRY" "GREENSF" "H1" "H2" "HAPERT" "HARMON" "HARMONIC_NUMBER" "HGAP" "HKICK" "HYPERVOLREFERENCE" "IDEALIZED" "INITIAL_OPTIMIZATION" "IMAGENAME" "INFO" "INITIALPOPULATION" "INPUT" "INPUTMOUNITS" "INTENSITYCUT" "INTERPL" "IS_CLOSED" "ITSOLVER" "K0" "K0S" "K1" "K1S" "K2" "K2S" "K3" "K3S" "KEYWORD" "KICK" "KN" "KS" "L" "LAG" "LASERPROFFN" "LATTICE_PHIINIT" "LATTICE_RINIT" "LATTICE_THETAINIT" "LENGTH" "LEVEL" "LOGBENDTRAJECTORY" "LOWER" "LOWERBOUND" "MASS" "MATERIAL" "MAXGENERATIONS" "MAXITERS" "MAXR" "MAXSTEPS" "MAXZ" "METHOD" "MINR" "MINZ" "MODE" "MREX" "MREY" "MSCALX" "MSCALY" "MT" "MUTATION_PROBABILITY" "MX" "MY" "NBIN" "NFREQ" "NLEFT" "NO" "NPART" "NPOINTS" "NRIGHT" "NUMCELLS" "NUM_COWORKERS" "NUM_IND_GEN" "NUM_MASTERS" "OBJECTIVES" "OFFSETPX" "OFFSETPY" "OFFSETPZ" "OFFSETT" "OFFSETX" "OFFSETY" "OFFSETZ" "ONE_PILOT_CONVERGE" "OPCHARGE" "OPMASS" "OPYIELD" "ORDER" "ORIENTATION" "ORIGIN" "OUTDIR" "OUTFN" "OUTPUT" "P1" "P2" "P3" "P4" "PARFFTT" "PARFFTX" "PARFFTY" "PARTICLE" "PARTICLEMATTERINTERACTION" "PATTERN" "PC" "PDIS" "PHI" "PHI0" "PHIINIT" "PHIMAX" "PHIMIN" "PLANE" "POLYORDER" "PRECMODE" "PRINIT" "PSDUMPEACHTURN" "PSDUMPFRAME" "PSDUMPFREQ" "PSI" "PTC" "PYMULT" "PZINIT" "PZMULT" "R51" "R52" "R61" "R62" "RADIUS" "RANDOM" "RANGE" "RASTER" "RECOMBINATION_PROBABILITY" "REFER" "REFPOS" "REPARTFREQ" "RESET" "RFFREQ" "RFMAPFN" "RFPHI" "RINIT" "RMAX" "RMIN" "ROTATION" "ROW" "S" "SAMPLINGS" "SCALABLE" "SEED" "SELECTED" "SEQUENCE" "SIGMA" "SIGMAPX" "SIGMAPY" "SIGMAPZ" "SIGMAR" "SIGMAT" "SIGMAX" "SIGMAY" "SIGMAZ" "SIGX" "SIGY" "SIMBIN_CROSSOVER_NU" "SIMTMPDIR" "SLPTC" "SOL_SYNCH" "SPLIT" "SPTDUMPFREQ" "STARTPOPULATION" "STATDUMPFREQ" "STEP" "STEPSPERTURN" "STOP" "STRING" "SUPERPOSE" "SYMMETRY" "T0" "TABLE" "TAU" "TELL" "TEMPLATEDIR" "TFALL" "THETA" "THIN" "THRESHOLD" "TIME" "TIMEINTEGRATOR" "TMULT" "TO" "TOL" "TOLERANCE" "TPULSEFWHM" "TRACE" "TRACKBACK" "TRIMCOILTHRESHOLD" "TRISE" "TURNS" "TYPE" "UPPER" "UPPERBOUND" "VARIABLE" "VERIFY" "VERSION" "VKICK" "VMAX" "VMIN" "VOLT" "W" "WAKEF" "WARN" "WARP" "WEIGHT" "WIDTH" "WRITETOFILE" "X" "XEND" "XMA" "XMULT" "XSIZE" "XSTART" "Y" "YEND" "YMA" "YMULT" "YSIZE" "YSTART" "Z" "Z0" "ZEND" "ZINIT" "ZSTART" "ZSTOP") t) "\\>")
(defconst opal-font-lock-keywords-parameters
(list
- '("\\<\\(A\\(?:DD\\|MPLITUDE_MODEL\\|NGLE\\|P\\(?:ERTURE\\|VETO\\)\\|T\\|UTOPHASE\\)\\|B\\(?:BOXINCR\\|C\\(?:FFT[TXYZ]\\|URRENT\\)\\|EAM_\\(?:\\(?:P\\(?:HI\\|R\\)\\|R\\)INIT\\)\\|FREQ\\|IRTH_CONTROL\\|MAX\\|OUNDPDESTROYFQ\\|SCALE\\|Y\\)\\|C\\(?:A\\(?:LLS\\|THTEMP\\)\\|ENTRE\\|HARGE\\|L\\(?:ASS\\|EAR\\)\\|MD\\|O\\(?:EF\\(?:DENOM\\(?:PHI\\)?\\|NUM\\(?:PHI\\)?\\)\\|LUMN\\|N\\(?:DUCT\\|ST\\(?:RAINTS\\|_LENGTH\\)\\|V_HVOL_PROG\\)\\|RR[XYZ]\\)\\|ROSSOVER\\|SRDUMP\\|UTOFF\\(?:LONG\\|P[XYZ]\\|[RXY]\\)\\|YHARMON\\|ZERO\\)\\|D\\(?:ESIGNENERGY\\|K\\(?:1S\\|2S\\|[NS]R\\|[12NS]\\)\\|LAG\\|P\\(?:[HS]I\\)\\|THETA\\|UMP\\(?:_\\(?:DAT\\|FREQ\\|OFFSPRING\\)\\)?\\|V\\(?:ARS\\|OLT\\)\\|[STXYZ]\\)\\|E\\(?:BDUMP\\|CHO\\|KIN\\|L\\(?:ASER\\|EME\\(?:DGE\\|NT\\)\\)\\|MI\\(?:SSION\\(?:MODEL\\|STEPS\\)\\|TTED\\)\\|N\\(?:ABLE\\(?:HDF5\\|RUTHERFORD\\)\\|D\\(?:SEQUENCE\\|_\\(?:NORMAL_[XY]\\|POSITION_[XY]\\)\\)\\|ERGY\\)\\|PSILON\\|SCALE\\|VERYSTEP\\|XP\\(?:ECTED_HYPERVOL\\|R\\)\\|[12TXY]\\)\\|F\\(?:E\\|GEOM\\|I\\(?:ELDMAPDIR\\|LE\\|NT\\)\\|MAPFN\\|NAME\\|ORM\\|R\\(?:EQ\\(?:UENCY_MODEL\\)?\\|OM\\)\\|STYPE\\|TOSC\\(?:AMPLITUDE\\|PERIODS\\)\\|ULL\\)\\|G\\(?:A\\(?:MMA\\|P\\(?:WIDTH\\)?\\)\\|E\\(?:\\(?:NE_MUTATION_PROBABILIT\\|OMETR\\)Y\\)\\|REENSF\\)\\|H\\(?:A\\(?:PERT\\|RMON\\(?:IC_NUMBER\\)?\\)\\|GAP\\|KICK\\|YPERVOLREFERENCE\\|[12]\\)\\|I\\(?:DEALIZED\\|MAGENAME\\|N\\(?:FO\\|ITIAL\\(?:\\(?:POPUL\\|_OPTIMIZ\\)ATION\\)\\|PUT\\(?:MOUNITS\\)?\\|TE\\(?:NSITYCUT\\|RPL\\)\\)\\|S_CLOSED\\|TSOLVER\\)\\|K\\(?:0S\\|1S\\|2S\\|3S\\|EYWORD\\|ICK\\|[0-3NS]\\)\\|L\\(?:A\\(?:G\\|SERPROFFN\\|TTICE_\\(?:\\(?:PHI\\|R\\|THETA\\)INIT\\)\\)\\|E\\(?:NGTH\\|VEL\\)\\|O\\(?:GBENDTRAJECTORY\\|WER\\(?:BOUND\\)?\\)\\)\\|M\\(?:A\\(?:SS\\|TERIAL\\|X\\(?:GENERATIONS\\|ITERS\\|STEPS\\|[RZ]\\)\\)\\|ETHOD\\|IN[RZ]\\|ODE\\|RE[XY]\\|SCAL[XY]\\|UTATION_PROBABILITY\\|[TXY]\\)\\|N\\(?:BIN\\|FREQ\\|LEFT\\|O\\|P\\(?:ART\\|OINTS\\)\\|RIGHT\\|UM\\(?:CELLS\\|_\\(?:COWORKERS\\|IND_GEN\\|MASTERS\\)\\)\\)\\|O\\(?:BJECTIVES\\|FFSET\\(?:P[XYZ]\\|[TXYZ]\\)\\|NE_PILOT_CONVERGE\\|P\\(?:CHARGE\\|MASS\\|YIELD\\)\\|R\\(?:DER\\|I\\(?:\\(?:ENTATIO\\|GI\\)N\\)\\)\\|UT\\(?:DIR\\|FN\\|PUT\\)\\)\\|P\\(?:A\\(?:R\\(?:FFT[TXY]\\|TICLE\\(?:MATTERINTERACTION\\)?\\)\\|TTERN\\)\\|DIS\\|HI\\(?:0\\|INIT\\|M\\(?:AX\\|IN\\)\\)?\\|LANE\\|OLYORDER\\|R\\(?:ECMODE\\|INIT\\)\\|S\\(?:DUMP\\(?:EACHTURN\\|FR\\(?:AME\\|EQ\\)\\)\\|I\\)\\|TC\\|\\(?:YMUL\\|Z\\(?:INI\\|MUL\\)\\)T\\|[1-4C]\\)\\|R\\(?:5[12]\\|6[12]\\|A\\(?:DIUS\\|N\\(?:DOM\\|GE\\)\\|STER\\)\\|E\\(?:COMBINATION_PROBABILITY\\|F\\(?:ER\\|POS\\)\\|PARTFREQ\\|SET\\)\\|F\\(?:FREQ\\|MAPFN\\|PHI\\)\\|INIT\\|M\\(?:AX\\|IN\\)\\|O\\(?:TATION\\|W\\)\\)\\|S\\(?:AMPLINGS\\|CALABLE\\|E\\(?:ED\\|LECTED\\|QUENCE\\)\\|I\\(?:G\\(?:MA\\(?:P[XYZ]\\|[RTXYZ]\\)?\\|[XY]\\)\\|M\\(?:BIN_CROSSOVER_NU\\|TMPDIR\\)\\)\\|LPTC\\|OL_SYNCH\\|P\\(?:LIT\\|TDUMPFREQ\\)\\|T\\(?:A\\(?:RTPOPULATION\\|TDUMPFREQ\\)\\|EP\\(?:SPERTURN\\)?\\|OP\\|RING\\)\\|UPERPOSE\\|YMMETRY\\)\\|T\\(?:A\\(?:BLE\\|U\\)\\|E\\(?:LL\\|MPLATEDIR\\)\\|FALL\\|H\\(?:ETA\\|IN\\|RESHOLD\\)\\|IME\\(?:INTEGRATOR\\)?\\|MULT\\|OL\\(?:ERANCE\\)?\\|PULSEFWHM\\|R\\(?:ACE\\|I\\(?:MCOILTHRESHOLD\\|SE\\)\\)\\|URNS\\|YPE\\|[0O]\\)\\|UPPER\\(?:BOUND\\)?\\|V\\(?:ARIABLE\\|ER\\(?:IFY\\|SION\\)\\|KICK\\|M\\(?:AX\\|IN\\)\\|OLT\\)\\|W\\(?:A\\(?:KEF\\|R[NP]\\)\\|EIGHT\\|IDTH\\|RITETOFILE\\)\\|X\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Y\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Z\\(?:0\\|END\\|INIT\\|ST\\(?:ART\\|OP\\)\\)\\|[ABLSW-Z]\\)\\>"
+ '("\\<\\(A\\(?:DD\\|MPLITUDE_MODEL\\|NGLE\\|P\\(?:ERTURE\\|VETO\\)\\|T\\|UTOPHASE\\)\\|B\\(?:BOXINCR\\|C\\(?:FFT[TXYZ]\\|URRENT\\)\\|EAM_\\(?:\\(?:P\\(?:HI\\|R\\)\\|R\\)INIT\\)\\|FREQ\\|IRTH_CONTROL\\|MAX\\|OUNDPDESTROYFQ\\|SCALE\\|Y\\)\\|C\\(?:A\\(?:LLS\\|THTEMP\\)\\|ENTRE\\|HARGE\\|L\\(?:ASS\\|EAR\\)\\|MD\\|O\\(?:EF\\(?:DENOM\\(?:PHI\\)?\\|NUM\\(?:PHI\\)?\\)\\|LUMN\\|N\\(?:DUCT\\|ST\\(?:RAINTS\\|_LENGTH\\)\\|V_HVOL_PROG\\)\\|RR[XYZ]\\)\\|ROSSOVER\\|SRDUMP\\|UTOFF\\(?:LONG\\|P[XYZ]\\|[RXY]\\)\\|YHARMON\\|ZERO\\)\\|D\\(?:ESIGNENERGY\\|K\\(?:1S\\|2S\\|[NS]R\\|[12NS]\\)\\|LAG\\|P\\(?:[HS]I\\)\\|THETA\\|UMP\\(?:_\\(?:DAT\\|FREQ\\|OFFSPRING\\)\\)?\\|V\\(?:ARS\\|OLT\\)\\|[STXYZ]\\)\\|E\\(?:BDUMP\\|CHO\\|KIN\\|L\\(?:ASER\\|EME\\(?:DGE\\|NT\\)\\)\\|MI\\(?:SSION\\(?:MODEL\\|STEPS\\)\\|TTED\\)\\|N\\(?:ABLE\\(?:HDF5\\|RUTHERFORD\\)\\|D\\(?:SEQUENCE\\|_\\(?:NORMAL_[XY]\\|POSITION_[XY]\\)\\)\\|ERGY\\)\\|PSILON\\|SCALE\\|VERYSTEP\\|XP\\(?:ECTED_HYPERVOL\\|R\\)\\|[12TXY]\\)\\|F\\(?:E\\|GEOM\\|I\\(?:ELDMAPDIR\\|LE\\|NT\\)\\|MAPFN\\|NAME\\|ORM\\|R\\(?:EQ\\(?:UENCY_MODEL\\)?\\|OM\\)\\|STYPE\\|TOSC\\(?:AMPLITUDE\\|PERIODS\\)\\|ULL\\)\\|G\\(?:A\\(?:MMA\\|P\\(?:WIDTH\\)?\\)\\|E\\(?:\\(?:NE_MUTATION_PROBABILIT\\|OMETR\\)Y\\)\\|REENSF\\)\\|H\\(?:A\\(?:PERT\\|RMON\\(?:IC_NUMBER\\)?\\)\\|GAP\\|KICK\\|YPERVOLREFERENCE\\|[12]\\)\\|I\\(?:DEALIZED\\|MAGENAME\\|N\\(?:FO\\|ITIAL\\(?:\\(?:POPUL\\|_OPTIMIZ\\)ATION\\)\\|PUT\\(?:MOUNITS\\)?\\|TE\\(?:NSITYCUT\\|RPL\\)\\)\\|S_CLOSED\\|TSOLVER\\)\\|K\\(?:0S\\|1S\\|2S\\|3S\\|EYWORD\\|ICK\\|[0-3NS]\\)\\|L\\(?:A\\(?:G\\|SERPROFFN\\|TTICE_\\(?:\\(?:PHI\\|R\\|THETA\\)INIT\\)\\)\\|E\\(?:NGTH\\|VEL\\)\\|O\\(?:GBENDTRAJECTORY\\|WER\\(?:BOUND\\)?\\)\\)\\|M\\(?:A\\(?:SS\\|TERIAL\\|X\\(?:GENERATIONS\\|ITERS\\|STEPS\\|[RZ]\\)\\)\\|ETHOD\\|IN[RZ]\\|ODE\\|RE[XY]\\|SCAL[XY]\\|UTATION_PROBABILITY\\|[TXY]\\)\\|N\\(?:BIN\\|FREQ\\|LEFT\\|O\\|P\\(?:ART\\|OINTS\\)\\|RIGHT\\|UM\\(?:CELLS\\|_\\(?:COWORKERS\\|IND_GEN\\|MASTERS\\)\\)\\)\\|O\\(?:BJECTIVES\\|FFSET\\(?:P[XYZ]\\|[TXYZ]\\)\\|NE_PILOT_CONVERGE\\|P\\(?:CHARGE\\|MASS\\|YIELD\\)\\|R\\(?:DER\\|I\\(?:\\(?:ENTATIO\\|GI\\)N\\)\\)\\|UT\\(?:DIR\\|FN\\|PUT\\)\\)\\|P\\(?:A\\(?:R\\(?:FFT[TXY]\\|TICLE\\(?:MATTERINTERACTION\\)?\\)\\|TTERN\\)\\|DIS\\|HI\\(?:0\\|INIT\\|M\\(?:AX\\|IN\\)\\)?\\|LANE\\|OLYORDER\\|R\\(?:ECMODE\\|INIT\\)\\|S\\(?:DUMP\\(?:EACHTURN\\|FR\\(?:AME\\|EQ\\)\\)\\|I\\)\\|TC\\|\\(?:YMUL\\|Z\\(?:INI\\|MUL\\)\\)T\\|[1-4C]\\)\\|R\\(?:5[12]\\|6[12]\\|A\\(?:DIUS\\|N\\(?:DOM\\|GE\\)\\|STER\\)\\|E\\(?:COMBINATION_PROBABILITY\\|F\\(?:ER\\|POS\\)\\|PARTFREQ\\|SET\\)\\|F\\(?:FREQ\\|MAPFN\\|PHI\\)\\|INIT\\|M\\(?:AX\\|IN\\)\\|O\\(?:TATION\\|W\\)\\)\\|S\\(?:AMPLINGS\\|CALABLE\\|E\\(?:ED\\|LECTED\\|QUENCE\\)\\|I\\(?:G\\(?:MA\\(?:P[XYZ]\\|[RTXYZ]\\)?\\|[XY]\\)\\|M\\(?:BIN_CROSSOVER_NU\\|TMPDIR\\)\\)\\|LPTC\\|OL_SYNCH\\|P\\(?:LIT\\|TDUMPFREQ\\)\\|T\\(?:A\\(?:RTPOPULATION\\|TDUMPFREQ\\)\\|EP\\(?:SPERTURN\\)?\\|OP\\|RING\\)\\|UPERPOSE\\|YMMETRY\\)\\|T\\(?:A\\(?:BLE\\|U\\)\\|E\\(?:LL\\|MPLATEDIR\\)\\|FALL\\|H\\(?:ETA\\|IN\\|RESHOLD\\)\\|IME\\(?:INTEGRATOR\\)?\\|MULT\\|OL\\(?:ERANCE\\)?\\|PULSEFWHM\\|R\\(?:AC\\(?:E\\|KBACK\\)\\|I\\(?:MCOILTHRESHOLD\\|SE\\)\\)\\|URNS\\|YPE\\|[0O]\\)\\|UPPER\\(?:BOUND\\)?\\|V\\(?:ARIABLE\\|ER\\(?:IFY\\|SION\\)\\|KICK\\|M\\(?:AX\\|IN\\)\\|OLT\\)\\|W\\(?:A\\(?:KEF\\|R[NP]\\)\\|EIGHT\\|IDTH\\|RITETOFILE\\)\\|X\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Y\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Z\\(?:0\\|END\\|INIT\\|ST\\(?:ART\\|OP\\)\\)\\|[ABLSW-Z]\\)\\>"
. font-lock-variable-name-face)
)
"Highlighting expressions for OPAL mode (parameters).")
diff --git a/tools/opal2sdds/main.cpp b/tools/opal2sdds/main.cpp
index 0eb589d675d18ac507279ece13892ee527b5f681..ad8e4637c997381dd9c423821ab29f759d106638 100644
--- a/tools/opal2sdds/main.cpp
+++ b/tools/opal2sdds/main.cpp
@@ -37,9 +37,13 @@ enum FORMAT {
typedef h5_file_t file_t;
-data_t readStepData(file_t file);
+data_t readStepData(file_t file, bool temporal);
attributes_t readStepAttributes(file_t file);
-void readH5HutFile(const std::string &fname, size_t step, data_t &data, attributes_t &attr);
+void readH5HutFile(const std::string &fname,
+ size_t step,
+ bool temporal,
+ data_t &data,
+ attributes_t &attr);
void convertToElegantUnits(data_t &data);
void writeSDDSFile(const std::string &fname, const data_t &data, const attributes_t &attr, FORMAT form);
void printInfo(const std::string &input);
@@ -50,6 +54,7 @@ int main(int argc, char **argv) {
std::string inputFile(""), outputFile("/dev/stdout");
bool doPrintInfo = false;
+ bool temporal = false;
int step = -1;
FORMAT form = ASCII;
@@ -75,6 +80,8 @@ int main(int argc, char **argv) {
return 0;
} else if (std::string(argv[i]) == "--info") {
doPrintInfo = true;
+ } else if (std::string(argv[i]) == "--temporal") {
+ temporal = true;
} else {
inputFile = std::string(argv[i]);
}
@@ -94,7 +101,7 @@ int main(int argc, char **argv) {
data_t data;
attributes_t attr;
- readH5HutFile(inputFile, step, data, attr);
+ readH5HutFile(inputFile, step, temporal, data, attr);
convertToElegantUnits(data);
writeSDDSFile(outputFile, data, attr, form);
@@ -106,7 +113,11 @@ void reportOnError(h5_int64_t rc, const char* file, int line) {
std::cerr << "H5 rc= " << rc << " in " << file << " @ line " << line << std::endl;
}
-void readH5HutFile(const std::string &fname, size_t step, data_t &data, attributes_t &attr) {
+void readH5HutFile(const std::string &fname,
+ size_t step,
+ bool temporal,
+ data_t &data,
+ attributes_t &attr) {
h5_prop_t props = H5CreateFileProp ();
MPI_Comm comm = MPI_COMM_WORLD;
H5SetPropFileMPIOCollective (props, &comm);
@@ -118,12 +129,12 @@ void readH5HutFile(const std::string &fname, size_t step, data_t &data, attribut
REPORTONERROR(H5SetStep(file, readStep));
- data = readStepData(file);
+ data = readStepData(file, temporal);
attr = readStepAttributes(file);
REPORTONERROR(H5CloseFile(file));
}
-data_t readStepData(file_t file) {
+data_t readStepData(file_t file, bool temporal) {
data_t data;
h5_ssize_t numParticles = H5PartGetNumParticles(file);
@@ -144,11 +155,19 @@ data_t readStepData(file_t file) {
}
data.insert(std::make_pair(std::string("y"), dData));
- READDATA(Float64, file, "z", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
- dData[n] = f64buffer[n];
+ if (temporal) {
+ READDATA(Float64, file, "time", f64buffer);
+ for(long int n = 0; n < numParticles; ++ n) {
+ dData[n] = f64buffer[n];
+ }
+ data.insert(std::make_pair(std::string("time"), dData));
+ } else {
+ READDATA(Float64, file, "z", f64buffer);
+ for(long int n = 0; n < numParticles; ++ n) {
+ dData[n] = f64buffer[n];
+ }
+ data.insert(std::make_pair(std::string("z"), dData));
}
- data.insert(std::make_pair(std::string("z"), dData));
READDATA(Float64, file, "px", f64buffer);
for(long int n = 0; n < numParticles; ++ n) {
@@ -201,22 +220,32 @@ attributes_t readStepAttributes(file_t file) {
}
void convertToElegantUnits(data_t &data) {
- const size_t size = data["x"].size();
+ size_t size;
+ bool temporal = false;
+
+ try {
+ size = data["z"].size();
+ } catch (...) {
+ temporal = true;
+ size = data["x"].size();
+ }
for (size_t i = 0; i < size; ++ i) {
data["px"][i] /= data["pz"][i];
data["py"][i] /= data["pz"][i];
+ if (temporal) continue;
+
data["z"][i] /= -299792458.0;
}
}
void writeSDDSFile(const std::string &fname, const data_t &data, const attributes_t &attr, FORMAT form) {
- const std::map<std::string, std::string> nameConversion {{"x", "x"},
- {"y", "y"},
- {"t", "z"},
- {"xp", "px"},
- {"yp", "py"},
- {"p", "pz"}};
+ std::map<std::string, std::string> nameConversion {{"x", "x"},
+ {"y", "y"},
+ {"t", "z"},
+ {"xp", "px"},
+ {"yp", "py"},
+ {"p", "pz"}};
const std::map<std::string, std::string> nameUnits {{"x", "m"},
{"y", "m"},
{"t", "s"},
@@ -237,6 +266,13 @@ void writeSDDSFile(const std::string &fname, const data_t &data, const attribute
char buffer1[64];
char buffer2[8];
+ for (const auto entry: data) {
+ if (entry.first == "time") {
+ nameConversion["t"] = "time";
+ break;
+ }
+ }
+
strcpy(buffer0, fname.c_str());
strcpy(buffer1, "extracted OPAL data");
if (SDDS_InitializeOutput(&SDDS_dataset, form, 1, buffer1, NULL, buffer0 ) != 1) {
@@ -368,6 +404,7 @@ void printUsage(char **argv) {
<< indent << std::setw(width) << std::left << "--binary" << "= whether SDDS output should be in \n"
<< indent << std::setw(width + 2) << std::left << " " << "binary format. Default: ASCII\n"
<< indent << std::setw(width) << std::left << "--info" << "= printing path length corresponding to step\n"
+ << indent << std::setw(width) << std::left << "--time-of-arrival" << "= extract and save the time of arrival instead of the longitudinal position\n"
<< indent << std::setw(width) << std::left << "--help" << "= this help\n"
<< std::endl;
-}
+}
\ No newline at end of file
diff --git a/tools/sdds2opal/main.cpp b/tools/sdds2opal/main.cpp
index fe8017c6cae27aa1dc07442083ab211062550429..91d7097dedb0c00a0d80430953a7ffb357714ee0 100644
--- a/tools/sdds2opal/main.cpp
+++ b/tools/sdds2opal/main.cpp
@@ -11,11 +11,12 @@
#include <cstring>
#include <unistd.h>
-std::vector<std::vector<double> > readSDDSFile(std::string fname);
+std::vector<std::vector<double> > readSDDSFile(std::string fname, bool temporal);
void printUsage(char **argv);
int main(int argc, char **argv) {
std::string inputFile(""), outputFile("/dev/stdout");
+ bool temporal = false;
if (argc == 1) {
printUsage(argv);
@@ -30,6 +31,8 @@ int main(int argc, char **argv) {
inputFile = std::string(argv[++ i]);
} else if (std::string(argv[i]) == "--output" && i + 1 < argc) {
outputFile = std::string(argv[++ i]);
+ } else if (std::string(argv[i]) == "--temporal") {
+ temporal = true;
} else if (std::string(argv[i]) == "--help") {
printUsage(argv);
return 0;
@@ -47,7 +50,7 @@ int main(int argc, char **argv) {
}
}
- auto data = readSDDSFile(inputFile);
+ auto data = readSDDSFile(inputFile, temporal);
std::ofstream out(outputFile);
out.precision(8);
out << data.size() << "\n";
@@ -61,7 +64,7 @@ int main(int argc, char **argv) {
return 0;
}
-std::vector<std::vector<double> > readSDDSFile(std::string fname) {
+std::vector<std::vector<double> > readSDDSFile(std::string fname, bool temporal) {
SDDS_DATASET SDDS_dataset;
void *columnData;
@@ -122,6 +125,8 @@ std::vector<std::vector<double> > readSDDSFile(std::string fname) {
for (std::vector<double> &row: fileData) {
row[1] *= row[5];
row[3] *= row[5];
+ if (temporal) continue;
+
row[4] *= -299792458.0;
}
@@ -147,6 +152,7 @@ void printUsage(char **argv) {
<< "Options\n\n"
<< indent << std::setw(width) << std::left << "--input <path to input>" << "= path to input file\n"
<< indent << std::setw(width) << std::left << "--output <file name>" << "= name of output. If omitted, directed\n"
+ << indent << std::setw(width) << std::left << "--time-of-arrival" << "= extract and save the time of arrival instead of the longitudinal position\n"
<< indent << std::setw(width + 2) << " " << "to stdout\n"
<< indent << std::setw(width) << std::left << "--help" << "= this help\n"
<< std::endl;