MultipoleT.h 16.3 KB
Newer Older
ext-rogers_c's avatar
ext-rogers_c committed
1 2
/*
 *  Copyright (c) 2017, Titus Dascalu
3
 *  Copyright (c) 2018, Martin Duy Tat
ext-rogers_c's avatar
ext-rogers_c committed
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 29 30 31 32 33 34 35 36 37 38 39 40
 *  All rights reserved.
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *  1. Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *  3. Neither the name of STFC nor the names of its contributors may be used to
 *     endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 */


#ifndef CLASSIC_MULTIPOLET_H
#define CLASSIC_MULTIPOLET_H

/** ---------------------------------------------------------------------
  *
  * MultipoleT defines a straight or curved combined function magnet (up 
  * to arbitrary multipole component) with Fringe Fields
  *
  * ---------------------------------------------------------------------
  *
  * Class category: AbsBeamline \n
41
  * $Author: Titus Dascalu, Martin Duy Tat, Chris Rogers
ext-rogers_c's avatar
ext-rogers_c committed
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 78 79 80 81 82 83
  *
  * ---------------------------------------------------------------------
  *
  * The field is obtained from the scalar potential \n
  *     @f[ V = f_0(x,s) z + f_1 (x,s) \frac{z^3}{3!} + f_2 (x,s) \frac{z^5}{5!} +
  *     ...  @f] \n
  *     (x,z,s) -> Frenet-Serret local coordinates along the magnet \n
  *     z -> vertical component \n
  *     assume mid-plane symmetry \n 
  *     set field on mid-plane -> @f$ B_z = f_0(x,s) = T(x) \cdot S(s) @f$ \n
  *     T(x) -> transverse profile; this is a polynomial describing
  *             the field expansion on the mid-plane inside the magnet
  *             (not in the fringe field);
  *             1st term is the dipole strength, 2nd term is the 
  *             quadrupole gradient * x, etc. \n
  *          -> when setting the magnet, one gives the multipole
  *             coefficients of this polynomial (i.e. dipole strength,  
  *             quadrupole gradient, etc.) \n
  * \n
  * ------------- example ----------------------------------------------- \n
  *     Setting a combined function magnet with dipole, quadrupole and 
  *     sextupole components: \n
  *     @f$ T(x) = B_0 + B_1 \cdot x + B_2 \cdot x^2 @f$\n
  *     user gives @f$ B_0, B_1, B_2 @f$ \n
  * ------------- example end ------------------------------------------- \n
  * \n
  *     S(s) -> fringe field \n
  *     recursion -> @f$ f_n (x,s) = (-1)^n \cdot \sum_{i=0}^{n} C_n^i 
  *     \cdot T^{(2i)} \cdot S^{(2n-2i)} @f$ \n
  *     for curved magnets the above recursion is more complicated \n
  *     @f$ C_n^i @f$ -> binomial coeff; 
  *     @f$ T^{(n)} @f$ -> n-th derivative
  *
  * ---------------------------------------------------------------------
  */

#include "AbsBeamline/EndFieldModel/Tanh.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
#include "Fields/BMultipoleField.h"
#include "Algorithms/Vektor.h"
#include "AbsBeamline/Component.h"
#include "AbsBeamline/BeamlineVisitor.h"
84 85
#include "AbsBeamline/MultipoleTFunctions/RecursionRelation.h"
#include "AbsBeamline/MultipoleTFunctions/RecursionRelationTwo.h"
ext-rogers_c's avatar
ext-rogers_c committed
86 87 88 89
#include "gsl/gsl_sf.h"
#include <vector>

