Commit 08556dc2 authored by Steve Russell's avatar Steve Russell
Browse files

Update of the RBend and SBend elements. Routines are now more robust and accurate. Input angles are

now in radians instead of degrees.
parent 9815e7b0
This diff is collapsed.
......@@ -31,12 +31,50 @@
class Fieldmap;
// Class RBend
// ------------------------------------------------------------------------
/// Interface for rectangular bend.
// Class RBend defines the abstract interface for rectangular bend magnets.
// A rectangular bend magnet has a rectilinear geometry about which its
// multipole components are specified.
/*
* Class RBend
*
* Interface for a rectangular bend magnet.
*
* A rectangular bend magnet physically has a rectangular shape.
*
* The standard rectangular magnet, for purposes of definitions, has a field in
* the y direction. This produces a bend in the horizontal (x) plane. Bends in
* other planes can be accomplished by rotating the magnet about the axes.
*
* A positive bend angle is defined as one that bends a beam to the right when
* looking down (in the negative y direction) so that the beam is bent in the
* negative x direction. (This definition of a positive bend is the same whether
* the charge is positive or negative.)
*
* A zero degree entrance edge angle is parallel to the x direction in an x/y/s
* coordinate system. A positive entrance edge angle is defined as one that
* rotates the positive edge (in x) of the angle toward the positive s axis.
*
* Since the magnet geometry is a fixed rectangle, the exit edge angle is
* defined by the bend angle of the magnet and the entrance edge angle. In
* general, the exit edge angle is equal to the bend angle minus the entrance
* edge angle.
*
* ------------------------------------------------------------------------
*
* This class defines three interfaces:
*
* 1) Interface for rectangular magnets for OPAL-MAP.
*
* Here we specify multipole components about the curved magnet trajectory.
*
*
* 2) Interface for rectangular magnets in OPAL-SLICE.
*
* Here we define a rectangular magnet for use in the OPAL
* slice model.
*
*
* 3) Interface for rectangular magnets for OPAL-T.
*
* Here we defined the magnet as a field map.
*/
class RBend: public Component {
......@@ -52,6 +90,11 @@ public:
/// Apply visitor to RBend.
virtual void accept(BeamlineVisitor &) const;
/*
* Methods for OPAL-MAP
* ====================
*/
/// Get dipole field of RBend.
virtual double getB() const = 0;
......@@ -121,159 +164,246 @@ public:
// Slices and stepsize used to determine integration step.
virtual double getStepsize() const = 0;
virtual void addKR(int i, double t, Vector_t &K);
/*
* Methods for OPAL-SLICE.
*/
virtual void addKR(int i, double t, Vector_t &K);
virtual void addKT(int i, double t, Vector_t &K);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
virtual void initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor);
virtual void finalise();
/*
* Methods for OPAL-T.
*/
/// Apply field to particles with coordinates in magnet frame.
virtual bool apply(const size_t &i,
const double &t,
double E[],
double B[]);
virtual bool apply(const size_t &i,
const double &t,
Vector_t &E,
Vector_t &B);
/// Apply field to particles in beam frame.
virtual bool apply(const Vector_t &R,
const Vector_t &centroid,
const double &t,
Vector_t &E,
Vector_t &B);
/// Indicates that element bends the beam.
virtual bool bends() const;
void setBendAngle(const double &angle);
void setAmplitudem(double vPeak);
void setFullGap(const double &gap);
void setLength(const double &length);
void setLongitudinalRotation(const double &rotation);
void setLongitudinalRotation(const double &k0, const double &k0s);
/// Set the name of the "field map";
/// For a bend this file contains the coefficients for the Enge function
void setFieldMapFN(std::string fmapfn);
std::string getFieldMapFN() const;
void setEngeCoefs(const std::vector<double> EngeCoefs);
void setAlpha(const double &alpha);
void setBeta(const double &beta);
void setDesignEnergy(const double &energy);
virtual void finalise();
virtual void getDimensions(double &sBegin, double &sEnd) const;
virtual const std::string &getType() const;
virtual void getDimensions(double &zBegin, double &zEnd) const;
double getEffectiveLength() const;
double getEffectiveCenter() const;
double getBendAngle() const;
double getStartElement() const;
double getR() const;
virtual void initialise(PartBunch *bunch,
double &startField,
double &endField,
const double &scaleFactor);
double GetBendAngle() const;
double GetBendRadius() const;
double GetEffectiveCenter() const;
double GetEffectiveLength() const;
std::string GetFieldMapFN() const;
double GetStartElement() const;
void SetAperture(double aperture);
void SetBendAngle(double angle);
void SetBeta(double beta);
void SetDesignEnergy(double energy);
void SetEntranceAngle(double entranceAngle);
void SetFieldAmplitude(double fieldAmplitude);
/*
* 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.
*/
void SetFieldMapFN(std::string fileName);
void SetFullGap(double gap);
/// Set quadrupole field component.
void SetK1(double k1);
void SetLength(double length);
/// Set rotation about z axis in bend frame.
void SetRotationAboutZ(double rotation);
void SetRotationAboutZ(double k0, double k0s);
private:
// Magnet length.
double length_m;
// Magnet gap.
double gap_m;
// Flag to reinitialize the field the first time the
// magnet is applied.
bool reinitialize_m;
// Not implemented.
void operator=(const RBend &);
std::string filename_m; /**< The name of the inputfile*/
Fieldmap *myFieldmap;
double scale_m; /**< scale multiplier*/
Fieldmap *fieldmap_m;
double amplitude_m;
Vektor<double, 2> field_orientation_m; // (cos(alpha), sin(alpha))
double startField_m;
double endField_m;
bool fast_m;
double sin_face_alpha_m; // alpha is the angle between the projection of the normal of the face onto the
double cos_face_alpha_m; // s-u plane
double tan_face_alpha_m;
double sin_face_beta_m; // beta is defined as arctan(|n_parallel|/|n_perpendicular|) where n_parallel is the component
double cos_face_beta_m; // of the normal of the face which is parallel to the s-u plane and n_perpendicular is
double tan_face_beta_m; // perpendicular to it
double design_energy_m;
double angle_m; // Bend angle of central trajectory particle (radians).
double *map_m;
int map_size_m;
double map_step_size_m;
BorisPusher pusher_m;
double startElement_m;
double R_m;
// Effective length of hard edge dipole approximation.
double effectiveLength_m;
// Effective center of hard edge dipole approximation.
//
// This is defined as the z point at which the field integral
// is half of the full field integral starting at the upstream
// end of the dipole.
double effectiveCenter_m;
static int RBend_counter_m;
void calculateEffectiveLength();
void calculateEffectiveCenter();
// Set the bend strength based on the desired bend angle
// and the reference energy.
void setBendStrength();
void AdjustFringeFields(double ratio);
double CalculateBendAngle();
void CalcCentralField(Vector_t R,
double deltaX,
double angle,
Vector_t &B);
void CalcDistFromRefOrbitCentralField(Vector_t R,
double &deltaX,
double &angle);
void CalcEngeFunction(double zNormalized,
std::vector<double> engeCoeff,
int polyOrder,
double &engeFunc,
double &engeFuncDeriv,
double &engeFuncSecDeriv);
void CalcEntranceFringeField(Vector_t REntrance,
double deltaX,
Vector_t &B);
void CalcExitFringeField(Vector_t RExit, double deltaX, Vector_t &B);
void CalculateMapField(Vector_t R, Vector_t &E, Vector_t &B);
void CalculateRefTrajectory(double &angleX, double &angleY);
double EstimateFieldAdjustmentStep(double actualBendAngle,
double mass,
double betaGamma);
void FindBendEffectiveLength(double startField, double endField);
bool FindBendLength(Inform &msg,
double &bendLength,
bool &bendLengthFromMap);
void FindBendStrength(double mass,
double gamma,
double betaGamma,
double charge);
bool FindIdealBendParameters(double bendLength);
void FindReferenceExitOrigin(double &x, double &z);
bool InitializeFieldMap(Inform &msg);
bool InMagnetCentralRegion(Vector_t R, double &deltaX, double &angle);
bool InMagnetEntranceRegion(Vector_t R, double &deltaX);
bool InMagnetExitRegion(Vector_t R, double &deltaX);
bool IsPositionInEntranceField(Vector_t R, Vector_t &REntrance);
bool IsPositionInExitField(Vector_t R, Vector_t &RExit);
void Print(Inform &msg, double bendAngleX, double bendAngle);
void ReadFieldMap(Inform &msg);
bool Reinitialize();
Vector_t RotateOutOfBendFrame(Vector_t X);
Vector_t RotateToBendFrame(Vector_t X);
void SetBendEffectiveLength(double startField, double endField);
void SetBendStrength();
void SetEngeOriginDelta(double delta);
void SetFieldCalcParam(bool lengthFromMap);
void SetGapFromFieldMap();
bool SetupBendGeometry(Inform &msg, double &startField, double &endField);
bool SetupDefaultFieldMap(Inform &msg);
void SetFieldBoundaries(double startField, double endField);
void SetupPusher(PartBunch *bunch);
bool TreatAsDrift(Inform &msg);
BorisPusher pusher_m; /// Pusher used to integrate reference particle
/// through the bend.
std::string fileName_m; /// Name of field map that defines magnet.
Fieldmap *fieldmap_m; /// Magnet field map.
bool fast_m; /// Flag to turn on fast field calculation.
/// (Not currently used.)
double angle_m; /// Bend angle for reference particle with bend
/// design energy (radians).
double aperture_m; /// Aperture of magnet in non-bend (horizontal)
/// plane.
double designEnergy_m; /// Bend design energy (eV).
double designRadius_m; /// Bend design radius (m).
double fieldAmplitude_m; /// Amplitude of magnet field (T).
double entranceAngle_m; /// Angle between incoming reference trajectory
/// 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 gradient_m; /// Quadrupole component of field.
double elementEdge_m; /// Physical start of magnet in s coordinates (m).
double startField_m; /// Start of magnet field map in s coordinates (m).
double endField_m; /// End of magnet field map in s coordinates (m).
/*
* Flag to reinitialize the bend the first time the magnet
* is called. This redefines the design energy of the bend
* to the current average beam energy, keeping the bend angle
* constant.
*/
bool reinitialize_m;
// Calculate bend angle given reference energy and field strength.
double calculateBendAngle(double bendLength);
bool recalcRefTraj_m; /// Re-calculate reference trajectory.
/*
* These two parameters are used to set up the bend geometry.
*
* When using the default, internal field map, the gap of the
* magnet must be set in the input file. Otherwise, this parameter
* is defined by the magnet field map.
*
* The magnet length is used to define the length of the magnet.
* This parameter must be set in the input file when using the default,
* internal field map. Otherwise defining it is optional.
*/
double length_m;
double gap_m;
// Calculate the reference particle trajectory map.
double calculateRefTrajectory(const double zBegin);
/// Map of reference particle trajectory.
std::vector<double> refTrajMapX_m;
std::vector<double> refTrajMapY_m;
std::vector<double> refTrajMapZ_m;
int refTrajMapSize_m;
double refTrajMapStepSize_m;
/*
* Enge function field map members.
*/
/*
* Entrance and exit position parameters. Ultimately they are used to
* determine the origins of the entrance and exit edge Enge functions and
* the extent of the field map. However, how they are used to do this
* depends on how the bend using the map is setup in the OPAL input file.
* So, we use generic terms to start.
*/
double entranceParameter1_m;
double entranceParameter2_m;
double entranceParameter3_m;
double exitParameter1_m;
double exitParameter2_m;
double exitParameter3_m;
/// Enge coefficients for map entry and exit regions.
std::vector<double> engeCoeffsEntry_m;
std::vector<double> engeCoeffsExit_m;
/*
* All coordinates are with respect to (x, z) = (0, 0). It is
* assumed the ideal reference trajectory passes through this point.
*/
double xOriginEngeEntry_m; /// x coordinate of entry Enge function origin.
double zOriginEngeEntry_m; /// z coordinate of entry Enge function origin.
double deltaBeginEntry_m; /// Perpendicular distance from entrance Enge
/// function origin where Enge function starts.
double deltaEndEntry_m; /// Perpendicular distance from entrance Enge
/// function origin that Enge function ends.
int polyOrderEntry_m; /// Enge function order for entry region.
/*
* The ideal reference trajectory passes through (xExit_m, zExit_m).
*/
double xExit_m;
double zExit_m;
double xOriginEngeExit_m; /// x coordinate of exit Enge function origin.
double zOriginEngeExit_m; /// z coordinate of exit Enge function origin.
double deltaBeginExit_m; /// Perpendicular distance from exit Enge
/// function origin that Enge function starts.
double deltaEndExit_m; /// Perpendicular distance from exit Enge
/// function origin that Enge function ends.
int polyOrderExit_m; /// Enge function order for entry region.
double cosEntranceAngle_m;
double sinEntranceAngle_m;
double exitEdgeAngle_m; /// Exit edge angle (radians.
double cosExitAngle_m;
double sinExitAngle_m;
// For OPAL-SLICE.
static int RBend_counter_m;
};
inline void RBend::setAlpha(const double &alpha) {
Orientation_m(0) = -alpha * Physics::pi / 180.0;
sin_face_alpha_m = sin(Orientation_m(0));
cos_face_alpha_m = cos(Orientation_m(0));
tan_face_alpha_m = tan(Orientation_m(0));
}
inline void RBend::setBeta(const double &beta) {
Orientation_m(1) = beta * Physics::pi / 180.0;
sin_face_beta_m = sin(Orientation_m(1));
cos_face_beta_m = cos(Orientation_m(1));
tan_face_beta_m = tan(Orientation_m(1));
}
inline void RBend::setDesignEnergy(const double &energy)
{ design_energy_m = energy; }
inline void RBend::setLongitudinalRotation(const double &rotation) {
Orientation_m(2) = rotation * Physics::pi / 180.0;
}
inline void RBend::setLongitudinalRotation(const double &k0, const double &k0s) {
Orientation_m(2) = atan2(k0s, k0);
amplitude_m = sqrt(pow(k0, 2.0) + pow(k0s, 2.0));
if(fabs(k0s) > 1.e-8) {
if(amplitude_m > 1.e-8) {
field_orientation_m(0) = k0 / amplitude_m;
field_orientation_m(1) = k0s / amplitude_m;
}
} else {
field_orientation_m(0) = 1.0;
field_orientation_m(1) = 0.0;
}
}
#endif // CLASSIC_RBend_HH
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -3,110 +3,85 @@
#include "Fields/Fieldmap.hh"
// Class FM1DProfile1
//---------------------------------------------------------------------------
/// Field definition for 1D representation of bending magnet.
//
// Class FM1DProfile1 defines a 1D field map for use in bending magnets.
/*
* Class FM1DProfile.
* --------------------------------------------------------------------------
* Field definition for 1D representation of bending magnet.
*
* Class FM1DProfile1 defines a 1D field map for us in bending magnets.
*/
class FM1DProfile1: public Fieldmap {
public:
/// Calculates the normalized (to 1) field strength on axis and the first and
/// second derivatives (strength[0], stength[1], strength[2]) of the profile
/// given the position, X.
virtual bool getFieldstrength(const Vector_t &X, Vector_t &strength, Vector_t &info) const;
virtual bool getFieldDerivative(const Vector_t &X, Vector_t &E, Vector_t &B, const DiffDirection &dir) const;
/// Get field dimensions.
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const;
virtual void getFieldDimensions(double &xIni, double &xFinal, double &yIni, double &yFinal, double &zIni, double &zFinal) const;
virtual void swap();
/// Get field information.
virtual void getInfo(Inform *);
virtual bool getFieldDerivative(const Vector_t &X,
Vector_t &E,
Vector_t &B,
const DiffDirection &dir) const;
virtual void Get1DProfile1EntranceParam(double &entranceParameter1,
double &entranceParameter2,
double &entranceParameter3);
virtual void Get1DProfile1ExitParam(double &exitParameter1,
double &exitParameter2,
double &exitParameter3);
virtual double GetFieldGap();
virtual void getFieldDimensions(double &zBegin,
double &zEnd,
double &rBegin,
double &rEnd) const;
virtual void getFieldDimensions(double &xIni,
double &xFinal,
double &yIni,
double &yFinal,
double &zIni,
double &zFinal) const;
virtual bool getFieldstrength(const Vector_t &X,
Vector_t &strength,
Vector_t &info) const;
virtual double getFrequency() const;
virtual void getInfo(Inform *);
virtual void setFrequency(double freq);
virtual void setExitFaceSlope(const double &);
/// Set exit edge constants.
virtual void setEdgeConstants(const double &bendAngle, const double &entranceAngle, const double &exitAngle);
/// Set magnet gap.
virtual void setFieldGap(const double &gap);
/// Set distance between Enge function zeros.
virtual void setFieldLength(const double &length);
virtual void Get1DProfile1EngeCoeffs(std::vector<double> &engeCoeffsEntry,
std::vector<double> &engeCoeffsExit);
virtual void swap();
/// Adjust the extend of the fringe field regions. This is used
/// to make sure we don't see funny jumps in the field as we
/// move into and out of the magnet.
virtual bool adjustFringeFields();
virtual void SetFieldGap(double gap);
private:
/// Constructor with field map file name.
FM1DProfile1(std::string aFilename);
FM1DProfile1(std::string Filename);
virtual ~FM1DProfile1();
/// Read field map.
virtual void freeMap();
virtual void readMap();
/// Free field map.
virtual void freeMap();
/*
* Entrance and exit position parameters. These are read in from the input
* input file. Ultimately they are used to determine the origin of the
* field Enge function and the extent of the field map. However, how they
* are used to do this depends on how the bend using the map is setup in the
* OPAL input file. So, we use generic terms to start.
*/
double entranceParameter1_m;
double entranceParameter2_m;
double entranceParameter3_m;
double exitParameter1_m;
double exitParameter2_m;
double exitParameter3_m;
/// Enge coefficients for map entry and exit regions.
double *EngeCoefs_entry_m;
double *EngeCoefs_exit_m;
/// Entry position of field map entry region.
double zbegin_entry_m;
/// End position of field map entry region.
double zend_entry_m;
/// Origin position for Enge function of field map entry region.
double polynomialOrigin_entry_m;
/// Enge function order for entry region.
int polynomialOrder_entry_m;
/// Entry position of field map exit region.
double zbegin_exit_m;
/// End position of field map exit region.
double zend_exit_m;
/// Origin position for Enge function of field map entry region.
double polynomialOrigin_exit_m;
/// Enge functioin order for entry region.
int polynomialOrder_exit_m;
/// Length of map (m).
double length_m;
/// Bend full gap height (m).
double gapHeight_m;
/// x position in local coordinate system where central trajectory intersects