Skip to content

GitLab

Commit fba88103 authored by snuverink_j's avatar snuverink_j
Browse files

include particle charge into account for geometry of Bend2D, RBend3D; refactor code

parent 45486554
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 234 additions and 276 deletions
src/Classic/AbsBeamline/Bend2D.cpp
View file @ fba88103
This diff is collapsed.
src/Classic/AbsBeamline/Bend2D.h
View file @ fba88103
......@@ -33,10 +33,10 @@
#endif
#include <array>
#include <cmath>
#include <string>
#include <vector>
class Fieldmap;
class MeshData;
/*
......@@ -60,7 +60,7 @@ public:
virtual ~Bend2D();
/// Apply visitor to Bend2D.
virtual void accept(BeamlineVisitor &) const = 0;
virtual void accept(BeamlineVisitor &) const override = 0;
/*
* Methods for OPAL-T.
......@@ -71,51 +71,48 @@ public:
virtual bool apply(const size_t &i,
const double &t,
Vector_t &E,
Vector_t &B);
Vector_t &B) override;
/// Apply field to particles in beam frame.
virtual bool apply(const Vector_t &R,
const Vector_t &P,
const double &t,
Vector_t &E,
Vector_t &B);
Vector_t &B) override;
virtual bool applyToReferenceParticle(const Vector_t &R,
const Vector_t &P,
const double &t,
Vector_t &E,
Vector_t &B);
Vector_t &B) override;
virtual void goOnline(const double &kineticEnergy);
virtual void goOnline(const double &kineticEnergy) override;
virtual void finalise();
virtual void getDimensions(double &sBegin, double &sEnd) const;
virtual ElementBase::ElementType getType() const = 0;
virtual void finalise() override;
virtual void getDimensions(double &sBegin, double &sEnd) const override;
virtual ElementBase::ElementType getType() const override = 0;
virtual void initialise(PartBunchBase<double, 3> *bunch,
double &startField,
double &endField);
double &endField) override;
double getBendRadius() const;
double getEffectiveCenter() const;
double getEffectiveLength() const;
double getStartElement() const;
/// Set quadrupole field component.
void setK1(double k1);
virtual void setEntranceAngle(double entranceAngle) override;
void setExitAngle(double exitAngle);
virtual double getExitAngle() const;
double getMapLength() const;
virtual double getExitAngle() const override;
std::vector<Vector_t> getOutline() const;
MeshData getSurfaceMesh() const;
virtual CoordinateSystemTrafo getEdgeToEnd() const;
virtual CoordinateSystemTrafo getEdgeToEnd() const override;
CoordinateSystemTrafo getBeginToEnd_local() const;
virtual bool isInside(const Vector_t &r) const;
virtual bool isInside(const Vector_t &r) const override;
//set number of slices for map tracking
......@@ -135,7 +132,6 @@ public:
protected:
void setMessageHeader(const std::string & header);
double getStartField() const;
Fieldmap* getFieldmap();
private:
......@@ -164,18 +160,13 @@ private:
Vector_t &B);
void calculateRefTrajectory(double &angleX,
double &angleY);
double estimateFieldAdjustmentStep(double actualBendAngle,
double mass,
double betaGamma);
double estimateFieldAdjustmentStep(double actualBendAngle);
void findBendEffectiveLength(double startField,
double endField);
void findBendStrength(double mass,
double gamma,
double betaGamma,
double charge);
void findBendStrength();
virtual bool findChordLength(double &chordLength) = 0;
bool findIdealBendParameters(double chordLength);
bool initializeFieldMap(Inform &msg);
bool initializeFieldMap();
bool inMagnetCentralRegion(const Vector_t &R) const;
bool inMagnetEntranceRegion(const Vector_t &R) const;
bool inMagnetExitRegion(const Vector_t &R) const;
......@@ -188,11 +179,11 @@ private:
void setEngeOriginDelta(double delta);
void setFieldCalcParam();
void setGapFromFieldMap();
bool setupBendGeometry(Inform &msg, double &startField, double &endField);
bool setupDefaultFieldMap(Inform &msg);
bool setupBendGeometry(double &startField, double &endField);
bool setupDefaultFieldMap();
void setFieldBoundaries(double startField, double endField);
void setupPusher(PartBunchBase<double, 3> *bunch);
bool treatAsDrift(Inform &msg, double chordlength);
bool isFieldZero();
void setCSTrafoToEntranceRegion(const CoordinateSystemTrafo &trafo);
void setCSTrafoToExitRegion(const CoordinateSystemTrafo &trafo);
......@@ -204,12 +195,8 @@ private:
BorisPusher pusher_m; /// Pusher used to integrate reference particle
/// through the bend.
Fieldmap *fieldmap_m; /// Magnet field map.
const bool fast_m = false; /// Flag to turn on fast field calculation.
double designRadius_m; /// Bend design radius (m).
/// and the entrance face of the magnet (radians).
double exitAngle_m; /// Angle between outgoing reference trajectory
/// and the exit face of the magnet (radians).
double fieldIndex_m; /// Dipole field index.
......@@ -313,11 +300,6 @@ double Bend2D::getEffectiveLength() const {
return designRadius_m * angle_m;
}
inline
double Bend2D::getStartElement() const {
return elementEdge_m;
}
inline
void Bend2D::setK1(double k1) {
if (std::abs(k1) > 0.0) {
......@@ -339,12 +321,6 @@ double Bend2D::getStartField() const
return startField_m;
}
inline
Fieldmap* Bend2D::getFieldmap()
{
return fieldmap_m;
}
inline
double Bend2D::getExitAngle() const
{
......@@ -352,16 +328,19 @@ double Bend2D::getExitAngle() const
}
inline
void Bend2D::setExitAngle(double angle)
void Bend2D::setEntranceAngle(double angle)
{
exitAngle_m = angle;
tanExitAngle_m = tan(exitAngle_m);
BendBase::setEntranceAngle(angle);
cosEntranceAngle_m = std::cos(entranceAngle_m);
sinEntranceAngle_m = std::sin(entranceAngle_m);
tanEntranceAngle_m = std::tan(entranceAngle_m);
}
inline
double Bend2D::getMapLength() const
void Bend2D::setExitAngle(double angle)
{
return exitParameter2_m - entranceParameter2_m;
exitAngle_m = angle;
tanExitAngle_m = std::tan(exitAngle_m);
}
inline
......
src/Classic/AbsBeamline/BendBase.cpp
View file @ fba88103
#include "AbsBeamline/BendBase.h"
#include "Algorithms/PartBunchBase.h"
#include <cmath>
BendBase::BendBase():
BendBase("")
{}
......@@ -10,12 +14,13 @@ BendBase::BendBase(const BendBase &right):
chordLength_m(right.chordLength_m),
angle_m(right.angle_m),
entranceAngle_m(right.entranceAngle_m),
fieldmap_m(right.fieldmap_m),
gap_m(right.gap_m),
designEnergy_m(right.designEnergy_m),
designEnergyChangeable_m(true),
refTrajMap_m(right.refTrajMap_m),
bX_m(right.bX_m),
bY_m(right.bY_m),
fieldAmplitudeX_m(right.fieldAmplitudeX_m),
fieldAmplitudeY_m(right.fieldAmplitudeY_m),
fieldAmplitude_m(right.fieldAmplitude_m),
fileName_m(right.fileName_m)
{}
......@@ -26,11 +31,12 @@ BendBase::BendBase(const std::string &name):
chordLength_m(0.0),
angle_m(0.0),
entranceAngle_m(0.0),
fieldmap_m(NULL),
gap_m(0.0),
designEnergy_m(0.0),
designEnergyChangeable_m(true),
bX_m(0.0),
bY_m(0.0),
fieldAmplitudeX_m(0.0),
fieldAmplitudeY_m(0.0),
fieldAmplitude_m(0.0),
fileName_m("")
{}
......@@ -47,4 +53,49 @@ std::vector<Vector_t> BendBase::getDesignPath() const {
}
return designPath;
}
void BendBase::setFieldAmplitude(double k0, double k0s) {
fieldAmplitudeY_m = k0;
fieldAmplitudeX_m = k0s;
}
double BendBase::calcDesignRadius(double fieldAmplitude) const
{
double mass = RefPartBunch_m->getM();
double betaGamma = calcBetaGamma();
double charge = RefPartBunch_m->getQ();
// Lorentz force: condition for a circular orbit
return std::abs(betaGamma * mass / (Physics::c * fieldAmplitude * charge));
}
double BendBase::calcFieldAmplitude(double radius) const
{
double mass = RefPartBunch_m->getM();
double betaGamma = calcBetaGamma();
double charge = RefPartBunch_m->getQ();
// Lorentz force: condition for a circular orbit
return betaGamma * mass / (Physics::c * radius * charge);
}
double BendBase::calcBendAngle(double chordLength, double radius) const
{
return 2.0 * std::asin(chordLength / (2.0 * radius));
}
double BendBase::calcDesignRadius(double chordLength, double angle) const
{
return chordLength / (2.0 * std::sin(angle / 2.0));
}
double BendBase::calcGamma() const
{
double mass = RefPartBunch_m->getM();
return designEnergy_m / mass + 1.0;
}
double BendBase::calcBetaGamma() const
{
double gamma = calcGamma();
return std::sqrt(std::pow(gamma, 2.0) - 1.0);
}
\ No newline at end of file
src/Classic/AbsBeamline/BendBase.h
View file @ fba88103
......@@ -6,6 +6,8 @@
#include <vector>
#include <string>
class Fieldmap;
class BendBase: public Component {
public:
BendBase();
......@@ -37,21 +39,40 @@ public:
std::string getFieldMapFN() const;
protected:
/// Calculate design radius from design energy and field amplitude
double calcDesignRadius(double fieldAmplitude) const;
/// Calculate field amplitude from design energy and radius
double calcFieldAmplitude(double radius) const;
/// Calculate bend angle from chord length and design radius
double calcBendAngle(double chordLength, double radius) const;
/// Calculate design radius from chord length and bend angle
double calcDesignRadius(double chordLength, double angle) const;
/// Calculate gamma from design energy
double calcGamma() const;
/// Calculate beta*gamma from design energy
double calcBetaGamma() const;
double length_m;
double chordLength_m;
double angle_m;
double entranceAngle_m; /// Angle between incoming reference trajectory
double angle_m; ///< Bend angle
double entranceAngle_m; ///< Angle between incoming reference trajectory
///< and the entrance face of the magnet (radians).
Fieldmap *fieldmap_m; ///< Magnet field map.
const bool fast_m = false; ///< Flag to turn on fast field calculation.
double gap_m;
double gap_m; ///< Full vertical gap of the magnets.
double designEnergy_m; /// Bend design energy (eV).
double designEnergy_m; ///< Bend design energy (eV).
bool designEnergyChangeable_m;
/// Map of reference particle trajectory.
std::vector<Vector_t> refTrajMap_m;
double bX_m;
double bY_m;
double fieldAmplitude_m;
double fieldAmplitudeX_m; ///< Field amplitude in x direction.
///< Value not updated if user defines strength with angle
double fieldAmplitudeY_m; ///< Field amplitude in y direction.
///< Value not updated if user defines strength with angle
double fieldAmplitude_m; ///< Field amplitude.
std::string fileName_m;
};
......@@ -121,12 +142,6 @@ double BendBase::getDesignEnergy() const {
return designEnergy_m;
}
inline
void BendBase::setFieldAmplitude(double k0, double k0s) {
bY_m = k0;
bX_m = k0s;
}
inline
double BendBase::getFieldAmplitude() const
{
......
src/Classic/AbsBeamline/RBend.cpp
View file @ fba88103
......@@ -138,13 +138,10 @@ bool RBend::findChordLength(double &chordLength) {
chordLength = 2 * getLength() * std::sin(0.5 * std::abs(angle)) /
(std::sin(E1) + std::sin(std::abs(angle) - E1));
} else {
double refMass = RefPartBunch_m->getM();
double refGamma = designEnergy_m / refMass + 1.0;
double refBetaGamma = std::sqrt(std::pow(refGamma, 2.0) - 1.0);
double refCharge = RefPartBunch_m->getQ();
double amplitude = (std::abs(bY_m) > 0.0? bY_m: bX_m);
double fieldAmplitude = refCharge * std::abs(amplitude / refCharge);
double designRadius = std::abs(refBetaGamma * refMass / (Physics::c * fieldAmplitude));
double amplitude = (fieldAmplitudeY_m != 0.0) ? fieldAmplitudeY_m : fieldAmplitudeX_m;
double fieldAmplitude = std::copysign(1.0, refCharge) * std::abs(amplitude);
double designRadius = calcDesignRadius(fieldAmplitude);
chordLength = std::sin(0.5 * (std::asin(getLength() / designRadius - std::sin(getEntranceAngle())) + getEntranceAngle())) * 2 * designRadius;
}
......
src/Classic/AbsBeamline/RBend3D.cpp
View file @ fba88103
......@@ -39,11 +39,9 @@ RBend3D::RBend3D():
RBend3D::RBend3D(const RBend3D &right):
BendBase(right),
myFieldmap_m(right.myFieldmap_m),
fieldAmplitudeError_m(right.fieldAmplitudeError_m),
startField_m(right.startField_m),
lengthField_m(right.lengthField_m),
fast_m(right.fast_m),
geometry_m(right.geometry_m),
dummyField_m() {
}
......@@ -51,37 +49,21 @@ RBend3D::RBend3D(const RBend3D &right):
RBend3D::RBend3D(const std::string &name):
BendBase(name),
myFieldmap_m(NULL),
fieldAmplitudeError_m(0.0),
startField_m(0.0),
lengthField_m(0.0),
fast_m(false),
geometry_m(),
dummyField_m() {
}
RBend3D::~RBend3D() {
// Fieldmap::deleteFieldmap(filename_m);
}
void RBend3D::accept(BeamlineVisitor &visitor) const {
visitor.visitRBend3D(*this);
}
void RBend3D::setFieldMapFN(std::string fn) {
fileName_m = fn;
}
void RBend3D::setFast(bool fast) {
fast_m = fast;
}
bool RBend3D::getFast() const {
return fast_m;
}
bool RBend3D::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return apply(RefPartBunch_m->R[i], RefPartBunch_m->P[i], t, E, B);
}
......@@ -90,7 +72,7 @@ bool RBend3D::apply(const Vector_t &R, const Vector_t &/*P*/, const double &/*t
const Vector_t tmpR(R(0), R(1), R(2) - startField_m);
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
myFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
B += (fieldAmplitude_m + fieldAmplitudeError_m) * tmpB;
......@@ -101,7 +83,7 @@ bool RBend3D::applyToReferenceParticle(const Vector_t &R, const Vector_t &/*P*/,
const Vector_t tmpR(R(0), R(1), R(2) - startField_m);
Vector_t tmpE(0.0), tmpB(0.0);
myFieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
fieldmap_m->getFieldstrength(tmpR, tmpE, tmpB);
B += fieldAmplitude_m * tmpB;
return false;
......@@ -112,14 +94,14 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
RefPartBunch_m = bunch;
myFieldmap_m = Fieldmap::getFieldmap(fileName_m, fast_m);
if(myFieldmap_m != NULL) {
fieldmap_m = Fieldmap::getFieldmap(fileName_m, fast_m);
if(fieldmap_m != NULL) {
msg << level2 << getName() << " using file ";
myFieldmap_m->getInfo(&msg);
fieldmap_m->getInfo(&msg);
goOnline(0.0);
double zBegin = 0.0, zEnd = 0.0, rBegin = 0.0, rEnd = 0.0;
myFieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);
fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);
if (length_m == 0.0) {
chordLength_m = 0.0;
......@@ -127,14 +109,14 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
double z = 0.0, dz = fieldLength / 1000;
Vector_t E(0.0), B(0.0);
while (z < fieldLength && B(1) < 0.5) {
myFieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
fieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
z += dz;
}
double zEntryEdge = z;
z = fieldLength;
B(1) = 0.0;
while (z > 0.0 && B(1) < 0.5) {
myFieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
fieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
z -= dz;
}
chordLength_m = z - zEntryEdge;
......@@ -147,16 +129,7 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
lengthField_m = zEnd - zBegin;
endField = startField + lengthField_m;
if (bX_m * bX_m + bY_m * bY_m > 0.0) {
double refCharge = bunch->getQ();
rotationZAxis_m += std::atan2(bX_m, bY_m);
if (refCharge < 0.0) {
rotationZAxis_m -= Physics::pi;
}
fieldAmplitude_m = (refCharge *
std::abs(std::sqrt(std::pow(bY_m, 2.0) + std::pow(bX_m, 2.0)) / refCharge));
angle_m = trackRefParticleThrough(bunch->getdT(), Options::writeBendTrajectories);
} else {
if (angle_m != 0.0) {
if (angle_m < 0.0) {
// Negative angle is a positive bend rotated 180 degrees.
entranceAngle_m = std::copysign(1, angle_m) * entranceAngle_m;
......@@ -164,21 +137,16 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
rotationZAxis_m += Physics::pi;
}
const double refCharge = RefPartBunch_m->getQ();
const double refMass = RefPartBunch_m->getM();
const double refGamma = designEnergy_m / refMass + 1.0;
const double refBetaGamma = std::sqrt(std::pow(refGamma, 2.0) - 1.0);
Vector_t B(0.0), E(0.0);
double z = 0.0, dz = lengthField_m / 999;
double integratedBy = 0.0;
myFieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
fieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
integratedBy += 0.5 * B(1);
z = dz;
while (z < lengthField_m) {
B = 0.0; E = 0.0;
myFieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
fieldmap_m->getFieldstrength(Vector_t(0.0, 0.0, z), E, B);
integratedBy += B(1);
z += dz;
}
......@@ -186,7 +154,7 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
integratedBy *= lengthField_m / 1000;
// estimate magnitude of field
fieldAmplitude_m = refCharge * refMass * refBetaGamma * angle_m / (Physics::c * integratedBy);
fieldAmplitude_m = calcFieldAmplitude(integratedBy / angle_m);
double angle = trackRefParticleThrough(bunch->getdT());
double error = (angle_m - angle) / angle_m;
......@@ -201,7 +169,16 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
if (Options::writeBendTrajectories)
trackRefParticleThrough(bunch->getdT(), true);
} else {
double refCharge = bunch->getQ();
rotationZAxis_m += std::atan2(fieldAmplitudeX_m, fieldAmplitudeY_m);
if (refCharge < 0.0) {
rotationZAxis_m -= Physics::pi;
}
fieldAmplitude_m = std::copysign(1.0, refCharge) * std::hypot(fieldAmplitudeX_m, fieldAmplitudeY_m);
angle_m = trackRefParticleThrough(bunch->getdT(), Options::writeBendTrajectories);
}
} else {
endField = startField;
}
......@@ -210,10 +187,6 @@ void RBend3D::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
void RBend3D::finalise()
{}
bool RBend3D::bends() const {
return true;
}
void RBend3D::goOnline(const double &) {
Fieldmap::readMap(fileName_m);
online_m = true;
......@@ -236,7 +209,7 @@ ElementBase::ElementType RBend3D::getType() const {
bool RBend3D::isInside(const Vector_t &r) const {
Vector_t tmpR(r(0), r(1), r(2) - startField_m);
return myFieldmap_m->isInside(tmpR);
return fieldmap_m->isInside(tmpR);
}
ElementBase* RBend3D::clone() const {
......@@ -261,9 +234,9 @@ const EMField &RBend3D::getField() const {
MeshData RBend3D::getSurfaceMesh() const {
Vector_t XIni, XFinal;
myFieldmap_m->getFieldDimensions(XIni(0), XFinal(0),
XIni(1), XFinal(1),
XIni(2), XFinal(2));
fieldmap_m->getFieldDimensions(XIni(0), XFinal(0),
XIni(1), XFinal(1),
XIni(2), XFinal(2));
MeshData mesh;
mesh.vertices_m.push_back(Vector_t(XIni(0), XIni(1), XIni(2)));
......@@ -292,10 +265,9 @@ MeshData RBend3D::getSurfaceMesh() const {
}
double RBend3D::trackRefParticleThrough(double dt, bool print) {
const double refMass = RefPartBunch_m->getM();
const double refGamma = designEnergy_m / refMass + 1.0;
const double refBetaGamma = std::sqrt(std::pow(refGamma, 2.0) - 1.0);
const double stepSize = refBetaGamma / refGamma * Physics::c * dt;
const double refGamma = calcGamma();
const double refBetaGamma = calcBetaGamma();
const double stepSize = refBetaGamma / refGamma * Physics::c * dt;
const Vector_t scaleFactor(Physics::c * dt);
print = print && (Ippl::myNode() == 0);
......
src/Classic/AbsBeamline/RBend3D.h
View file @ fba88103
......@@ -75,19 +75,10 @@ public:
virtual void finalise() override;
virtual bool bends() const override;
virtual void goOnline(const double &kineticEnergy) override;
virtual void goOffline() override;
// Assign the field filename.
void setFieldMapFN(std::string fn);
void setFast(bool fast);
bool getFast() const;
virtual ElementBase::ElementType getType() const override;
virtual void getDimensions(double &zBegin, double &zEnd) const override;
......@@ -100,15 +91,11 @@ public:
private:
double trackRefParticleThrough(double dt, bool print = false);
// std::string name; /**< The name of the object*/
Fieldmap *myFieldmap_m;
double fieldAmplitudeError_m; /**< scale multiplier error*/
double startField_m; /**< startingpoint of field, m*/