class MultipoleT: public Component {
90 91 92 93
public: 
    /** Constructor
     *  \param name -> User-defined name
     */
ext-rogers_c's avatar
ext-rogers_c committed
94 95 96 97 98 99
    explicit MultipoleT(const std::string &name);
    /** Copy constructor */
    MultipoleT(const MultipoleT &right);
    /** Destructor */ 
    ~MultipoleT();
    /** Inheritable copy constructor */
100
    ElementBase* clone() const override;
101
    /** Return a dummy field value */
102
    EMField &getField() override;
103
    /** Return a dummy field value */
104
    const EMField &getField() const override;
ext-rogers_c's avatar
ext-rogers_c committed
105
    /** Not implemented */
106
    void getDimensions(double &zBegin, double &zEnd) const override;
ext-rogers_c's avatar
ext-rogers_c committed
107
    /** Calculate the field at some arbitrary position
108 109 110 111 112 113 114
     *  If particle is outside field map true is returned,
     *  otherwise false is returned
     *  \param R -> Position in the lab coordinate system of the multipole
     *  \param P -> Not used
     *  \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
ext-rogers_c's avatar
ext-rogers_c committed
115 116
     */
    bool apply(const Vector_t &R, const Vector_t &P, const double &t,
117
               Vector_t &E, Vector_t &B) override;
ext-rogers_c's avatar
ext-rogers_c committed
118
    /** Calculate the field at the position of the ith particle
119
     *  \param i -> Index of the particle event override; field is calculated at this
120 121 122 123 124 125
     *  position
     *  If particle is outside field map true is returned,
     *  otherwise false is returned
     *  \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
ext-rogers_c's avatar
ext-rogers_c committed
126
     */
127
    bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) override;
128 129 130 131 132 133 134
    /** Initialise the MultipoleT
     *  \param bunch -> Bunch the global bunch object
     *  \param startField -> Not used
     *  \param endField -> Not used
     */
    void initialise(PartBunchBase<double, 3>*,
                    double &startField,
135
                    double &endField) override;
136 137 138
    /** Initialises the geometry */
    void initialise();

139
    /** Finalise the MultipoleT - sets bunch to NULL */
140
    void finalise() override;
141
    /** Return true if dipole component not zero */
142
    bool bends() const override;
ext-rogers_c's avatar
ext-rogers_c committed
143
    /** Return the cell geometry */
144
    PlanarArcGeometry& getGeometry() override;
ext-rogers_c's avatar
ext-rogers_c committed
145
    /** Return the cell geometry */
146
    const PlanarArcGeometry& getGeometry() const override;
ext-rogers_c's avatar
ext-rogers_c committed
147
    /** Accept a beamline visitor */
148
    void accept(BeamlineVisitor& visitor) const override;
ext-rogers_c's avatar
ext-rogers_c committed
149 150 151 152
    /** Get the dipole constant B_0 */
    double getDipoleConstant() const;
    /** Set the dipole constant B_0 */
    void setDipoleConstant(double B0);
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
    /** Get the number of terms used in calculation of field components */
    std::size_t getMaxOrder() const;
    /** Set the number of terms used in calculation of field components \n
     *  Maximum power of z in Bz is 2 * maxOrder_m \n
     *  \param maxOrder -> Number of terms in expansion in z
     */
    void setMaxOrder(std::size_t maxOrder);
    /** Get highest power of x in polynomial expansions */
    std::size_t getMaxXOrder() const;
    /** Set the number of terms used in polynomial expansions
     *  \param maxXOrder -> Number of terms in expansion in z
     */
    void setMaxXOrder(std::size_t maxXOrder);
    /** Get the maximum order in the given transverse profile */
    std::size_t getTransMaxOrder() const;
ext-rogers_c's avatar
ext-rogers_c committed
168
    /** Set the maximum order in the given transverse profile
169 170 171
     *  \param transMaxOrder -> Highest power of x in field expansion 
     */
    void setTransMaxOrder(std::size_t transMaxOrder);
ext-rogers_c's avatar
ext-rogers_c committed
172 173
    /** Set transverse profile T(x)
      * T(x) = B_0 + B1 x + B2 x^2 + B3 x^3 + ...
174 175
      * \param n -> Order of the term (d^n/dx^n) to be set
      * \param Bn -> Value of transverse profile coefficient
ext-rogers_c's avatar
ext-rogers_c committed
176
      */
