Bend.h 11.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
#ifndef CLASSIC_BEND_H
#define CLASSIC_BEND_H

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

24
#include "AbsBeamline/BendBase.h"
25
#include "Algorithms/PartPusher.h"
26 27 28 29
#include "Utilities/GeneralClassicException.h"

#include "gsl/gsl_spline.h"
#include "gsl/gsl_interp.h"
30

ext-rogers_c's avatar
ext-rogers_c committed
31 32
#include <gtest/gtest_prod.h>

33
#include <vector>
34 35

class Fieldmap;
36
class MeshData;
37 38 39 40 41 42 43 44 45 46

/*
 * Class Bend
 *
 * Interface for general bend magnet.
 *
 * ------------------------------------------------------------------------
 *
 */

47
class Bend: public BendBase {
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 78 79 80

public:

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

    Bend();
    Bend(const Bend &);
    virtual ~Bend();

    /// Apply visitor to Bend.
    virtual void accept(BeamlineVisitor &) const = 0;

    /*
     * Methods for OPAL-SLICE.
     */
    virtual void addKR(int i, double t, Vector_t &K) { };
    virtual void addKT(int i, double t, Vector_t &K) { };


    /*
     * Methods for OPAL-T.
     * ===================
     */

    /// Apply field to particles with coordinates in magnet frame.
    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,
81
                       const Vector_t &P,
82 83 84 85
                       const double &t,
                       Vector_t &E,
                       Vector_t &B);

86 87 88 89 90
    virtual bool applyToReferenceParticle(const Vector_t &R,
                                          const Vector_t &P,
                                          const double &t,
                                          Vector_t &E,
                                          Vector_t &B);
91 92 93 94 95 96 97 98

    virtual void goOnline(const double &kineticEnergy);

    virtual void finalise();
    virtual void getDimensions(double &sBegin, double &sEnd) const;
    virtual ElementBase::ElementType getType() const = 0;
    virtual void initialise(PartBunch *bunch,
                            double &startField,
99 100 101 102 103 104 105 106 107
                            double &endField);

    // double getBendAngle() const;
    double getBendRadius() const;
    double getEffectiveCenter() const;
    double getEffectiveLength() const;

    double getStartElement() const;

108 109 110 111 112 113 114 115 116
    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.
117
    void setK1(double k1);
118 119 120 121

    void resetReinitializeFlag();
    void resetRecalcRefTrajFlag();

122 123 124
    double getExitAngle() const;
    double getMapLength() const;

125
    std::pair<Vector_t, Vector_t> getDesignPathSecant(double startsAtDistFromEdge, double length) const;
126 127 128 129 130 131 132

    std::vector<Vector_t> getOutline() const;
    MeshData getSurfaceMesh() const;

    virtual CoordinateSystemTrafo getBeginToEnd() const;

    virtual bool isInside(const Vector_t &r) const;
133 134 135 136 137 138 139
protected:
    void setMessageHeader(const std::string & header);
    double getStartField() const;
    Fieldmap* getFieldmap();

private:

ext-rogers_c's avatar
ext-rogers_c committed
140 141
    FRIEND_TEST(Maxwell, Zeros);
    
142 143 144
    // Not implemented.
    void operator=(const Bend &);

145 146 147
    void adjustFringeFields(double ratio);
    double calculateBendAngle();
    void calcEngeFunction(double zNormalized,
148 149 150 151 152
                          const std::vector<double> &engeCoeff,
                          int polyOrder,
                          double &engeFunc,
                          double &engeFuncDeriv,
                          double &engeFuncSecDerivNorm);
153 154 155 156 157 158 159 160 161 162 163
    Vector_t calcCentralField(const Vector_t &R,
                              double deltaX);
    Vector_t calcEntranceFringeField(const Vector_t &R,
                                     double deltaX);
    Vector_t calcExitFringeField(const Vector_t &R,
                                 double deltaX);
    bool calculateMapField(const Vector_t &R,
                           Vector_t &B);
    void calculateRefTrajectory(double &angleX,
                                double &angleY);
    double estimateFieldAdjustmentStep(double actualBendAngle,
164 165
                                       double mass,
                                       double betaGamma);
