From dff9c1a56a406ecd87608c7575591989dd039a30 Mon Sep 17 00:00:00 2001 From: Christof Kraus <christof.kraus@psi.ch> Date: Thu, 9 May 2019 14:16:08 +0200 Subject: [PATCH] Let OrbitThreader find next element even if zstop is inbetween two elements. For this try to estimate length of longest implicit drift. Then search for twice this distance for further elements. --- src/Algorithms/IndexMap.cpp | 6 +- src/Algorithms/OrbitThreader.cpp | 97 +++++++++++++--- src/Algorithms/OrbitThreader.h | 7 +- .../AbsBeamline/{Bend.cpp => Bend2D.cpp} | 106 +++++++++--------- src/Classic/AbsBeamline/{Bend.h => Bend2D.h} | 74 ++++++------ src/Classic/AbsBeamline/BendBase.h | 4 +- src/Classic/AbsBeamline/CMakeLists.txt | 6 +- src/Classic/AbsBeamline/RBend.cpp | 10 +- src/Classic/AbsBeamline/RBend.h | 6 +- src/Classic/AbsBeamline/RBend3D.h | 8 +- src/Classic/AbsBeamline/SBend.cpp | 6 +- src/Classic/AbsBeamline/SBend.h | 6 +- src/Classic/Solvers/CSRIGFWakeFunction.cpp | 15 +-- src/Classic/Solvers/CSRWakeFunction.cpp | 6 +- src/Classic/Structure/MeshGenerator.cpp | 4 +- src/Elements/OpalBeamline.cpp | 10 +- src/Elements/OpalElement.cpp | 2 +- .../AbsBeamline/DipoleFieldTest.cpp | 16 +-- 18 files changed, 228 insertions(+), 161 deletions(-) rename src/Classic/AbsBeamline/{Bend.cpp => Bend2D.cpp} (96%) rename src/Classic/AbsBeamline/{Bend.h => Bend2D.h} (88%) diff --git a/src/Algorithms/IndexMap.cpp b/src/Algorithms/IndexMap.cpp index f1bb31b09..330b424c9 100644 --- a/src/Algorithms/IndexMap.cpp +++ b/src/Algorithms/IndexMap.cpp @@ -6,7 +6,7 @@ #include "Algorithms/IndexMap.h" #include "AbstractObjects/OpalData.h" #include "AbsBeamline/Multipole.h" -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "Utilities/Util.h" #include "Physics/Physics.h" #include "OPALconfig.h" @@ -24,6 +24,8 @@ IndexMap::IndexMap(): { } void IndexMap::print(std::ostream &out) const { + if (mapRange2Element_m.size() == 0) return; + out << "Size of map " << mapRange2Element_m.size() << " sections " << std::endl; out << std::fixed << std::setprecision(6); auto mapIti = mapRange2Element_m.begin(); @@ -252,7 +254,7 @@ void IndexMap::saveSDDS(double startS) const { case ElementBase::RBEND: case ElementBase::SBEND: { - const Bend* bend = static_cast<const Bend*>(element.get()); + const Bend2D* bend = static_cast<const Bend2D*>(element.get()); if (bend->getRotationAboutZ() > 0.5 * Physics::pi && bend->getRotationAboutZ() < 1.5 * Physics::pi) { items[DIPOLE] = -1; diff --git a/src/Algorithms/OrbitThreader.cpp b/src/Algorithms/OrbitThreader.cpp index f2ea1d95d..a66bbb56b 100644 --- a/src/Algorithms/OrbitThreader.cpp +++ b/src/Algorithms/OrbitThreader.cpp @@ -2,6 +2,7 @@ #include "Algorithms/CavityAutophaser.h" #include "AbsBeamline/RFCavity.h" +#include "AbsBeamline/BendBase.h" #include "AbsBeamline/TravelingWave.h" #include "AbstractObjects/OpalData.h" @@ -80,12 +81,17 @@ OrbitThreader::OrbitThreader(const PartData &ref, void OrbitThreader::execute() { double initialPathLength = pathLength_m; - trackBack(); + double maxDistance = computeMaximalImplicitDrift(); + + auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY); + std::set<std::string> visitedElements; + + trackBack(2 * maxDistance); + Vector_t nextR = r_m / (Physics::c * dt_m); integrator_m.push(nextR, p_m, dt_m); nextR *= Physics::c * dt_m; - auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY); - std::set<std::string> visitedElements; + setDesignEnergy(allElements, visitedElements); auto elementSet = itsOpalBeamline_m.getElements(nextR); @@ -98,7 +104,7 @@ void OrbitThreader::execute() { double initialS = pathLength_m; Vector_t initialR = r_m; Vector_t initialP = p_m; - integrate(elementSet, maxIntegSteps_m); + integrate(elementSet, maxIntegSteps_m, 2 * maxDistance); registerElement(elementSet, initialS, initialR, initialP); @@ -123,19 +129,21 @@ void OrbitThreader::execute() { elementSet = itsOpalBeamline_m.getElements(nextR); } - } while ((time_m - initialTime_m) / dt_m < maxIntegSteps_m && - errorFlag_m != HITMATERIAL && + } while (errorFlag_m != HITMATERIAL && errorFlag_m != EOL); - // imap_m.tidyUp(); + imap_m.tidyUp(); *gmsg << level1 << "\n" << imap_m << endl; imap_m.saveSDDS(initialPathLength); processElementRegister(); } -void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxSteps) { +void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxSteps, double maxDrift) { static size_t step = 0; + CoordinateSystemTrafo labToBeamline = itsOpalBeamline_m.getCSTrafoLab2Local(); + const double oldPathLength = pathLength_m; + double stepLength; Vector_t nextR; do { errorFlag_m = EVERYTHINGFINE; @@ -165,7 +173,10 @@ void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxStep Bf += itsOpalBeamline_m.rotateFromLocalCS(*it, localB); } - if (step % loggingFrequency_m == 0 && Ippl::myNode() == 0 && !OpalData::getInstance()->isOptimizerRun()) { + if (pathLength_m > 0.0 && + pathLength_m < zstop_m && + step % loggingFrequency_m == 0 && Ippl::myNode() == 0 && + !OpalData::getInstance()->isOptimizerRun()) { logger_m << std::setw(18) << std::setprecision(8) << pathLength_m + std::copysign(euclidean_norm(r_m - oldR), dt_m) << std::setw(18) << std::setprecision(8) << r_m(0) << std::setw(18) << std::setprecision(8) << r_m(1) @@ -198,12 +209,15 @@ void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t maxStep integrator_m.push(nextR, p_m, dt_m); nextR *= Physics::c * dt_m; - if (pathLength_m > zstop_m) { + if ((activeSet.size() == 0 && std::abs(pathLength_m - oldPathLength) > maxDrift) || + (activeSet.size() > 0 && pathLength_m > zstop_m)) { errorFlag_m = EOL; return; } - } while (((time_m - initialTime_m) / dt_m < maxSteps) && + 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))); } @@ -263,7 +277,7 @@ double OrbitThreader::getMaxDesignEnergy(const IndexMap::value_t &elementSet) co return designEnergy; } -void OrbitThreader::trackBack() { +void OrbitThreader::trackBack(double maxDrift) { dt_m *= -1; double initialPathLength = pathLength_m; @@ -271,10 +285,11 @@ void OrbitThreader::trackBack() { 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) { auto elementSet = itsOpalBeamline_m.getElements(nextR); - integrate(elementSet, 1000); + integrate(elementSet, 1000, maxDrift); nextR = r_m / (Physics::c * dt_m); integrator_m.push(nextR, p_m, dt_m); @@ -368,3 +383,59 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std:: } } } + +double OrbitThreader::computeMaximalImplicitDrift() { + FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY); + double maxDrift = 0.0; + FieldList::iterator it = allElements.begin(); + const FieldList::iterator end = allElements.end(); + + for (; it != end; ++ it) { + auto element1 = it->getElement(); + const auto &toEdge = element1->getCSTrafoGlobal2Local(); + auto toEnd = element1->getEdgeToEnd() * toEdge; + Vector_t end1 = toEnd.transformFrom(Vector_t(0, 0, 0)); + Vector_t directionEnd = toEnd.rotateFrom(Vector_t(0, 0, 1)); + if (element1->getType() == ElementBase::RBEND || + element1->getType() == ElementBase::SBEND || + element1->getType() == ElementBase::RBEND3D ) { + auto toBegin = element1->getEdgeToBegin() * toEdge; + + BendBase *bend = static_cast<BendBase*>(element1.get()); + double angleRelativeToFace = bend->getEntranceAngle() - bend->getBendAngle(); + directionEnd = toBegin.rotateFrom(Vector_t(sin(angleRelativeToFace), 0, cos(angleRelativeToFace))); + } + + double minDistanceLocal = std::numeric_limits<double>::max(); + FieldList::iterator it2 = allElements.begin(); + for (; it2 != end; ++ it2) { + if (it == it2) continue; + + auto element2 = it2->getElement(); + const auto &toEdge = element2->getCSTrafoGlobal2Local(); + auto toBegin = element2->getEdgeToBegin() * toEdge; + Vector_t begin2 = toBegin.transformFrom(Vector_t(0, 0, 0)); + Vector_t directionBegin = toBegin.rotateFrom(Vector_t(0, 0, 1)); + if (element2->getType() == ElementBase::RBEND || + element2->getType() == ElementBase::SBEND || + element2->getType() == ElementBase::RBEND3D ) { + BendBase *bend = static_cast<BendBase*>(element2.get()); + double E1 = bend->getEntranceAngle(); + directionBegin = toBegin.rotateFrom(Vector_t(sin(E1), 0, cos(E1))); + } + + + double distance = euclidean_norm(begin2 - end1); + double directionProjection = dot(directionEnd, directionBegin); + if (directionProjection > 0.999 && minDistanceLocal > distance) { + minDistanceLocal = distance; + } + } + + if (maxDrift < minDistanceLocal) { + maxDrift = minDistanceLocal; + } + } + + return maxDrift; +} \ No newline at end of file diff --git a/src/Algorithms/OrbitThreader.h b/src/Algorithms/OrbitThreader.h index b9f51da5c..6a78739ed 100644 --- a/src/Algorithms/OrbitThreader.h +++ b/src/Algorithms/OrbitThreader.h @@ -69,8 +69,8 @@ private: std::multimap<std::string, elementPosition> elementRegistry_m; - void trackBack(); - void integrate(const IndexMap::value_t &activeSet, size_t maxSteps); + void trackBack(double maxDrift); + void integrate(const IndexMap::value_t &activeSet, size_t maxSteps, double maxDrift = 10.0); bool containsCavity(const IndexMap::value_t &activeSet); void autophaseCavities(const IndexMap::value_t &activeSet, const std::set<std::string> &visitedElements); double getMaxDesignEnergy(const IndexMap::value_t &elementSet) const; @@ -78,6 +78,7 @@ private: void registerElement(const IndexMap::value_t &elementSet, double, const Vector_t &r, const Vector_t &p); void processElementRegister(); void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements); + double computeMaximalImplicitDrift(); }; inline @@ -97,4 +98,4 @@ IndexMap::value_t OrbitThreader::getTouchingElements(const std::pair<double, dou return imap_m.getTouchingElements(range); } -#endif +#endif \ No newline at end of file diff --git a/src/Classic/AbsBeamline/Bend.cpp b/src/Classic/AbsBeamline/Bend2D.cpp similarity index 96% rename from src/Classic/AbsBeamline/Bend.cpp rename to src/Classic/AbsBeamline/Bend2D.cpp index db872e392..eb85aa78a 100644 --- a/src/Classic/AbsBeamline/Bend.cpp +++ b/src/Classic/AbsBeamline/Bend2D.cpp @@ -1,12 +1,12 @@ // ------------------------------------------------------------------------ -// $RCSfile: Bend.cpp,v $ +// $RCSfile: Bend2D.cpp,v $ // ------------------------------------------------------------------------ // $Revision: 1.1.1.1 $ // ------------------------------------------------------------------------ // Copyright: see Copyright.readme // ------------------------------------------------------------------------ // -// Definitions for class: Bend +// Definitions for class: Bend2D // Defines the abstract interface for a sector bend magnet. // // ------------------------------------------------------------------------ @@ -18,7 +18,7 @@ // // ------------------------------------------------------------------------ -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "Algorithms/PartBunchBase.h" #include "AbsBeamline/BeamlineVisitor.h" #include "Utilities/Options.h" @@ -34,10 +34,10 @@ extern Inform *gmsg; -// Class Bend +// Class Bend2D // ------------------------------------------------------------------------ -Bend::Bend(): +Bend2D::Bend2D(): BendBase(), messageHeader_m(" * "), pusher_m(), @@ -78,7 +78,7 @@ Bend::Bend(): } -Bend::Bend(const Bend &right): +Bend2D::Bend2D(const Bend2D &right): BendBase(right), messageHeader_m(" * "), pusher_m(right.pusher_m), @@ -119,7 +119,7 @@ Bend::Bend(const Bend &right): } -Bend::Bend(const std::string &name): +Bend2D::Bend2D(const std::string &name): BendBase(name), messageHeader_m(" * "), pusher_m(), @@ -160,7 +160,7 @@ Bend::Bend(const std::string &name): } -Bend::~Bend() { +Bend2D::~Bend2D() { if (entryFieldAccel_m != NULL) { for (unsigned int i = 0; i < 3u; ++ i) { gsl_spline_free(entryFieldValues_m[i]); @@ -188,7 +188,7 @@ Bend::~Bend() { * This function merely repackages the field arrays as type Vector_t and calls * the equivalent method but with the Vector_t data types. */ -bool Bend::apply(const size_t &i, +bool Bend2D::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) { @@ -212,7 +212,7 @@ bool Bend::apply(const size_t &i, return false; } -bool Bend::apply(const Vector_t &R, +bool Bend2D::apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, @@ -235,7 +235,7 @@ bool Bend::apply(const Vector_t &R, } -bool Bend::applyToReferenceParticle(const Vector_t &R, +bool Bend2D::applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, @@ -257,11 +257,11 @@ bool Bend::applyToReferenceParticle(const Vector_t &R, return false; } -void Bend::goOnline(const double &) { +void Bend2D::goOnline(const double &) { online_m = true; } -void Bend::initialise(PartBunchBase<double, 3> *bunch, +void Bend2D::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) { @@ -305,7 +305,7 @@ void Bend::initialise(PartBunchBase<double, 3> *bunch, } } -void Bend::adjustFringeFields(double ratio) { +void Bend2D::adjustFringeFields(double ratio) { findChordLength(*gmsg, chordLength_m); double delta = std::abs(entranceParameter1_m - entranceParameter2_m); @@ -323,7 +323,7 @@ void Bend::adjustFringeFields(double ratio) { setupFringeWidths(); } -double Bend::calculateBendAngle() { +double Bend2D::calculateBendAngle() { const double mass = RefPartBunch_m->getM(); const double gamma = designEnergy_m / mass + 1.0; @@ -365,7 +365,7 @@ double Bend::calculateBendAngle() { } -void Bend::calcEngeFunction(double zNormalized, +void Bend2D::calcEngeFunction(double zNormalized, const std::vector<double> &engeCoeff, int polyOrder, double &engeFunc, @@ -429,7 +429,7 @@ void Bend::calcEngeFunction(double zNormalized, } } -Vector_t Bend::calcCentralField(const Vector_t &R, +Vector_t Bend2D::calcCentralField(const Vector_t &R, double deltaX) { Vector_t B(0, 0, 0); @@ -448,7 +448,7 @@ Vector_t Bend::calcCentralField(const Vector_t &R, return B; } -Vector_t Bend::calcEntranceFringeField(const Vector_t &R, +Vector_t Bend2D::calcEntranceFringeField(const Vector_t &R, double deltaX) { const CoordinateSystemTrafo toEntranceRegion(Vector_t(0, 0, entranceParameter2_m), @@ -485,7 +485,7 @@ Vector_t Bend::calcEntranceFringeField(const Vector_t &R, return toEntranceRegion.rotateFrom(B); } -Vector_t Bend::calcExitFringeField(const Vector_t &R, +Vector_t Bend2D::calcExitFringeField(const Vector_t &R, double deltaX) { const CoordinateSystemTrafo fromEndToExitRegion(Vector_t(0, 0, exitParameter2_m), @@ -521,7 +521,7 @@ Vector_t Bend::calcExitFringeField(const Vector_t &R, return toExitRegion.rotateFrom(B); } -bool Bend::calculateMapField(const Vector_t &R, Vector_t &B) { +bool Bend2D::calculateMapField(const Vector_t &R, Vector_t &B) { B = Vector_t(0.0); bool verticallyInside = (std::abs(R(1)) < 0.5 * gap_m); @@ -584,7 +584,7 @@ bool Bend::calculateMapField(const Vector_t &R, Vector_t &B) { return hitMaterial; } -void Bend::calculateRefTrajectory(double &angleX, double &angleY) { +void Bend2D::calculateRefTrajectory(double &angleX, double &angleY) { const double mass = RefPartBunch_m->getM(); const double gamma = designEnergy_m / mass + 1.; @@ -651,7 +651,7 @@ void Bend::calculateRefTrajectory(double &angleX, double &angleY) { angleX = -atan2(P(0), P(2)) + entranceAngle_m; } -double Bend::estimateFieldAdjustmentStep(double actualBendAngle, +double Bend2D::estimateFieldAdjustmentStep(double actualBendAngle, double mass, double betaGamma) { @@ -669,7 +669,7 @@ double Bend::estimateFieldAdjustmentStep(double actualBendAngle, } -void Bend::findBendEffectiveLength(double startField, double endField) { +void Bend2D::findBendEffectiveLength(double startField, double endField) { /* * Use an iterative procedure to set the width of the @@ -759,7 +759,7 @@ void Bend::findBendEffectiveLength(double startField, double endField) { } } -void Bend::findBendStrength(double mass, +void Bend2D::findBendStrength(double mass, double gamma, double betaGamma, double charge) { @@ -840,7 +840,7 @@ void Bend::findBendStrength(double mass, } } -bool Bend::findIdealBendParameters(double chordLength) { +bool Bend2D::findIdealBendParameters(double chordLength) { double refMass = RefPartBunch_m->getM(); double refGamma = designEnergy_m / refMass + 1.0; @@ -882,7 +882,7 @@ bool Bend::findIdealBendParameters(double chordLength) { return reinitialize; } -void Bend::findReferenceExitOrigin(double &x, double &z) { +void Bend2D::findReferenceExitOrigin(double &x, double &z) { /* * Find x,z coordinates of reference trajectory as it passes exit edge @@ -903,7 +903,7 @@ void Bend::findReferenceExitOrigin(double &x, double &z) { } } -bool Bend::initializeFieldMap(Inform &msg) { +bool Bend2D::initializeFieldMap(Inform &msg) { fieldmap_m = Fieldmap::getFieldmap(fileName_m, fast_m); @@ -918,7 +918,7 @@ bool Bend::initializeFieldMap(Inform &msg) { } -bool Bend::inMagnetCentralRegion(const Vector_t &R) const { +bool Bend2D::inMagnetCentralRegion(const Vector_t &R) const { Vector_t rotationCenter(-designRadius_m * cosEntranceAngle_m, R(1), designRadius_m * sinEntranceAngle_m); double distFromRotCenter = euclidean_norm(R - rotationCenter); @@ -936,14 +936,14 @@ bool Bend::inMagnetCentralRegion(const Vector_t &R) const { return false; } -bool Bend::inMagnetEntranceRegion(const Vector_t &R) const { +bool Bend2D::inMagnetEntranceRegion(const Vector_t &R) const { return (R(2) >= entranceParameter1_m && R(2) < 0.0 && std::abs(R(0)) < aperture_m.second[0]); } -bool Bend::inMagnetExitRegion(const Vector_t &R) const { +bool Bend2D::inMagnetExitRegion(const Vector_t &R) const { Vector_t Rprime = getBeginToEnd_local().transformTo(R); @@ -952,14 +952,14 @@ bool Bend::inMagnetExitRegion(const Vector_t &R) const { std::abs(Rprime(0)) < aperture_m.second[0]); } -bool Bend::isPositionInEntranceField(const Vector_t &R) const { +bool Bend2D::isPositionInEntranceField(const Vector_t &R) const { return (polyOrderEntry_m >= 0 && R(2) >= entranceParameter1_m && R(2) < entranceParameter3_m); } -bool Bend::isPositionInExitField(const Vector_t &R) const { +bool Bend2D::isPositionInExitField(const Vector_t &R) const { Vector_t Rprime = getBeginToEnd_local().transformTo(R); return (polyOrderExit_m >= 0 && @@ -967,7 +967,7 @@ bool Bend::isPositionInExitField(const Vector_t &R) const { Rprime(2) < exitParameter3_m); } -void Bend::print(Inform &msg, double bendAngleX, double bendAngleY) { +void Bend2D::print(Inform &msg, double bendAngleX, double bendAngleY) { msg << level2 << "\n" << "Start of field map: " << startField_m @@ -1067,7 +1067,7 @@ void Bend::print(Inform &msg, double bendAngleX, double bendAngleY) { } -void Bend::readFieldMap(Inform &msg) { +void Bend2D::readFieldMap(Inform &msg) { msg << level2 << getName() << " using file "; fieldmap_m->getInfo(&msg); @@ -1147,7 +1147,7 @@ void Bend::readFieldMap(Inform &msg) { } } -bool Bend::reinitialize() { +bool Bend2D::reinitialize() { setBendStrength(); double bendAngleX = 0.0; @@ -1163,7 +1163,7 @@ bool Bend::reinitialize() { return false; } -void Bend::setBendEffectiveLength(double startField, double endField) { +void Bend2D::setBendEffectiveLength(double startField, double endField) { // Find initial angle. double actualBendAngle = calculateBendAngle(); @@ -1175,7 +1175,7 @@ void Bend::setBendEffectiveLength(double startField, double endField) { } -void Bend::setBendStrength() { +void Bend2D::setBendStrength() { // Estimate bend field magnitude. double mass = RefPartBunch_m->getM(); @@ -1191,7 +1191,7 @@ void Bend::setBendStrength() { findBendStrength(mass, gamma, betaGamma, charge); } -void Bend::setEngeOriginDelta(double delta) { +void Bend2D::setEngeOriginDelta(double delta) { /* * This function is used to shift the perpendicular distance of the * entrance and exit Enge function origins with respect to the entrance @@ -1213,7 +1213,7 @@ void Bend::setEngeOriginDelta(double delta) { setupFringeWidths(); } -void Bend::setFieldCalcParam() { +void Bend2D::setFieldCalcParam() { cosEntranceAngle_m = cos(entranceAngle_m); sinEntranceAngle_m = sin(entranceAngle_m); @@ -1307,7 +1307,7 @@ void Bend::setFieldCalcParam() { maxAngleExitAperture -= rotationCenter; } -void Bend::setGapFromFieldMap() { +void Bend2D::setGapFromFieldMap() { if(gap_m <= 0.0) gap_m = fieldmap_m->getFieldGap(); @@ -1316,7 +1316,7 @@ void Bend::setGapFromFieldMap() { } -bool Bend::setupBendGeometry(Inform &msg, double &startField, double &endField) { +bool Bend2D::setupBendGeometry(Inform &msg, double &startField, double &endField) { chordLength_m = 0.0; if(!findChordLength(msg, chordLength_m)) @@ -1368,7 +1368,7 @@ bool Bend::setupBendGeometry(Inform &msg, double &startField, double &endField) } -bool Bend::setupDefaultFieldMap(Inform &msg) { +bool Bend2D::setupDefaultFieldMap(Inform &msg) { if(length_m <= 0.0) { ERRORMSG("If using \"1DPROFILE1-DEFAULT\" field map you must set the " @@ -1384,7 +1384,7 @@ bool Bend::setupDefaultFieldMap(Inform &msg) { } -void Bend::setFieldBoundaries(double startField, double endField) { +void Bend2D::setFieldBoundaries(double startField, double endField) { startField_m = startField - deltaBeginEntry_m / cos(entranceAngle_m); endField_m = (startField + angle_m * designRadius_m + @@ -1392,14 +1392,14 @@ void Bend::setFieldBoundaries(double startField, double endField) { } -void Bend::setupPusher(PartBunchBase<double, 3> *bunch) { +void Bend2D::setupPusher(PartBunchBase<double, 3> *bunch) { RefPartBunch_m = bunch; pusher_m.initialise(bunch->getReference()); } -bool Bend::treatAsDrift(Inform &msg, double chordLength) { +bool Bend2D::treatAsDrift(Inform &msg, double chordLength) { if(designEnergy_m <= 0.0) { WARNMSG("Warning: bend design energy is zero. Treating as drift." << endl); @@ -1436,7 +1436,7 @@ bool Bend::treatAsDrift(Inform &msg, double chordLength) { return false; } -std::vector<Vector_t> Bend::getOutline() const { +std::vector<Vector_t> Bend2D::getOutline() const { std::vector<Vector_t> outline; Vector_t rotationCenter = Vector_t(-designRadius_m * cosEntranceAngle_m, 0.0, designRadius_m * sinEntranceAngle_m); unsigned int numSteps = 2; @@ -1539,7 +1539,7 @@ std::vector<Vector_t> Bend::getOutline() const { return outline; } -MeshData Bend::getSurfaceMesh() const { +MeshData Bend2D::getSurfaceMesh() const { MeshData mesh; const Vector_t hgap(0, 0.5 * getFullGap(), 0); std::vector<Vector_t> outline = getOutline(); @@ -1767,7 +1767,7 @@ MeshData Bend::getSurfaceMesh() const { return mesh; } -bool Bend::isInside(const Vector_t &R) const { +bool Bend2D::isInside(const Vector_t &R) const { if (std::abs(R(1)) > gap_m) return false; if (inMagnetEntranceRegion(R)) { @@ -1781,25 +1781,25 @@ bool Bend::isInside(const Vector_t &R) const { return (std::abs(R(1)) < 0.5 * gap_m && inMagnetCentralRegion(R)); } -void Bend::setupFringeWidths() +void Bend2D::setupFringeWidths() { widthEntranceFringe_m = 2 * std::min(entranceParameter3_m - entranceParameter1_m, aperture_m.second[0]) + aperture_m.second[0]; widthExitFringe_m = 2 * std::min(exitParameter3_m - exitParameter1_m, aperture_m.second[0]) + aperture_m.second[0]; } //set the number of slices for map tracking -void Bend::setNSlices(const std::size_t& nSlices) { +void Bend2D::setNSlices(const std::size_t& nSlices) { nSlices_m = nSlices; } //get the number of slices for map tracking -std::size_t Bend::getNSlices() const { +std::size_t Bend2D::getNSlices() const { return nSlices_m; } /// Get entrance fringe field length. // Used to create fringe fields in ThickTracker, (before edge[m], after edge[m]) -std::array<double,2> Bend::getEntranceFringeFieldLength() const { +std::array<double,2> Bend2D::getEntranceFringeFieldLength() const { double entrParam1, entrParam2, entrParam3; fieldmap_m->get1DProfile1EntranceParam(entrParam1, @@ -1819,7 +1819,7 @@ std::array<double,2> Bend::getEntranceFringeFieldLength() const { /// Get exit fringe field length. // Used to create fringe fields in ThickTracker, (after edge[m], before edge[m]) -std::array<double,2> Bend::getExitFringeFieldLength() const { +std::array<double,2> Bend2D::getExitFringeFieldLength() const { std::array<double,2> extFFL; double exitParam1, exitParam2, exitParam3; diff --git a/src/Classic/AbsBeamline/Bend.h b/src/Classic/AbsBeamline/Bend2D.h similarity index 88% rename from src/Classic/AbsBeamline/Bend.h rename to src/Classic/AbsBeamline/Bend2D.h index 9f58ebdb3..aa87c1628 100644 --- a/src/Classic/AbsBeamline/Bend.h +++ b/src/Classic/AbsBeamline/Bend2D.h @@ -9,7 +9,7 @@ // Copyright: see Copyright.readme // ------------------------------------------------------------------------ // -// Definitions for class: Bend +// Definitions for class: Bend2D // Defines the abstract interface for a general bend magnet. // // ------------------------------------------------------------------------ @@ -39,7 +39,7 @@ class Fieldmap; class MeshData; /* - * Class Bend + * Class Bend2D * * Interface for general bend magnet. * @@ -47,18 +47,18 @@ class MeshData; * */ -class Bend: public BendBase { +class Bend2D: public BendBase { public: /// Constructor with given name. - explicit Bend(const std::string &name); + explicit Bend2D(const std::string &name); - Bend(); - Bend(const Bend &); - virtual ~Bend(); + Bend2D(); + Bend2D(const Bend2D &); + virtual ~Bend2D(); - /// Apply visitor to Bend. + /// Apply visitor to Bend2D. virtual void accept(BeamlineVisitor &) const = 0; /* @@ -107,21 +107,15 @@ public: double getStartElement() const; - void setExitAngle(double exitAngle); - - /* - * Set the name of the field map. - * - * For now this means a file that contains Enge function coefficients - * that describe the fringe fields at the entrance and exit. - */ /// Set quadrupole field component. void setK1(double k1); void resetReinitializeFlag(); void resetRecalcRefTrajFlag(); - double getExitAngle() const; + void setExitAngle(double exitAngle); + virtual double getExitAngle() const; + double getMapLength() const; std::vector<Vector_t> getOutline() const; @@ -158,7 +152,7 @@ private: FRIEND_TEST(Maxwell, Zeros); #endif // Not implemented. - void operator=(const Bend &); + void operator=(const Bend2D &); void adjustFringeFields(double ratio); double calculateBendAngle(); @@ -314,120 +308,120 @@ private: inline -void Bend::finalise() { +void Bend2D::finalise() { online_m = false; } inline -void Bend::getDimensions(double &sBegin, double &sEnd) const { +void Bend2D::getDimensions(double &sBegin, double &sEnd) const { sBegin = startField_m; sEnd = endField_m; } inline -double Bend::getBendRadius() const { +double Bend2D::getBendRadius() const { return designRadius_m; } inline -double Bend::getEffectiveCenter() const { +double Bend2D::getEffectiveCenter() const { return elementEdge_m + designRadius_m * angle_m / 2.0; } inline -double Bend::getEffectiveLength() const { +double Bend2D::getEffectiveLength() const { return designRadius_m * angle_m; } inline -double Bend::getStartElement() const { +double Bend2D::getStartElement() const { return elementEdge_m; } inline -void Bend::setK1(double k1) { +void Bend2D::setK1(double k1) { if (std::abs(k1) > 0.0) { - throw GeneralClassicException("Bend::setK1", + throw GeneralClassicException("Bend2D::setK1", "Quadrupole field temporarily not supported"); } fieldIndex_m = k1; } inline -void Bend::resetReinitializeFlag() { +void Bend2D::resetReinitializeFlag() { reinitialize_m = true; } inline -void Bend::resetRecalcRefTrajFlag() { +void Bend2D::resetRecalcRefTrajFlag() { recalcRefTraj_m = true; } inline -void Bend::setMessageHeader(const std::string & header) +void Bend2D::setMessageHeader(const std::string & header) { messageHeader_m = header; } inline -double Bend::getStartField() const +double Bend2D::getStartField() const { return startField_m; } inline -Fieldmap* Bend::getFieldmap() +Fieldmap* Bend2D::getFieldmap() { return fieldmap_m; } inline -double Bend::getExitAngle() const +double Bend2D::getExitAngle() const { return exitAngle_m; } inline -void Bend::setExitAngle(double angle) +void Bend2D::setExitAngle(double angle) { exitAngle_m = angle; } inline -double Bend::getMapLength() const +double Bend2D::getMapLength() const { return exitParameter2_m - entranceParameter2_m; } inline -CoordinateSystemTrafo Bend::getEdgeToEnd() const +CoordinateSystemTrafo Bend2D::getEdgeToEnd() const { return beginToEnd_m; } inline -CoordinateSystemTrafo Bend::getBeginToEnd_local() const +CoordinateSystemTrafo Bend2D::getBeginToEnd_local() const { return beginToEnd_lcs_m; } inline -void Bend::setCSTrafoToEntranceRegion(const CoordinateSystemTrafo &trafo) { +void Bend2D::setCSTrafoToEntranceRegion(const CoordinateSystemTrafo &trafo) { toEntranceRegion_m = trafo; } inline -void Bend::setCSTrafoToExitRegion(const CoordinateSystemTrafo &trafo) { +void Bend2D::setCSTrafoToExitRegion(const CoordinateSystemTrafo &trafo) { toExitRegion_m = trafo; } inline -Vector_t Bend::transformToEntranceRegion(const Vector_t &R) const { +Vector_t Bend2D::transformToEntranceRegion(const Vector_t &R) const { return toEntranceRegion_m.transformTo(R); } inline -Vector_t Bend::transformToExitRegion(const Vector_t &R) const { +Vector_t Bend2D::transformToExitRegion(const Vector_t &R) const { return toExitRegion_m.transformTo(R); } diff --git a/src/Classic/AbsBeamline/BendBase.h b/src/Classic/AbsBeamline/BendBase.h index c309dc867..fee02884d 100644 --- a/src/Classic/AbsBeamline/BendBase.h +++ b/src/Classic/AbsBeamline/BendBase.h @@ -22,7 +22,7 @@ public: double getBendAngle() const; virtual void setEntranceAngle(double entranceAngle); double getEntranceAngle() const; - + virtual double getExitAngle() const = 0; void setFullGap(double); double getFullGap() const; @@ -144,4 +144,4 @@ std::string BendBase::getFieldMapFN() const { } -#endif +#endif \ No newline at end of file diff --git a/src/Classic/AbsBeamline/CMakeLists.txt b/src/Classic/AbsBeamline/CMakeLists.txt index a7cb13c19..67ebd3db1 100644 --- a/src/Classic/AbsBeamline/CMakeLists.txt +++ b/src/Classic/AbsBeamline/CMakeLists.txt @@ -6,7 +6,7 @@ set (_SRCS AttributeSet.cpp BeamBeam.cpp BeamlineVisitor.cpp - Bend.cpp + Bend2D.cpp BendBase.cpp CCollimator.cpp Component.cpp @@ -63,7 +63,7 @@ set (HDRS AttributeSet.h BeamBeam.h BeamlineVisitor.h - Bend.h + Bend2D.h BendBase.h CCollimator.h Component.h @@ -115,4 +115,4 @@ install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline" #Turn on gcov #set(CMAKE_CXX_FLAGS "-g -O0 -Wall -fprofile-arcs -ftest-coverage") -#set(CMAKE_CXX_OUTPUT_EXTENSION_REPLACE 1) +#set(CMAKE_CXX_OUTPUT_EXTENSION_REPLACE 1) \ No newline at end of file diff --git a/src/Classic/AbsBeamline/RBend.cpp b/src/Classic/AbsBeamline/RBend.cpp index c4facb937..6629b5ddd 100644 --- a/src/Classic/AbsBeamline/RBend.cpp +++ b/src/Classic/AbsBeamline/RBend.cpp @@ -33,19 +33,19 @@ extern Inform *gmsg; // ------------------------------------------------------------------------ RBend::RBend(): - Bend() + Bend2D() { setMessageHeader("RBend "); } RBend::RBend(const RBend &right): - Bend(right) + Bend2D(right) { setMessageHeader("RBend "); } RBend::RBend(const std::string &name): - Bend(name) + Bend2D(name) { setMessageHeader("RBend "); } @@ -136,12 +136,12 @@ ElementBase::ElementType RBend::getType() const { } void RBend::setBendAngle(double angle) { - Bend::setBendAngle(angle); + Bend2D::setBendAngle(angle); setExitAngle(angle - getEntranceAngle()); } void RBend::setEntranceAngle(double entranceAngle) { - Bend::setEntranceAngle(entranceAngle); + Bend2D::setEntranceAngle(entranceAngle); setExitAngle(getBendAngle() - entranceAngle); } diff --git a/src/Classic/AbsBeamline/RBend.h b/src/Classic/AbsBeamline/RBend.h index 3965a86aa..f35200527 100644 --- a/src/Classic/AbsBeamline/RBend.h +++ b/src/Classic/AbsBeamline/RBend.h @@ -21,7 +21,7 @@ // // ------------------------------------------------------------------------ -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "BeamlineGeometry/RBendGeometry.h" #include "Fields/BMultipoleField.h" @@ -70,7 +70,7 @@ * Here we defined the magnet as a field map. */ -class RBend: public Bend { +class RBend: public Bend2D { public: @@ -183,4 +183,4 @@ private: double &chordLength) override; }; -#endif // CLASSIC_RBend_HH +#endif // CLASSIC_RBend_HH \ No newline at end of file diff --git a/src/Classic/AbsBeamline/RBend3D.h b/src/Classic/AbsBeamline/RBend3D.h index 96689edfa..02a69ac94 100644 --- a/src/Classic/AbsBeamline/RBend3D.h +++ b/src/Classic/AbsBeamline/RBend3D.h @@ -100,6 +100,7 @@ public: MeshData getSurfaceMesh() const; + virtual double getExitAngle() const; private: double trackRefParticleThrough(double dt, bool print = false); @@ -120,4 +121,9 @@ private: void operator=(const RBend3D &); }; -#endif // CLASSIC_RBend3D_HH +inline +double RBend3D::getExitAngle() const { + return getBendAngle() - getEntranceAngle(); +} + +#endif // CLASSIC_RBend3D_HH \ No newline at end of file diff --git a/src/Classic/AbsBeamline/SBend.cpp b/src/Classic/AbsBeamline/SBend.cpp index afb7577de..ba04adda1 100644 --- a/src/Classic/AbsBeamline/SBend.cpp +++ b/src/Classic/AbsBeamline/SBend.cpp @@ -33,19 +33,19 @@ extern Inform *gmsg; // ------------------------------------------------------------------------ SBend::SBend(): - Bend() + Bend2D() { setMessageHeader("SBend "); } SBend::SBend(const SBend &right): - Bend(right) + Bend2D(right) { setMessageHeader("SBend "); } SBend::SBend(const std::string &name): - Bend(name) + Bend2D(name) { setMessageHeader("SBend "); } diff --git a/src/Classic/AbsBeamline/SBend.h b/src/Classic/AbsBeamline/SBend.h index 9f0df6315..61e466997 100644 --- a/src/Classic/AbsBeamline/SBend.h +++ b/src/Classic/AbsBeamline/SBend.h @@ -21,7 +21,7 @@ // // ------------------------------------------------------------------------ -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "BeamlineGeometry/PlanarArcGeometry.h" #include "Fields/BMultipoleField.h" #include <string> @@ -65,7 +65,7 @@ * Here we defined the magnet as a field map. */ -class SBend: public Bend { +class SBend: public Bend2D { public: @@ -174,4 +174,4 @@ private: }; -#endif // CLASSIC_SBend_HH +#endif // CLASSIC_SBend_HH \ No newline at end of file diff --git a/src/Classic/Solvers/CSRIGFWakeFunction.cpp b/src/Classic/Solvers/CSRIGFWakeFunction.cpp index 474f08d84..bb15837df 100644 --- a/src/Classic/Solvers/CSRIGFWakeFunction.cpp +++ b/src/Classic/Solvers/CSRIGFWakeFunction.cpp @@ -108,27 +108,18 @@ void CSRIGFWakeFunction::apply(PartBunchBase<double, 3> *bunch) { } void CSRIGFWakeFunction::initialize(const ElementBase *ref) { - double End; if(ref->getType() == ElementBase::RBEND || ref->getType() == ElementBase::SBEND) { - const Bend *bend = static_cast<const Bend *>(ref); - // const RBend *bend = dynamic_cast<const RBend *>(ref); + const Bend2D *bend = static_cast<const Bend2D *>(ref); + double End; + bendRadius_m = bend->getBendRadius(); bend->getDimensions(Begin_m, End); Length_m = bend->getEffectiveLength(); FieldBegin_m = bend->getEffectiveCenter() - Length_m / 2.0; totalBendAngle_m = std::abs(bend->getBendAngle()); bendName_m = bend->getName(); - - // } else if(dynamic_cast<const SBend *>(ref)) { - // const SBend *bend = dynamic_cast<const SBend *>(ref); - // bendRadius_m = bend->getBendRadius(); - // bend->getDimensions(Begin_m, End); - // Length_m = bend->getEffectiveLength(); - // FieldBegin_m = bend->getEffectiveCenter() - Length_m / 2.0; - // totalBendAngle_m = bend->getBendAngle(); - // bendName_m = bend->getName(); } } diff --git a/src/Classic/Solvers/CSRWakeFunction.cpp b/src/Classic/Solvers/CSRWakeFunction.cpp index 2d792c55d..4c6c4e06e 100644 --- a/src/Classic/Solvers/CSRWakeFunction.cpp +++ b/src/Classic/Solvers/CSRWakeFunction.cpp @@ -113,10 +113,12 @@ void CSRWakeFunction::apply(PartBunchBase<double, 3> *bunch) { } void CSRWakeFunction::initialize(const ElementBase *ref) { - double End; if(ref->getType() == ElementBase::RBEND || ref->getType() == ElementBase::SBEND) { - const Bend *bend = static_cast<const Bend *>(ref); + + const Bend2D *bend = static_cast<const Bend2D *>(ref); + double End; + bendRadius_m = bend->getBendRadius(); bend->getDimensions(Begin_m, End); Length_m = bend->getEffectiveLength(); diff --git a/src/Classic/Structure/MeshGenerator.cpp b/src/Classic/Structure/MeshGenerator.cpp index 4248da385..bf70f5eea 100644 --- a/src/Classic/Structure/MeshGenerator.cpp +++ b/src/Classic/Structure/MeshGenerator.cpp @@ -1,7 +1,7 @@ #include "Structure/MeshGenerator.h" #include "Physics/Physics.h" #include "Utilities/Util.h" -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "AbsBeamline/RBend3D.h" #include "AbsBeamline/Multipole.h" @@ -22,7 +22,7 @@ void MeshGenerator::add(const ElementBase &element) { if (element.getType() == ElementBase::SBEND || element.getType() == ElementBase::RBEND) { - const Bend* dipole = static_cast<const Bend*>(&element); + const Bend2D* dipole = static_cast<const Bend2D*>(&element); mesh = dipole->getSurfaceMesh(); mesh.type_m = DIPOLE; } else if (element.getType() == ElementBase::RBEND3D) { diff --git a/src/Elements/OpalBeamline.cpp b/src/Elements/OpalBeamline.cpp index 20157a18b..216aceb9d 100644 --- a/src/Elements/OpalBeamline.cpp +++ b/src/Elements/OpalBeamline.cpp @@ -3,7 +3,7 @@ #include "Utilities/Options.h" #include "Utilities/Util.h" #include "AbstractObjects/OpalData.h" -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "Structure/MeshGenerator.h" #include <boost/filesystem.hpp> @@ -400,7 +400,7 @@ void OpalBeamline::plot3DLattice() { std::vector<Vector_t> corners; if (element->getType() == ElementBase::RBEND || element->getType() == ElementBase::SBEND) { - std::vector<Vector_t> outline = static_cast<const Bend*>(element.get())->getOutline(); + std::vector<Vector_t> outline = static_cast<const Bend2D*>(element.get())->getOutline(); for (auto point: outline) { corners.push_back(rotDiagonal.rotate(toBegin.transformFrom(point))); @@ -495,7 +495,7 @@ void OpalBeamline::save3DLattice() { if (element->getType() == ElementBase::SBEND || element->getType() == ElementBase::RBEND) { - Bend * bendElement = static_cast<Bend*>(element.get()); + Bend2D * bendElement = static_cast<Bend2D*>(element.get()); std::vector<Vector_t> designPath = bendElement->getDesignPath(); CoordinateSystemTrafo toEnd = bendElement->getBeginToEnd_local() * (*it).getCoordTransformationTo(); Vector_t exit3D = toEnd.getOrigin(); @@ -662,7 +662,7 @@ void OpalBeamline::save3DInput() { if ((*it).getElement()->getType() == ElementBase::RBEND || (*it).getElement()->getType() == ElementBase::SBEND) { - const Bend* dipole = static_cast<const Bend*>((*it).getElement().get()); + const Bend2D* dipole = static_cast<const Bend2D*>((*it).getElement().get()); double angle = dipole->getBendAngle(); double E1 = dipole->getEntranceAngle(); double E2 = dipole->getExitAngle(); @@ -699,7 +699,7 @@ void OpalBeamline::activateElements() { std::shared_ptr<Component> element = (*it).getElement(); if (element->getType() == ElementBase::SBEND || element->getType() == ElementBase::RBEND) { - Bend * bendElement = static_cast<Bend*>(element.get()); + Bend2D * bendElement = static_cast<Bend2D*>(element.get()); designEnergy = bendElement->getDesignEnergy() * 1e-6; } (*it).setOn(designEnergy); diff --git a/src/Elements/OpalElement.cpp b/src/Elements/OpalElement.cpp index 76b524d52..4a1a535c9 100644 --- a/src/Elements/OpalElement.cpp +++ b/src/Elements/OpalElement.cpp @@ -21,7 +21,7 @@ #include "Elements/OpalElement.h" #include "AbsBeamline/AlignWrapper.h" #include "AbsBeamline/ElementImage.h" -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "AbstractObjects/Attribute.h" #include "AbstractObjects/Expressions.h" #include "AbstractObjects/OpalData.h" diff --git a/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp b/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp index 47110b76f..268c9cf6d 100644 --- a/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp +++ b/tests/classic_src/AbsBeamline/DipoleFieldTest.cpp @@ -1,6 +1,6 @@ #include "gtest/gtest.h" #include <gtest/gtest_prod.h> -#include "AbsBeamline/Bend.h" +#include "AbsBeamline/Bend2D.h" #include "Fields/Fieldmap.h" #include "AbsBeamline/BeamlineVisitor.h" #include "BeamlineCore/RBendRep.h" @@ -102,7 +102,7 @@ TEST(Maxwell, Zeros) bunch->resetM(0.938); bunch->setdT(1.0e-12);//time step double startField = 2.0, endField = 10.0 ; - myMagnet->Bend::initialise(bunch, startField, endField); + myMagnet->Bend2D::initialise(bunch, startField, endField); delete partData; delete bunch; Vector_t R(0. ,0., 0.); @@ -120,9 +120,9 @@ TEST(Maxwell, Zeros) //std::cout<<"Step #"<<step<<endl; counter ++; Vector_t B(0.0); - R(0) = (myMagnet->Bend::designRadius_m + x) * cos(phi); + R(0) = (myMagnet->Bend2D::designRadius_m + x) * cos(phi); R(1) = z; - R(2) = (myMagnet->Bend::designRadius_m + x) * sin(phi); + R(2) = (myMagnet->Bend2D::designRadius_m + x) * sin(phi); double t = 0; myMagnet->apply(R, P, t , E, B); //B /= myMagnet->fieldAmplitude_m; //normalisation @@ -147,14 +147,14 @@ TEST(Maxwell, Zeros) } //fout.close(); - cout<<"bending radius: "<<myMagnet->Bend::designRadius_m<<endl; + cout<<"bending radius: "<<myMagnet->Bend2D::designRadius_m<<endl; cout<<"field amplitude: "<<myMagnet->fieldAmplitude_m<<endl; /** Vector_t B(0.0); double phi_new = Physics::pi / 20 ; - R[0] = (myMagnet->Bend::designRadius_m ) * cos(phi_new); + R[0] = (myMagnet->Bend2D::designRadius_m ) * cos(phi_new); R[1] = 0.005; - R[2] = (myMagnet->Bend::designRadius_m ) * sin(phi_new); + R[2] = (myMagnet->Bend2D::designRadius_m ) * sin(phi_new); myMagnet->apply(R, P, 0 , E, B); cout<<"Derivative vs stepSize:"<<endl; cout<<partialsDerivB(R, B,1.0, myMagnet)[1][2]<<' '<<1.0<<endl; @@ -206,4 +206,4 @@ TEST(Quad, Quadrupole) gout.close(); cout<<"length: "<<quad->getElementLength()<<endl; -} \ No newline at end of file +} -- GitLab