177 178 179 180
    void setTransProfile(std::size_t n, double Bn);
    /** Get transverse profile
     *  \param n -> Power of x
     */
ext-rogers_c's avatar
ext-rogers_c committed
181
    double getTransProfile(int n) const;
182
    /** Get all terms of transverse profile */
ext-rogers_c's avatar
ext-rogers_c committed
183
    std::vector<double> getTransProfile() const;
184 185 186 187 188 189 190 191
    /** Set fringe field model \n
     *  Tanh model used here \n
     *  @f[ 1/2 * \left [tanh \left( \frac{s + s_0}{\lambda_{left}} \right) 
     *  - tanh \left( \frac{s - s_0}{\lambda_{right}} \right) \right] @f] 
     *  \param s0 -> Centre field length
     *  \param \lambda_{left} -> Left end field length
     *  \param \lambda_{right} -> Right end field length
     */
ext-rogers_c's avatar
ext-rogers_c committed
192 193 194 195 196 197 198 199 200
    bool setFringeField(double s0, double lambda_left, double lambda_right);
    /** Return vector of 2 doubles
      * [left fringe length, right fringelength]  
      */
    std::vector<double> getFringeLength() const;
    /** Set the bending angle of the magnet */
    void setBendAngle(double angle);
    /** Get the bending angle of the magnet */
    double getBendAngle() const;
201 202 203
    /** Set the entrance angle
     *  \param entranceAngle -> Entrance angle
     */
ext-rogers_c's avatar
ext-rogers_c committed
204 205 206
    void setEntranceAngle(double entranceAngle);
    /** Get the entrance angle */
    double getEntranceAngle() const;
207 208
    /** Get the bending radius \n 
      * Not used, when needed radius is found from length_m / angle_m
ext-rogers_c's avatar
ext-rogers_c committed
209 210
      */
    double getBendRadius() const;
211 212 213
    /** Set the length of the magnet \n
      * If straight-> Actual length \n
      * If curved -> Arc length \n
ext-rogers_c's avatar
ext-rogers_c committed
214 215 216 217
      */
    void setLength(double length);
    /** Get the length of the magnet */
    double getLength() const;
218
    /** Not used */
ext-rogers_c's avatar
ext-rogers_c committed
219 220
    double getChordLength() const;
    /** Set the aperture dimensions
221 222 223
      * This element only supports a rectangular aperture
      * \param vertAp -> Vertical aperture length
      * \param horizAp -> Horisontal aperture length
ext-rogers_c's avatar
ext-rogers_c committed
224 225 226
      */
    void setAperture(double vertAp, double horizAp);
    /** Get the aperture dimensions
227
      * Returns a vector of 2 doubles
ext-rogers_c's avatar
ext-rogers_c committed
228 229
      */
    std::vector<double> getAperture() const;
230 231 232 233
    /** Set the angle of rotation of the magnet around its axis \n
     *  To make skew components
     *  \param rot -> Angle of rotation
     */
ext-rogers_c's avatar
ext-rogers_c committed
234 235 236 237 238 239 240
    void setRotation(double rot);
    /** Get the angle of rotation of the magnet around its axis */
    double getRotation() const;
    /** Set variable radius flag to true */
    void setVarRadius();
    /** Get the value of variableRadius_m */
    bool getVarRadius() const;
241 242 243 244 245 246
    /** Get distance between centre of magnet and entrance */
    double getBoundingBoxLength() const;
    /** Set distance between centre of magnet and enctrance
     *  \param boundingBoxLength -> Distance between centre and entrance
     */
    void setBoundingBoxLength(const double &boundingBoxLength);
ext-rogers_c's avatar
ext-rogers_c committed
247

248
private:
ext-rogers_c's avatar
ext-rogers_c committed
249 250
    MultipoleT operator=(const MultipoleT& rhs);
    // End fields