166 167 168
    void findBendEffectiveLength(double startField,
                                 double endField);
    void findBendStrength(double mass,
169 170 171
                          double gamma,
                          double betaGamma,
                          double charge);
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
    virtual bool findChordLength(Inform &msg,
                                 double &chordLength) = 0;
    bool findIdealBendParameters(double chordLength);
    void findReferenceExitOrigin(double &x, double &z);
    bool initializeFieldMap(Inform &msg);
    bool inMagnetCentralRegion(const Vector_t &R) const;
    bool inMagnetEntranceRegion(const Vector_t &R) const;
    bool inMagnetExitRegion(const Vector_t &R) const;
    bool isPositionInEntranceField(const Vector_t &R) const;
    bool isPositionInExitField(const Vector_t &R) const;
    void print(Inform &msg, double bendAngleX, double bendAngle);
    void readFieldMap(Inform &msg);
    bool reinitialize();
    void setBendEffectiveLength(double startField, double endField);
    void setBendStrength();
    void setEngeOriginDelta(double delta);
    void setFieldCalcParam();
    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, double chordlength);
195 196
    void retrieveDesignEnergy(double startField);

197 198
    CoordinateSystemTrafo getBeginToEnd_local() const;

199 200 201 202
    std::string messageHeader_m;

    BorisPusher pusher_m;       /// Pusher used to integrate reference particle
    /// through the bend.
203

204 205
    Fieldmap *fieldmap_m;       /// Magnet field map.
    bool fast_m;                /// Flag to turn on fast field calculation.
206

207
    double designRadius_m;      /// Bend design radius (m).
208

209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
    /// 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.
    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;

    bool recalcRefTraj_m;       /// Re-calculate reference trajectory.

    /*
     * 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;

248 249 250 251 252
    gsl_spline** entryFieldValues_m;
    gsl_spline** exitFieldValues_m;
    gsl_interp_accel *entryFieldAccel_m;
    gsl_interp_accel *exitFieldAccel_m;

253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
    /*
     * All coordinates are with respect to (x, z) = (0, 0). It is
     * assumed the ideal reference trajectory passes through this point.
     */
    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 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;
277 278 279 280 281 282 283 284
    double tanEntranceAngle_m;
    double tanExitAngle_m;

    CoordinateSystemTrafo beginToEnd_m;
    CoordinateSystemTrafo beginToEnd_lcs_m; // local coordinate system

    CoordinateSystemTrafo computeAngleTrafo_m;
    double maxAngle_m;
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 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327

inline
void Bend::finalise() {
    online_m = false;
}

inline
void Bend::getDimensions(double &sBegin, double &sEnd) const {
    sBegin = startField_m;
    sEnd = endField_m;
}

inline
double Bend::getBendRadius() const {
    return designRadius_m;
}

inline
double Bend::getEffectiveCenter() const {
    return elementEdge_m + designRadius_m * angle_m / 2.0;
}

inline
double Bend::getEffectiveLength() const {
    return designRadius_m * angle_m;
}

inline
double Bend::getStartElement() const {
    return elementEdge_m;
}

inline
void Bend::setK1(double k1) {
    if (std::abs(k1) > 0.0) {
        throw GeneralClassicException("Bend::setK1",
                                      "Quadrupole field temporarily not supported");
    }
    fieldIndex_m = k1;
}

328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
inline
void Bend::resetReinitializeFlag() {
    reinitialize_m = true;
}

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

inline
void Bend::setMessageHeader(const std::string & header)
{
    messageHeader_m = header;
}

inline
double Bend::getStartField() const
{
    return startField_m;
}

inline
Fieldmap* Bend::getFieldmap()
{
    return fieldmap_m;
}

inline
357
double Bend::getExitAngle() const
358
{
359
    return exitAngle_m;
360 361 362
}

inline
363
void Bend::setExitAngle(double angle)
364
{
365
    exitAngle_m = angle;
366 367 368
}

inline
369
double Bend::getMapLength() const
370
{
371
    return exitParameter2_m - entranceParameter2_m;
372 373 374
}

inline
375
CoordinateSystemTrafo Bend::getBeginToEnd() const
376
{
377
    return beginToEnd_m;
378 379 380
}

inline
381
CoordinateSystemTrafo Bend::getBeginToEnd_local() const
382
{
383
    return beginToEnd_lcs_m;
384 385
}

ext-rogers_c's avatar
ext-rogers_c committed
386
#endif // CLASSIC_BEND_H