RBend.h 15.5 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
#ifndef CLASSIC_RBend_HH
#define CLASSIC_RBend_HH

// ------------------------------------------------------------------------
// $RCSfile: RBend.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1.2.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: RBend
//   Defines the abstract interface for a rectangular bend magnet.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 18:57:53 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/Component.h"
#include "BeamlineGeometry/RBendGeometry.h"
#include "Fields/BMultipoleField.h"
#include "Physics/Physics.h"
#include "Algorithms/PartPusher.h"
29
#include "Fields/Fieldmap.h"
gsell's avatar
gsell committed
30 31 32 33
#include <vector>

class Fieldmap;

34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
/*
 * 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.
 */
gsell's avatar
gsell committed
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92

class RBend: public Component {

public:

    /// Constructor with given name.
    explicit RBend(const std::string &name);

    RBend();
    RBend(const RBend &);
    virtual ~RBend();

    /// Apply visitor to RBend.
    virtual void accept(BeamlineVisitor &) const;

93 94 95 96 97
    /*
     * Methods for OPAL-MAP
     * ====================
     */

gsell's avatar
gsell committed
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
    /// Get dipole field of RBend.
    virtual double getB() const = 0;

    /// Get RBend geometry.
    //  Version for non-constant object.
    virtual RBendGeometry &getGeometry() = 0;

    /// Get RBend geometry
    //  Version for constant object.
    virtual const RBendGeometry &getGeometry() const = 0;

    /// Get multipole expansion of field.
    //  Version for non-constant object.
    virtual BMultipoleField &getField() = 0;

    /// Get multipole expansion of field.
    //  Version for constant object.
    virtual const BMultipoleField &getField() const = 0;

    /// Get normal component.
    //  Return the normal component of order [b]n[/b] in T/m**(n-1).
    //  If [b]n[/b] is larger than the maximum order, the return value is zero.
    double getNormalComponent(int) const;

    /// Get skew component.
    //  Return the skew component of order [b]n[/b] in T/m**(n-1).
    //  If [b]n[/b] is larger than the maximum order, the return value is zero.
    double getSkewComponent(int) const;

    /// Set normal component.
    //  Set the normal component of order [b]n[/b] in T/m**(n-1).
    //  If [b]n[/b] is larger than the maximum order, the component is created.
    void setNormalComponent(int, double);

    /// Set skew component.
    //  Set the skew component of order [b]n[/b] in T/m**(n-1).
    //  If [b]n[/b] is larger than the maximum order, the component is created.
    void setSkewComponent(int, double);

    /// Get pole entry face rotation.
    //  Return the rotation of the entry pole face with respect to the x-axis.
    //  A positive angle rotates the pole face normal away from the centre
    //  of the machine.
    virtual double getEntryFaceRotation() const = 0;

    /// Get exit pole face rotation.
    //  Return the rotation of the exit pole face with respect to the x-axis.
    //  A positive angle rotates the pole face normal away from the centre
    //  of the machine.
    virtual double getExitFaceRotation() const = 0;

    /// Get entry pole face curvature.
    //  Return the curvature of the entry pole face.
    //  A positive curvature creates a convex pole face.
    virtual double getEntryFaceCurvature() const = 0;

    /// Get exit pole face curvature.
    //  Return the curvature of the exit pole face.
    //  A positive curvature creates a convex pole face.
    virtual double getExitFaceCurvature() const = 0;

    /// Get number of slices.
    //  Slices and stepsize used to determine integration step.
    virtual double getSlices() const = 0;

    /// Get stepsize.
    //  Slices and stepsize used to determine integration step.
    virtual double getStepsize() const = 0;


168 169 170 171
    /*
     * Methods for OPAL-SLICE.
     */
    virtual void addKR(int i, double t, Vector_t &K);
gsell's avatar
gsell committed
172 173 174
    virtual void addKT(int i, double t, Vector_t &K);


175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
    /*
     * 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.
gsell's avatar
gsell committed
197 198
    virtual bool bends() const;

kraus's avatar
kraus committed
199 200
    virtual void goOnline(const double &kineticEnergy);

201 202
    virtual void finalise();
    virtual void getDimensions(double &sBegin, double &sEnd) const;
203
    virtual ElementBase::ElementType getType() const;
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219
    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);
220
    void SetFieldAmplitude(double k0, double k0s);
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236

    /*
     * 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);
gsell's avatar
gsell committed
237

238 239
    void resetReinitializeFlag();
    void resetRecalcRefTrajFlag();
240

gsell's avatar
gsell committed
241 242 243 244
private:

    // Not implemented.
    void operator=(const RBend &);
Steve Russell's avatar
Steve Russell committed
245

246 247 248 249 250 251 252 253 254 255 256 257 258 259
    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,
260
                          double &engeFuncSecDerivNorm);
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
    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).
310
    double aperture_m;          /// Aperture of magnet in non-bend plane.
311 312 313
    double designEnergy_m;      /// Bend design energy (eV).
    double designRadius_m;      /// Bend design radius (m).
    double fieldAmplitude_m;    /// Amplitude of magnet field (T).
314 315
    double bX_m;                /// Amplitude of Bx field (T).
    double bY_m;                /// Amplitude of By field (T).
316 317 318 319
    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).
320
    double fieldIndex_m;        /// Dipole field index.
321 322 323 324 325 326 327 328 329 330 331
    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;
Steve Russell's avatar
Steve Russell committed
332

333 334 335 336 337 338 339 340 341 342 343 344 345 346 347
    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;
Steve Russell's avatar
Steve Russell committed
348

349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
    /// 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;
gsell's avatar
gsell committed
412 413
};

414 415 416 417 418 419 420 421 422 423
inline
void RBend::resetReinitializeFlag() {
    reinitialize_m = true;
}

inline
void RBend::resetRecalcRefTrajFlag() {
    recalcRefTraj_m = true;
}

424
#endif // CLASSIC_RBend_HH