Commit ae23e85d authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Merge branch '402-vffa' of gitlab.psi.ch:OPAL/src into origin/402-vffa -...

Merge branch '402-vffa' of gitlab.psi.ch:OPAL/src into origin/402-vffa - implement VFFA analytical field model.
parent 37a2ca48
......@@ -61,6 +61,7 @@
#include "AbsBeamline/Stripper.h"
#include "AbsBeamline/VariableRFCavity.h"
#include "AbsBeamline/VariableRFCavityFringeField.h"
#include "AbsBeamline/VerticalFFAMagnet.h"
#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"
......@@ -893,6 +894,15 @@ void ParallelCyclotronTracker::visitVariableRFCavityFringeField
"Need to define a RINGDEFINITION to use VariableRFCavity element");
}
void ParallelCyclotronTracker::visitVerticalFFAMagnet(const VerticalFFAMagnet &mag) {
*gmsg << "Adding Vertical FFA Magnet" << endl;
if (opalRing_m != NULL)
opalRing_m->appendElement(mag);
else
throw OpalException("ParallelCyclotronTracker::visitVerticalFFAMagnet",
"Need to define a RINGDEFINITION to use VerticalFFAMagnet element");
}
/**
*
*
......
......@@ -181,6 +181,9 @@ public:
virtual void visitVariableRFCavityFringeField
(const VariableRFCavityFringeField &cav);
/// Apply the algorithm to a VerticalFFAMagnet.
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &bend);
/// Apply the algorithm to the top-level beamline.
// overwrite the execute-methode from DefaultVisitor
virtual void execute();
......
......@@ -73,6 +73,7 @@ class Stripper;
class TravelingWave;
class VariableRFCavity;
class VariableRFCavityFringeField;
class VerticalFFAMagnet;
// Integrators.
class Integrator;
......@@ -261,6 +262,9 @@ public:
/// Apply the algorithm to an integrator capable of mapping.
virtual void visitMapIntegrator(const MapIntegrator &) = 0;
/// Apply the algorithm to a vertical FFA magnet
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &) = 0;
private:
// Not implemented.
......
......@@ -51,6 +51,7 @@ set (_SRCS
TravelingWave.cpp
VariableRFCavity.cpp
VariableRFCavityFringeField.cpp
VerticalFFAMagnet.cpp
)
include_directories (
......@@ -111,6 +112,7 @@ set (HDRS
TravelingWave.h
VariableRFCavity.h
VariableRFCavityFringeField.h
VerticalFFAMagnet.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
......
......@@ -120,6 +120,23 @@ public:
Vector_t &E,
Vector_t &B);
/** Calculate the four-potential at some position relative to the component
*
* \param R position in the local coordinate system of the component
* \param t time
* \param A filled with the calculated magnetic vector potential
* \param phi filled with the calculated electric potential
* Note that any existing values in A and phi may be overwritten by this
* method.
*
* \returns true if particle is outside the field map, else false
* Default for component is to return false and make no change to A and phi
*/
virtual bool getPotential(const Vector_t &R,
const double &t,
Vector_t &A,
double &phi) {return false;}
virtual double getDesignEnergy() const;
virtual void setDesignEnergy(const double& energy, bool changeable);
......
......@@ -37,6 +37,7 @@
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/SBend3D.h"
#include "AbsBeamline/ScalingFFAMagnet.h"
#include "AbsBeamline/VerticalFFAMagnet.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
......@@ -190,6 +191,9 @@ public:
/// Apply the algorithm to a spiral sector magnet.
virtual void visitScalingFFAMagnet(const ScalingFFAMagnet &);
/// Apply the algorithm to a spiral sector magnet.
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &);
/// Apply the algorithm to a ParallelPlate.
virtual void visitParallelPlate(const ParallelPlate &);
......@@ -431,6 +435,11 @@ void SpecificElementVisitor<ELEM>::visitScalingFFAMagnet(const ScalingFFAMagnet
CastsTrait<ELEM, ScalingFFAMagnet>::apply(allElementsOfTypeE, element);
}
template<class ELEM>
void SpecificElementVisitor<ELEM>::visitVerticalFFAMagnet(const VerticalFFAMagnet &element) {
CastsTrait<ELEM, VerticalFFAMagnet>::apply(allElementsOfTypeE, element);
}
template<class ELEM>
void SpecificElementVisitor<ELEM>::visitSeparator(const Separator &element) {
CastsTrait<ELEM, Separator>::apply(allElementsOfTypeE, element);
......
//
// Source file for VerticalFFAMagnet Component
//
// Copyright (c) 2019 Chris Rogers
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
//
#include "AbsBeamline/BeamlineVisitor.h"
#include "AbsBeamline/EndFieldModel/EndFieldModel.h"
#include "AbsBeamline/VerticalFFAMagnet.h"
VerticalFFAMagnet::VerticalFFAMagnet(const std::string &name)
: Component(name), straightGeometry_m(1.) {
setElType(isDrift);
}
VerticalFFAMagnet::VerticalFFAMagnet(const VerticalFFAMagnet &right)
: Component(right),
straightGeometry_m(right.straightGeometry_m),
dummy(right.dummy),
maxOrder_m(right.maxOrder_m),
k_m(right.k_m),
Bz_m(right.Bz_m),
zNegExtent_m(right.zNegExtent_m),
zPosExtent_m(right.zPosExtent_m),
halfWidth_m(right.halfWidth_m),
bbLength_m(right.bbLength_m),
endField_m(right.endField_m->clone()),
dfCoefficients_m(right.dfCoefficients_m) {
RefPartBunch_m = right.RefPartBunch_m;
}
VerticalFFAMagnet::~VerticalFFAMagnet() {
}
ElementBase* VerticalFFAMagnet::clone() const {
VerticalFFAMagnet* magnet = new VerticalFFAMagnet(*this);
magnet->initialise();
return magnet;
}
EMField &VerticalFFAMagnet::getField() {
return dummy;
}
const EMField &VerticalFFAMagnet::getField() const {
return dummy;
}
void VerticalFFAMagnet::initialise() {
calculateDfCoefficients();
straightGeometry_m.setElementLength(bbLength_m); // length = phi r
}
void VerticalFFAMagnet::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
RefPartBunch_m = bunch;
initialise();
}
void VerticalFFAMagnet::finalise() {
RefPartBunch_m = NULL;
}
BGeometryBase& VerticalFFAMagnet::getGeometry() {
return straightGeometry_m;
}
const BGeometryBase& VerticalFFAMagnet::getGeometry() const {
return straightGeometry_m;
}
void VerticalFFAMagnet::accept(BeamlineVisitor& visitor) const {
visitor.visitVerticalFFAMagnet(*this);
}
bool VerticalFFAMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
if (abs(R[0]) > halfWidth_m ||
R[2] < 0. || R[2] > bbLength_m ||
R[1] < -zNegExtent_m || R[1] > zPosExtent_m) {
return true;
}
std::vector<double> fringeDerivatives(maxOrder_m+2, 0.);
double zRel = R[2]-bbLength_m/2.; // z relative to centre of magnet
for (size_t i = 0; i < fringeDerivatives.size(); ++i) {
fringeDerivatives[i] = endField_m->function(zRel, i); // d^i_phi f
}
std::vector<double> x_n(maxOrder_m+1); // x^n
x_n[0] = 1.; // x^0
for (size_t i = 1; i < x_n.size(); ++i) {
x_n[i] = x_n[i-1]*R[0];
}
// note that the last element is always 0, because dfCoefficients_m is
// of size maxOrder_m+1. This leads to better Maxwellianness in testing.
std::vector<double> f_n(maxOrder_m+2, 0.);
std::vector<double> dz_f_n(maxOrder_m+1, 0.);
for (size_t n = 0; n < dfCoefficients_m.size(); ++n) {
const std::vector<double>& coefficients = dfCoefficients_m[n];
for (size_t i = 0; i < coefficients.size(); ++i) {
f_n[n] += coefficients[i]*fringeDerivatives[i];
dz_f_n[n] += coefficients[i]*fringeDerivatives[i+1];
}
}
double bref = Bz_m*exp(k_m*R[1]);
B[0] = 0.;
B[1] = 0.;
B[2] = 0.;
for (size_t n = 0; n < x_n.size(); ++n) {
B[0] += bref*f_n[n+1]*(n+1)/k_m*x_n[n];
B[1] += bref*f_n[n]*x_n[n];
B[2] += bref*dz_f_n[n]/k_m*x_n[n];
}
return false;
}
void VerticalFFAMagnet::calculateDfCoefficients() {
dfCoefficients_m = std::vector< std::vector<double> >(maxOrder_m+1);
dfCoefficients_m[0] = std::vector<double>(1, 1.);
if (maxOrder_m > 0) {
dfCoefficients_m[1] = std::vector<double>();
}
// n indexes like the polynomial order of the midplane expansion
// e.g. Bz = exp(mz) f_n y^n
// where y is distance from the midplane and z is height
for (size_t n = 2; n < dfCoefficients_m.size(); n+=2) {
const std::vector<double>& oldCoefficients = dfCoefficients_m[n-2];
std::vector<double> coefficients(oldCoefficients.size()+2, 0);
// j indexes the derivative of f_0
for (size_t j = 0; j < oldCoefficients.size(); ++j) {
coefficients[j] += -1./(n)/(n-1)*k_m*k_m*oldCoefficients[j];
coefficients[j+2] += -1./(n)/(n-1)*oldCoefficients[j];
}
dfCoefficients_m[n] = coefficients;
}
}
void VerticalFFAMagnet::setEndField(endfieldmodel::EndFieldModel* endField) {
endField_m.reset(endField);
}
//
// Header file for VerticalFFAMagnet Component
//
// Copyright (c) 2019 Chris Rogers
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
//
#include "Fields/BMultipoleField.h"
#include "BeamlineGeometry/StraightGeometry.h"
#include "AbsBeamline/Component.h"
#include "Algorithms/PartBunch.h"
#ifndef ABSBEAMLINE_VerticalFFAMagnet_H
#define ABSBEAMLINE_VerticalFFAMagnet_H
namespace endfieldmodel {
class EndFieldModel;
}
/** Bending magnet with an exponential dependence on field in the vertical plane
*
* VerticalFFAMagnet makes a rectangular bending magnet with a dipole field
* that has a dependence like B0 exp(mz)
*/
class VerticalFFAMagnet : public Component {
public:
/** Construct a new VerticalFFAMagnet
*
* \param name User-defined name of the VerticalFFAMagnet
*/
explicit VerticalFFAMagnet(const std::string &name);
/** Destructor - deletes the field */
~VerticalFFAMagnet();
/** Inheritable copy constructor */
ElementBase* clone() const;
/** Calculate the field at the position of the ith particle
*
* \param i index of the particle event; field is calculated at this
* position
* \param t time at which the field is to be calculated
* \param E calculated electric field - always 0 (no E-field)
* \param B calculated magnetic field
* \returns true if particle is outside the field map
*/
inline bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
/** Calculate the field at some arbitrary position
*
* \param R position in the local coordinate system of the magnet
* \param P not used
* \param t not used
* \param E not used
* \param B calculated magnetic field
* \returns true if particle is outside the field map, else false
*/
inline bool apply(const Vector_t &R, const Vector_t &P, const double &t,
Vector_t &E, Vector_t &B);
/** Calculate the field at some arbitrary position in cartesian coordinates
*
* \param R position in the local coordinate system of the bend, in
* cartesian coordinates defined like (x, y, z)
* \param B calculated magnetic field defined like (Bx, By, Bz)
* \returns true if particle is outside the field map, else false
*/
bool getFieldValue(const Vector_t &R, Vector_t &B) const;
/** Initialise the VerticalFFAMagnet
*
* \param bunch the global bunch object (but not used)
* \param startField not used
* \param endField not used
*/
void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField);
/** Initialise the VerticalFFAMagnet
*
* Sets up the field expansion and the geometry; call after changing any
* field parameters
*/
void initialise();
/** Finalise the VerticalFFAMagnet - sets bunch to NULL */
void finalise();
/** Return false - VerticalFFAMagnet is a straight magnet
*
* Nb: the VerticalFFAMagnet geometry is straight even though trajectories
* are not
*/
inline bool bends() const {return false;}
/** Not implemented */
void getDimensions(double &zBegin, double &zEnd) const {}
/** Return the cell geometry */
BGeometryBase& getGeometry();
/** Return the cell geometry */
const BGeometryBase& getGeometry() const;
/** Return a dummy (0.) field value (what is this for?) */
EMField &getField();
/** Return a dummy (0.) field value (what is this for?) */
const EMField &getField() const;
/** Accept a beamline visitor */
void accept(BeamlineVisitor& visitor) const;
/** Get the fringe field
*
* Returns the fringe field model; VerticalFFAMagnet retains ownership of
* the returned memory.
*/
endfieldmodel::EndFieldModel* getEndField() const {return endField_m.get();}
/** Set the fringe field
*
* - endField: the new fringe field; VerticalFFAMagnet takes ownership of
* the memory associated with endField.
*/
void setEndField(endfieldmodel::EndFieldModel* endField);
/** Get the maximum power of x used in the off-midplane expansion;
*/
size_t getMaxOrder() const {return maxOrder_m;}
/** Set the maximum power of x used in the off-midplane expansion;
*/
void setMaxOrder(size_t maxOrder) {maxOrder_m = maxOrder;}
/** Get the centre field at z=0 */
double getB0() const {return Bz_m/Tesla;}
/** Set the centre field at z=0 */
void setB0(double Bz) {Bz_m = Bz*Tesla;}
/** Get the field index */
double getFieldIndex() const {return k_m*mm;} // units are [m^{-1}]
/** Set the field index */
void setFieldIndex(double index) {k_m = index/mm;}
/** Get the maximum extent below z = 0 */
double getNegativeVerticalExtent() const {return zNegExtent_m/mm;}
/** Set the maximum extent below z = 0 */
inline void setNegativeVerticalExtent(double negativeExtent);
/** Get the maximum extent above z = 0 */
double getPositiveVerticalExtent() const {return zPosExtent_m/mm;}
/** set the maximum extent above z = 0 */
inline void setPositiveVerticalExtent(double positiveExtent);
/** Get the length of the bounding box (centred on magnet centre) */
double getBBLength() const {return bbLength_m/mm;}
/** Set the length of the bounding box (centred on magnet centre) */
void setBBLength(double bbLength) {bbLength_m = bbLength*mm;}
/** Get the full width of the bounding box (centred on magnet centre) */
double getWidth() const {return halfWidth_m/mm*2.;}
/** Set the full width of the bounding box (centred on magnet centre) */
void setWidth(double width) {halfWidth_m = width/2*mm;}
/** Get the coefficients used for the field expansion
*
* B_y is given by
* sum_n B_0 exp(ky) f_n x^n
* where
* f_n = sum_k c_{nk} partial_k f_0
*
* Returns a vector of vectors, like c[n][k]. The expansion for the other
* field elements can be related back to c[n][k] (see elsewhere for details).
*/
inline std::vector<std::vector<double> > getDfCoefficients() const;
private:
void calculateDfCoefficients();
/** Copy constructor */
VerticalFFAMagnet(const VerticalFFAMagnet &right);
VerticalFFAMagnet& operator=(const VerticalFFAMagnet& rhs);
StraightGeometry straightGeometry_m;
BMultipoleField dummy;
size_t maxOrder_m = 0;
double k_m = 0.;
double Bz_m = 0.;
double zNegExtent_m = 0.; // extent downwards from the midplane
double zPosExtent_m = 0.; // extent upwards from the midplane
double halfWidth_m = 0.; // extent in either +x or -x
double bbLength_m = 0.;
std::unique_ptr<endfieldmodel::EndFieldModel> endField_m;
std::vector<std::vector<double> > dfCoefficients_m;
const double mm=1000.;
const double Tesla=10.;
};
void VerticalFFAMagnet::setNegativeVerticalExtent(double negativeExtent) {
zNegExtent_m = negativeExtent*mm;
}
void VerticalFFAMagnet::setPositiveVerticalExtent(double positiveExtent) {
zPosExtent_m = positiveExtent*mm;
}
bool VerticalFFAMagnet::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);
}
bool VerticalFFAMagnet::apply(const Vector_t &R, const Vector_t &P,
const double &t, Vector_t &E, Vector_t &B) {
return getFieldValue(R, B);
}
std::vector<std::vector<double> > VerticalFFAMagnet::getDfCoefficients() const {
return dfCoefficients_m;
}
#endif
......@@ -52,6 +52,7 @@
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/SBend3D.h"
#include "AbsBeamline/ScalingFFAMagnet.h"
#include "AbsBeamline/VerticalFFAMagnet.h"
#include "AbsBeamline/Separator.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
......@@ -232,6 +233,9 @@ void DefaultVisitor::visitScalingFFAMagnet(const ScalingFFAMagnet &spiral) {
applyDefault(spiral);
}
void DefaultVisitor::visitVerticalFFAMagnet(const VerticalFFAMagnet &mag) {
applyDefault(mag);
}
void DefaultVisitor::visitSeparator(const Separator &sep) {
applyDefault(sep);
......
......@@ -148,6 +148,9 @@ public:
/// Apply the algorithm to a scaling FFA magnet.
virtual void visitScalingFFAMagnet(const ScalingFFAMagnet &);
/// Apply the algorithm to a RF cavity.
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &);
/// Apply the algorithm to a separator.
virtual void visitSeparator(const Separator &);
......
......@@ -49,6 +49,7 @@ set (_SRCS
OpalSRot.cpp
OpalVariableRFCavity.cpp
OpalVariableRFCavityFringeField.cpp
OpalVerticalFFAMagnet.cpp
OpalVKicker.cpp
OpalVMonitor.cpp
OpalYRot.cpp
......@@ -114,6 +115,7 @@ set (HDRS
OpalTravelingWave.h
OpalVariableRFCavity.h
OpalVariableRFCavityFringeField.h
OpalVerticalFFAMagnet.h
OpalVKicker.h
OpalVMonitor.h
OpalYRot.h
......
//
// Source file for OpalVerticalFFAMagnet Element
//
// Copyright (c) 2019 Chris Rogers
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
#include "Attributes/Attributes.h"
#include "AbsBeamline/EndFieldModel/Tanh.h" // classic
#include "AbsBeamline/VerticalFFAMagnet.h" // classic
#include "Elements/OpalVerticalFFAMagnet.h"
std::string OpalVerticalFFAMagnet::docstring_m =
std::string("The \"VerticalFFAMagnet\" element defines a vertical FFA ")+
std::string("magnet, which has a field that increases in the vertical ")+
std::string("direction while maintaining similar orbits.");
OpalVerticalFFAMagnet::OpalVerticalFFAMagnet() :
OpalElement(SIZE, "VERTICALFFAMAGNET", docstring_m.c_str()) {
itsAttr[B0] = Attributes::makeReal
("B0", "The nominal dipole field of the magnet at zero height [T].");
itsAttr[FIELD_INDEX] = Attributes::makeReal("FIELD_INDEX",
"Exponent term in the field index [m^(-1)].");
itsAttr[WIDTH] = Attributes::makeReal("WIDTH",
"The full width of the magnet. Particles moving more than WIDTH/2 horizontally, in either direction, are out of the aperture.");
itsAttr[MAX_HORIZONTAL_POWER] = Attributes::makeReal("MAX_HORIZONTAL_POWER",
"The maximum power in horizontal coordinate that will be considered in the field expansion.");
itsAttr[END_LENGTH] = Attributes::makeReal("END_LENGTH",
"The end length of the FFA fringe field [m].");
itsAttr[CENTRE_LENGTH] = Attributes::makeReal("CENTRE_LENGTH",
"The centre length of the FFA (i.e. length of the flat top) [m].");
itsAttr[BB_LENGTH] = Attributes::makeReal("BB_LENGTH",
"Determines the length of the bounding box. Magnet is situated symmetrically in the bounding box. [m]");
itsAttr[HEIGHT_POS_EXTENT] = Attributes::makeReal("HEIGHT_POS_EXTENT",
"Height of the magnet above z=0. Particles moving upwards more than HEIGHT_POS_EXTENT are out of the aperture [m].");
itsAttr[HEIGHT_NEG_EXTENT] = Attributes::makeReal("HEIGHT_NEG_EXTENT",
"Height of the magnet below z=0. Particles moving downwards more than HEIGHT_NEG_EXTENT are out of the aperture [m].");
registerRealAttribute("B0");
registerRealAttribute("FIELD_INDEX");
registerRealAttribute("WIDTH");
registerRealAttribute("MAX_HORIZONTAL_POWER");
registerRealAttribute("END_LENGTH");
registerRealAttribute("CENTRE_LENGTH");
registerRealAttribute("BB_LENGTH");
registerRealAttribute("HEIGHT_NEG_EXTENT");
registerRealAttribute("HEIGHT_POS_EXTENT");
registerOwnership();
VerticalFFAMagnet* magnet = new VerticalFFAMagnet("VerticalFFAMagnet");
magnet->setEndField(new endfieldmodel::Tanh(1., 1., 1));