251 252 253 254 255 256 257 258 259 260 261 262 263 264
    endfieldmodel::Tanh fringeField_l; // Left
    endfieldmodel::Tanh fringeField_r; // Right
    /** Field expansion parameters \n
     *  Number of terms in z expansion used in calculating field components
     */
    std::size_t maxOrder_m = 0;
    /** Highest order of polynomial expansions in x */
    std::size_t maxOrderX_m = 0;
    /** Objects for storing differential operator acting on Fn */
    std::vector<polynomial::RecursionRelationTwo> recursion_VarRadius_m;
    std::vector<polynomial::RecursionRelation> recursion_ConstRadius_m;
    /** Highest power in given mid-plane field */
    std::size_t transMaxOrder_m = 0; 
    /** List of transverse profile coefficients */
ext-rogers_c's avatar
ext-rogers_c committed
265
    std::vector<double> transProfile_m;
266
    /** Geometry */
ext-rogers_c's avatar
ext-rogers_c committed
267 268
    PlanarArcGeometry planarArcGeometry_m;
    /** Rotate frame for skew elements
269 270 271 272 273
     *  Consecutive rotations:
     *  1st -> about central axis
     *  2nd -> azimuthal rotation
     *  \param R -> Vector to be rotated
     */
ext-rogers_c's avatar
ext-rogers_c committed
274
    Vector_t rotateFrame(const Vector_t &R);
275 276 277
    /** Inverse of the 1st rotation in rotateFrame() method \n
     *  Used to rotate B field back to global coordinate system
     */
ext-rogers_c's avatar
ext-rogers_c committed
278
    Vector_t rotateFrameInverse(Vector_t &B);
279
    /** Transform to Frenet-Serret coordinates for sector magnets */
ext-rogers_c's avatar
ext-rogers_c committed
280
    Vector_t transformCoords(const Vector_t &R);
281
    /** Magnet parameters */
ext-rogers_c's avatar
ext-rogers_c committed
282 283 284 285
    double length_m;
    double angle_m;
    double entranceAngle_m;
    double rotation_m;
286
    /** Variable radius flag */
ext-rogers_c's avatar
ext-rogers_c committed
287
    bool variableRadius_m;
288 289 290 291
    /** Distance between centre of magnet and entrance */
    double boundingBoxLength_m;
    /** Get field component methods
     */
ext-rogers_c's avatar
ext-rogers_c committed
292 293 294
    double getBx (const Vector_t &R);
    double getBz (const Vector_t &R);
    double getBs (const Vector_t &R);
295
    /** Assume rectangular aperture with these dimensions */
ext-rogers_c's avatar
ext-rogers_c committed
296 297
    double verticalApert_m;
    double horizApert_m;
298
    /** Not implemented */
ext-rogers_c's avatar
ext-rogers_c committed
299
    BMultipoleField dummy;
300 301 302 303
    /** Returns the value of the fringe field n-th derivative at s
     *  \param n -> nth derivative
     *  \param s -> Coordinate s
     */
ext-rogers_c's avatar
ext-rogers_c committed
304
    double getFringeDeriv(int n, double s);
305 306 307 308 309 310 311 312
    /** Returns the value of the transverse field n-th derivative at x
     *  \param n -> nth derivative
     *  \param x -> Coordinate x
     */
    double getTransDeriv(std::size_t n, double x);
    /** Tests if inside the magnet
     *  \param R -> Coordinate vector
     */
ext-rogers_c's avatar
ext-rogers_c committed
313
    bool insideAperture(const Vector_t &R);
314 315 316 317 318 319 320 321 322
    /** Radius of curvature \n
     *  If radius of curvature is infinite, -1 is returned \n
     *  If radius is constant, then \n
     *  @f$ \rho(s) = length_m / angle_m @f$ \n
     *  If radius is variable, then
     *  @f$ \rho(s) = rho(0) * S(0) / S(s) @f$
     *  where S(s) is the fringe field
     *  \param s -> Coordinate s
     */
ext-rogers_c's avatar
ext-rogers_c committed
323
    double getRadius(double s);
324 325 326 327
    /** Returns the scale factor @f$ h_s = 1 + x / \rho(s) @f$
     *  \param x -> Coordinate x
     *  \param s -> Coordinate s
     */
ext-rogers_c's avatar
ext-rogers_c committed
328
    double getScaleFactor(double x, double s);
329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351
    /** Calculate partial derivative of fn wrt x using a 5-point
     *  finite difference formula \n
     *  Error of order stepSize^4
     *  \param n -> nth derivative
     *  \param x -> Coordinate x
     *  \param s -> Coordinate s
     */
    double getFnDerivX(std::size_t n, double x, double s);
    /** Calculate partial derivative of fn wrt s using a 5-point
     *  finite difference formuln
     *  Error of order stepSize^4
     *  \param n -> nth derivative
     *  \param x -> Coordinate x
     *  \param s -> Coordinate s
     */
    double getFnDerivS(std::size_t n, double x, double s);
    /** Calculate fn(x, s) by expanding the differential operator
     *  (from Laplacian and scalar potential) in terms of polynomials
     *  \param n -> nth derivative
     *  \param x -> Coordinate x
     *  \param s -> Coordinate s
     */
    double getFn(std::size_t n, double x, double s);
ext-rogers_c's avatar
ext-rogers_c committed
352
};
353

ext-rogers_c's avatar
ext-rogers_c committed
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
inline
    void MultipoleT::setVarRadius() {
        variableRadius_m = true;
}
inline
    bool MultipoleT::getVarRadius() const {
        return variableRadius_m;
}
inline
    void MultipoleT::setEntranceAngle(double entranceAngle) {
        entranceAngle_m = entranceAngle;
}
inline
    double MultipoleT::getEntranceAngle() const {
        return entranceAngle_m;
}
inline 
    double MultipoleT::getTransProfile(int n) const {
        return transProfile_m[n];
}
inline
    std::vector<double> MultipoleT::getTransProfile() const {
        return transProfile_m;
}
inline 
    double MultipoleT::getDipoleConstant() const {
         return transProfile_m[0];
}
inline
383
    std::size_t MultipoleT::getMaxOrder() const {
ext-rogers_c's avatar
ext-rogers_c committed
384 385
         return maxOrder_m;
}
386

ext-rogers_c's avatar
ext-rogers_c committed
387
inline
388 389
    std::size_t MultipoleT::getMaxXOrder() const {
        return maxOrderX_m;
ext-rogers_c's avatar
ext-rogers_c committed
390
}
391

ext-rogers_c's avatar
ext-rogers_c committed
392
inline
393 394 395 396 397
    void MultipoleT::setMaxXOrder(std::size_t maxOrderX) {
        maxOrderX_m = maxOrderX;
}
inline
    std::size_t MultipoleT::getTransMaxOrder() const {
ext-rogers_c's avatar
ext-rogers_c committed
398 399 400
        return transMaxOrder_m;
}
inline 
401
    void MultipoleT::setTransMaxOrder(std::size_t transMaxOrder) {
ext-rogers_c's avatar
ext-rogers_c committed
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
        transMaxOrder_m = transMaxOrder;
	transProfile_m.resize(transMaxOrder + 1, 0.);
}
inline
    double MultipoleT::getRotation() const {
         return rotation_m;
}
inline 
    void MultipoleT::setRotation(double rot) {
         rotation_m = rot;
}
inline
    void MultipoleT::setBendAngle(double angle) {
        angle_m = angle;
}
inline
    double MultipoleT::getBendAngle() const {
        return angle_m;
}
inline
    void MultipoleT::setLength(double length) {
        length_m = std::abs(length);
}
inline
    double MultipoleT::getLength() const {
        return length_m;
}
429 430 431 432 433 434 435 436
inline
    double MultipoleT::getBoundingBoxLength() const {
        return boundingBoxLength_m;
}
inline
    void MultipoleT::setBoundingBoxLength(const double &boundingBoxLength) {
        boundingBoxLength_m = boundingBoxLength;
}
ext-rogers_c's avatar
ext-rogers_c committed
437 438

#endif