Commit 41047061 authored by Christof Metzger-Kraus's avatar Christof Metzger-Kraus
Browse files

fixing bug that prevented the OrbitThreader from setting the design energy of elements

parent 42d23d3c
......@@ -85,22 +85,23 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
optimizedPhase = status.first;
finalEnergy = status.second;
AstraPhase = std::fmod(optimizedPhase + Physics::pi / 2, Physics::two_pi);
newPhase = std::fmod(originalPhase + optimizedPhase + Physics::two_pi, Physics::two_pi);
element->setPhasem(newPhase);
element->setAutophaseVeto();
OpalData::getInstance()->setMaxPhase(itsCavity_m->getName(), optimizedPhase);
double frequency = element->getFrequencym();
double basePhase = std::fmod(frequency * (t + tErr), 1.0);
basePhase = std::fmod(optimizedPhase - basePhase + Physics::two_pi, Physics::two_pi);
INFOMSG(itsCavity_m->getName() << "_phi = " << optimizedPhase * Physics::rad2deg << " [deg], "
INFOMSG(itsCavity_m->getName() << "_phi = " << std::fmod(newPhase + basePhase, Physics::two_pi) * Physics::rad2deg << " [deg], "
<< "corresp. in Astra = " << AstraPhase * Physics::rad2deg << " [deg],\n"
<< "E = " << finalEnergy << " [MeV], " << "phi_nom = " << originalPhase * Physics::rad2deg << " [deg]\n"
<< "Ez_0 = " << amplitude << " [MV/m]" << "\n"
<< "time = " << (t + tErr) * 1e9 << " [ns], dt = " << dt * 1e12 << " [ps]" << endl);
OpalData::getInstance()->setMaxPhase(itsCavity_m->getName(), optimizedPhase);
return optimizedPhase;
}
......@@ -133,7 +134,11 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
double originalPhase = element->getPhasem();
if (element->getAutophaseVeto()) {
std::pair<double, double> status(originalPhase, Util::getEnergy(initialP_m, itsReference_m.getM() * 1e-6));
double frequency = element->getFrequencym();
double basePhase = std::fmod(frequency * t, 1.0);
double phase = std::fmod(originalPhase - basePhase + Physics::two_pi, Physics::two_pi);
double E = track(initialR_m, initialP_m, t, dt, phase);
std::pair<double, double> status(-basePhase, E);//Util::getEnergy(initialP_m, itsReference_m.getM() * 1e-6));
return status;
}
......@@ -186,6 +191,8 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
Phimax = std::fmod(initialPhase + Physics::two_pi, Physics::two_pi);
E = track(initialR_m, initialP_m, t, dt, Phimax + originalPhase);
*gmsg << level1 << __DBGMSG__ << element->getName() << ": phase = " << Phimax + originalPhase << ", energy = " << E << endl;
std::pair<double, double> status(Phimax, E);
return status;
}
......
......@@ -358,8 +358,8 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
for (; it != end; ++ it) {
std::shared_ptr<Component> element = (*it).getElement();
if (visitedElements.find(element->getName()) == visitedElements.end() &&
!(element->getType() != ElementBase::RFCAVITY ||
element->getType() != ElementBase::TRAVELINGWAVE)) {
!(element->getType() == ElementBase::RFCAVITY ||
element->getType() == ElementBase::TRAVELINGWAVE)) {
element->setDesignEnergy(kineticEnergyeV);
}
......
......@@ -298,7 +298,7 @@ private:
void computeWakefield(IndexMap::value_t &elements);
void computeParticleMatterInteraction(IndexMap::value_t elements, OrbitThreader &oth);
void computeSpaceChargeFields(unsigned long long step);
void prepareOpalBeamlineSections();
// void prepareOpalBeamlineSections();
void dumpStats(long long step, bool psDump, bool statDump);
void setOptionalVariables();
bool hasEndOfLineReached();
......
......@@ -315,7 +315,7 @@ void Bend::initialise(PartBunch *bunch,
startField_m -= elementEdge_m;
endField_m -= elementEdge_m;
elementEdge_m = 0.0;
// elementEdge_m = 0.0;
msg << level2
<< "======================================================================\n"
......@@ -1432,65 +1432,97 @@ void Bend::retrieveDesignEnergy(double startField) {
// }
}
std::pair<Vector_t, Vector_t> Bend::getDesignPathSecant(double startsAtDistFromEdge, double length) const {
std::vector<Vector_t> Bend::getPartialDesignPath(double startsAtDistFromEdge, double length) const {
length = std::abs(length);
Vector_t startPosition(0.0);
Vector_t tangent(0, 0, 1);
if (length <= 0.0 || startsAtDistFromEdge < 0.0)
throw GeneralClassicException("Bend::getDesignPath",
"both length and path length from edge have to be positive");
double pathLength = refTrajMap_m[0](2);
unsigned int size = refTrajMap_m.size();
const unsigned int size = refTrajMap_m.size();
for (unsigned int i = 1; i < size; ++ i) {
Vector_t step = refTrajMap_m[i] - refTrajMap_m[i-1];
double stepSize = euclidian_norm(step);
if (pathLength + stepSize > startsAtDistFromEdge) {
double diff = startsAtDistFromEdge - pathLength;
tangent = step / stepSize;
startPosition = refTrajMap_m[i-1] + diff * tangent;
unsigned int firstIdx = 1;
while (firstIdx < size && refTrajMap_m[firstIdx](2) < 0.0) {
++ firstIdx;
}
if (firstIdx == size)
throw GeneralClassicException("Bend::getDesignPath",
"never crossing edge");
if (length <= 1e-12) {
return std::make_pair(startPosition, tangent);
}
Vector_t step = refTrajMap_m[firstIdx] - refTrajMap_m[firstIdx - 1];
double stepSize = euclidian_norm(step);
double minPathLength = (refTrajMap_m[firstIdx](2) /
(refTrajMap_m[firstIdx](2) - refTrajMap_m[firstIdx - 1](2))) * stepSize;
double pathLength = 0.0;
std::vector<Vector_t> designPath;
for (unsigned j = i; j < size; ++ j) {
Vector_t position = refTrajMap_m[j];
unsigned int idx = firstIdx + 1;
if (startsAtDistFromEdge < minPathLength) {
pathLength = (minPathLength - startsAtDistFromEdge);
if (euclidian_norm(position - startPosition) >= length) {
step = refTrajMap_m[j-1] - refTrajMap_m[j];
double tau = (dot(startPosition - position, step) - sqrt(std::pow(dot(startPosition - position, step), 2) - dot(step, step) * (dot(startPosition - position, startPosition - position) - std::pow(length, 2)))) / dot(step, step);
if (pathLength > length) {
designPath.push_back(refTrajMap_m[firstIdx] -
pathLength / stepSize * step);
designPath.push_back(designPath[0] + length / stepSize * step);
tangent = position + tau * step - startPosition;
tangent /= euclidian_norm(tangent);
transformFrom(designPath);
return std::make_pair(startPosition, tangent);
}
}
return designPath;
} else {
designPath.push_back(refTrajMap_m[firstIdx] -
(minPathLength - startsAtDistFromEdge) / stepSize * step);
designPath.push_back(refTrajMap_m[firstIdx]);
}
} else {
pathLength = minPathLength;
const CoordinateSystemTrafo fromEndToExitRegion(Vector_t(0, 0, exitParameter2_m),
Quaternion(1, 0, 0, 0));
const CoordinateSystemTrafo toExitRegion = (fromEndToExitRegion *
getBeginToEnd_local());
while (idx < size && pathLength < startsAtDistFromEdge) {
step = refTrajMap_m[idx] - refTrajMap_m[idx - 1];
stepSize = euclidian_norm(step);
pathLength += stepSize;
++ idx;
}
if (idx == size || toExitRegion.transformTo(refTrajMap_m[idx - 1])(2) > 0.0)
throw GeneralClassicException("Bend::getDesignPath",
"start position outside dipole");
Vector_t position = refTrajMap_m[size - 1];
Vector_t step = position - refTrajMap_m[size - 2];
step /= euclidian_norm(step);
double tau = (dot(startPosition - position, step) + sqrt(std::pow(dot(startPosition - position, step), 2) - dot(step, step) * (dot(startPosition - position, startPosition - position) - std::pow(length, 2)))) / dot(step, step);
-- idx;
designPath.push_back(refTrajMap_m[idx] -
(pathLength - startsAtDistFromEdge) / stepSize * step);
designPath.push_back(refTrajMap_m[idx]);
pathLength -= startsAtDistFromEdge;
tangent = position + tau * step - startPosition;
tangent /= euclidian_norm(tangent);
if (pathLength > length) {
designPath.back() -= (pathLength - length) / stepSize * step;
return std::make_pair(startPosition, tangent);
transformFrom(designPath);
return designPath;
}
++ idx;
}
while (idx < size && pathLength < length) {
designPath.push_back(refTrajMap_m[idx]);
step = refTrajMap_m[idx] - refTrajMap_m[idx - 1];
stepSize = euclidian_norm(step);
pathLength += stepSize;
++ idx;
}
double diff = startsAtDistFromEdge - pathLength;
unsigned int i = size - 1;
Vector_t step = refTrajMap_m[i] - refTrajMap_m[i-1];
double stepSize = euclidian_norm(step);
tangent = step / stepSize;
startPosition = refTrajMap_m[i] + diff * tangent;
tangent = tangent;
designPath.back() -= (pathLength - length) / stepSize * step;
transformFrom(designPath);
return std::make_pair(startPosition, tangent);
return designPath;
}
std::vector<Vector_t> Bend::getOutline() const {
......
......@@ -120,7 +120,7 @@ public:
double getExitAngle() const;
double getMapLength() const;
std::pair<Vector_t, Vector_t> getDesignPathSecant(double startsAtDistFromEdge, double length) const;
std::vector<Vector_t> getPartialDesignPath(double startsAtDistFromEdge, double length) const;
std::vector<Vector_t> getOutline() const;
MeshData getSurfaceMesh() const;
......@@ -192,6 +192,8 @@ private:
CoordinateSystemTrafo getBeginToEnd_local() const;
void transformFrom(std::vector<Vector_t> & points) const;
std::string messageHeader_m;
BorisPusher pusher_m; /// Pusher used to integrate reference particle
......@@ -379,4 +381,11 @@ CoordinateSystemTrafo Bend::getBeginToEnd_local() const
return beginToEnd_lcs_m;
}
inline
void Bend::transformFrom(std::vector<Vector_t> & points) const {
for (auto it = points.begin(); it != points.end(); ++ it) {
*it = getCSTrafoGlobal2Local().transformFrom(*it);
}
}
#endif // CLASSIC_BEND_H
\ No newline at end of file
......@@ -26,6 +26,8 @@
#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
#include "Solvers/ParticleMaterInteractionHandler.hh"
#include "Utilities/Util.h"
#include <memory>
extern Inform *gmsg;
......@@ -68,7 +70,9 @@ Collimator::Collimator():
pitch_m(0.0),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
parmatint_m(NULL),
isWarping_m(true),
warpCurve_m(NULL)
{}
......@@ -107,7 +111,9 @@ Collimator::Collimator(const Collimator &right):
pitch_m(right.pitch_m),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
parmatint_m(NULL),
isWarping_m(right.isWarping_m),
warpCurve_m(NULL)
{
setGeom();
}
......@@ -148,12 +154,14 @@ Collimator::Collimator(const std::string &name):
pitch_m(0.0),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
parmatint_m(NULL),
isWarping_m(true),
warpCurve_m(NULL)
{}
Collimator::~Collimator() {
if(online_m)
if (online_m)
goOffline();
}
......@@ -170,8 +178,8 @@ inline bool Collimator::isInColl(Vector_t R, Vector_t P, double recpgamma) {
*/
const double z = R(2) + P(2) * recpgamma;
if((z > 0.0) && (z <= getElementLength())) {
if(isAPepperPot_m) {
if ((z > 0.0) && (z <= getElementLength())) {
if (isAPepperPot_m) {
/**
------------
......@@ -190,24 +198,24 @@ inline bool Collimator::isInColl(Vector_t R, Vector_t P, double recpgamma) {
particles they are not lost.
*/
const double h = pitch_m;
const double h = pitch_m;
const double xL = - 0.5 * h * (nHolesX_m - 1);
const double yL = - 0.5 * h * (nHolesY_m - 1);
bool alive = false;
for(unsigned int m = 0; (m < nHolesX_m && (!alive)); m++) {
for(unsigned int n = 0; (n < nHolesY_m && (!alive)); n++) {
double x_m = xL + (m * h);
double y_m = yL + (n * h);
for (unsigned int m = 0; (m < nHolesX_m && (!alive)); m++) {
for (unsigned int n = 0; (n < nHolesY_m && (!alive)); n++) {
double x_m = xL + (m * h);
double y_m = yL + (n * h);
/** are we in a) ? */
double rr = std::pow((R(0) - x_m) / rHole_m, 2) + std::pow((R(1) - y_m) / rHole_m, 2);
alive = (rr < 1.0);
}
}
return !alive;
} else if(isASlit_m || isARColl_m) {
} else if (isASlit_m || isARColl_m) {
return (std::abs(R(0)) >= getXsize() || std::abs(R(1)) >= getYsize());
} else if(isAWire_m) {
} else if (isAWire_m) {
ERRORMSG("Not yet implemented");
} else {
// case of an elliptic collimator
......@@ -221,10 +229,10 @@ bool Collimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &
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 / sqrt(1.0 + dot(P, P));
const double recpgamma = Physics::c * dt / sqrt(1.0 + dot(P, P));
bool pdead = isInColl(R, P, recpgamma);
if(pdead) {
if (pdead) {
if (lossDs_m) {
double frac = -R(2) / P(2) * recpgamma;
lossDs_m->addParticle(R, P,
......@@ -239,7 +247,7 @@ bool Collimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &
bool Collimator::applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
const double dt = RefPartBunch_m->getdT();
const double recpgamma = Physics::c * dt / sqrt(1.0 + dot(P, P));
const double recpgamma = Physics::c * dt / sqrt(1.0 + dot(P, P));
return isInColl(R, P, recpgamma);
}
......@@ -249,9 +257,9 @@ bool Collimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
bool isDead = false;
if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if(r(2) < zend_m && r(2) > zstart_m ) {
if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if ( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if (r(2) < zend_m && r(2) > zstart_m ) {
int pflag = checkPoint(r(0), r(1));
isDead = (pflag != 0);
}
......@@ -276,16 +284,16 @@ bool Collimator::checkCollimator(PartBunch &bunch, const int turnnumber, const d
boundingSphere.first = 0.5 * (rmax + rmin);
boundingSphere.second = euclidian_norm(rmax - boundingSphere.first);
if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
// if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if( r1 > r_start - 100.0 && r1 < r_end + 100.0 ){
if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
// if ( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if ( r1 > r_start - 100.0 && r1 < r_end + 100.0 ){
size_t tempnum = bunch.getLocalNum();
int pflag = 0;
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch.PType[i] == ParticleType::REGULAR && bunch.R[i](2) < zend_m && bunch.R[i](2) > zstart_m ) {
for (unsigned int i = 0; i < tempnum; ++i) {
if (bunch.PType[i] == ParticleType::REGULAR && bunch.R[i](2) < zend_m && bunch.R[i](2) > zstart_m ) {
pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
/// bunch.Bin[i] != -1 makes sure the partcile is not stored in more than one collimator
if ((pflag != 0) && (bunch.Bin[i] != -1)) {
if ((pflag != 0) && (bunch.Bin[i] != -1)) {
if (!parmatint_m)
lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
bunch.Bin[i] = -1;
......@@ -335,7 +343,7 @@ void Collimator::initialise(PartBunch *bunch) {
void Collimator::finalise()
{
if(online_m)
if (online_m)
goOffline();
*gmsg << "* Finalize probe" << endl;
}
......@@ -355,11 +363,11 @@ void Collimator::goOnline(const double &) {
}
void Collimator::print() {
if(RefPartBunch_m == NULL) {
if(!informed_m) {
if (RefPartBunch_m == NULL) {
if (!informed_m) {
std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
ERRORMSG(errormsg << endl);
if(Ippl::myNode() == 0) {
if (Ippl::myNode() == 0) {
ofstream omsg("errormsg.txt", ios_base::app);
omsg << errormsg << endl;
omsg.close();
......@@ -370,20 +378,20 @@ void Collimator::print() {
}
*gmsg << level3;
if(isAPepperPot_m)
if (isAPepperPot_m)
*gmsg << "* Pepperpot x= " << a_m << " y= " << b_m << " r= " << rHole_m
<< " nx= " << nHolesX_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
else if(isASlit_m)
else if (isASlit_m)
*gmsg << "* Slit x= " << getXsize() << " Slit y= " << getYsize()
<< " fn= " << filename_m << endl;
else if(isARColl_m)
else if (isARColl_m)
*gmsg << "* RCollimator a= " << getXsize() << " b= " << getYsize()
<< " fn= " << filename_m
<< " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
else if(isACColl_m)
else if (isACColl_m)
*gmsg << "* CCollimator angle start " << xstart_m << " (Deg) angle end " << ystart_m << " (Deg) "
<< "R start " << xend_m << " (mm) R rend " << yend_m << " (mm)" << endl;
else if(isAWire_m)
else if (isAWire_m)
*gmsg << "* Wire x= " << x0_m << " y= " << y0_m << endl;
else
*gmsg << "* ECollimator a= " << getXsize() << " b= " << b_m << " fn= " << filename_m
......@@ -422,15 +430,15 @@ ElementBase::ElementType Collimator::getType() const {
}
string Collimator::getCollimatorShape() {
if(isAPepperPot_m)
if (isAPepperPot_m)
return "PeperPot";
else if(isASlit_m)
else if (isASlit_m)
return "Slit";
else if(isARColl_m)
else if (isARColl_m)
return "RCollimator";
else if(isACColl_m)
else if (isACColl_m)
return "CCollimator";
else if(isAWire_m)
else if (isAWire_m)
return "Wire";
else
return "ECollimator";
......@@ -454,7 +462,7 @@ void Collimator::setGeom() {
geom_m[1].y = ystart_m - halfdist / coeff2;
geom_m[2].x = xend_m + halfdist * coeff1;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[3].x = xend_m - halfdist * coeff1;
geom_m[3].y = yend_m + halfdist / coeff2;
......@@ -471,16 +479,82 @@ void Collimator::setGeom() {
int Collimator::checkPoint(const double &x, const double &y) {
int cn = 0;
int cn = 0;
for(int i = 0; i < 4; i++) {
if(((geom_m[i].y <= y) && (geom_m[i+1].y > y))
|| ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
for (int i = 0; i < 4; i++) {
if (((geom_m[i].y <= y) && (geom_m[i+1].y > y))
|| ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
if(x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
if (x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
++cn;
}
}
return (cn & 1); // 0 if even (out), and 1 if odd (in)
}
void Collimator::setWarpCurve(const std::vector<Vector_t> & curve) {
const size_t size = curve.size();
unsigned int every = 1;
double error = 0.0;
gsl_interp_accel *accel = gsl_interp_accel_alloc();
do {
every *= 2;
bool addLast = ((size - 1) % every > 0 || size < every);
size_t reducedSize = (size - 1) / every + (addLast? 2: 1);
if (reducedSize < 3) break;
std::vector<double> zvalues(reducedSize);
std::vector<double> xvalues(reducedSize);
size_t j = 0;
for (size_t i = 0; i < size; i += every, ++ j) {
zvalues[j] = curve[i](2);
xvalues[j] = curve[i](0);
}
if (addLast) {
zvalues[reducedSize - 1] = curve.back()(2);
xvalues[reducedSize - 1] = curve.back()(0);
}
gsl_spline *interpolant = gsl_spline_alloc(gsl_interp_cspline, reducedSize);
gsl_spline_init(interpolant, zvalues.data(), xvalues.data(), reducedSize);
error = 0.0;
for (size_t i = 0; i < size; ++ i) {
double x = gsl_spline_eval(interpolant, curve[i](2), accel);
double localError = x - curve[i](0);
error += std::pow(localError, 2);
}
error = sqrt(error / size);
} while (error < 1e-9 && every * 2 < size);
every /= 2;
bool addLast = ((size - 1) % every > 0 || size < every);
size_t reducedSize = (size - 1) / every + (addLast? 2: 1);
std::vector<double> zvalues(reducedSize);
std::vector<double> xvalues(reducedSize);
size_t j = 0;
for (size_t i = 0; i < size; i += every, ++ j) {
zvalues[j] = curve[i](2);
xvalues[j] = curve[i](0);
}
if (addLast) {
zvalues[reducedSize - 1] = curve.back()(2);
xvalues[reducedSize - 1] = curve.back()(0);
}
gsl_spline *interpolant = gsl_spline_alloc(gsl_interp_cspline, reducedSize);
gsl_spline_init(interpolant, zvalues.data(), xvalues.data(), reducedSize);
}
\ No newline at end of file
......@@ -24,6 +24,9 @@
#include "AbsBeamline/Component.h"
#include "gsl/gsl_spline.h"
#include "gsl/gsl_interp.h"
#include <vector>
class BeamlineVisitor;
......@@ -148,7 +151,11 @@ public:
virtual bool isInColl(Vector_t R, Vector_t P, double recpgamma);
int checkPoint(const double & x, const double & y );
int checkPoint(const double & x, const double & y );
bool getWarpFlag() const;
void setWarpFlag(bool warp = true);
void setWarpCurve(const std::vector<Vector_t> &curve);
private:
// Not implemented.
......@@ -202,6 +209,8 @@ private: