Commit 410a1c3e authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Merge in MultipoleT routines; checked code compiles and passes tests

parents b13f3fc9 0d8ee91b
......@@ -39,6 +39,7 @@
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/MultipoleT.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
......@@ -61,6 +62,7 @@
#include "BeamlineGeometry/Euclid3D.h"
#include "BeamlineGeometry/PlanarArcGeometry.h"
#include "BeamlineGeometry/RBendGeometry.h"
#include "BeamlineGeometry/StraightGeometry.h"
#include "Beamlines/Beamline.h"
#include "Fields/BMultipoleField.h"
......@@ -98,6 +100,7 @@ class PartData;
using Physics::pi;
using Physics::q_e;
const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
......@@ -725,6 +728,22 @@ void ParallelCyclotronTracker::visitMultipole(const Multipole &mult) {
myElements.push_back(dynamic_cast<Multipole *>(mult.clone()));
}
/**
*
*
* @param multT
*/
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
*gmsg << "Adding MultipoleT" << endl;
if (opalRing_m != NULL) {
opalRing_m->appendElement(multT);
} else {
throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
"Need to define a RINGDEFINITION to use MultipoleT element");
}
myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}
/**
*
*
......
......@@ -121,6 +121,9 @@ public:
/// Apply the algorithm to a Multipole.
virtual void visitMultipole(const Multipole &);
/// Apply the algorithm to a MultipoleT
virtual void visitMultipoleT (const MultipoleT &);
/// Apply the algorithm to a Offset.
virtual void visitOffset(const Offset &);
......
......@@ -45,6 +45,7 @@ class Lambertson;
class Marker;
class Monitor;
class Multipole;
class MultipoleT;
class Offset;
class Patch;
class Probe;
......@@ -145,6 +146,9 @@ public:
/// Apply the algorithm to a multipole.
virtual void visitMultipole(const Multipole &) = 0;
/// Apply the algorithm to an arbitrary straight Multipole.
virtual void visitMultipoleT(const MultipoleT &) = 0;
/// Apply the algorithm to a patch.
virtual void visitPatch(const Patch &) = 0;
......
......@@ -192,7 +192,7 @@ bool Bend::apply(const size_t &i,
const double &t,
Vector_t &E,
Vector_t &B) {
if(designRadius_m > 0.0) {
// Shift position to magnet frame.
......@@ -221,7 +221,6 @@ bool Bend::apply(const Vector_t &R,
if(designRadius_m > 0.0) {
Vector_t X = R;
Vector_t bField(0.0);
if (!calculateMapField(X, bField)) {
......@@ -232,7 +231,6 @@ bool Bend::apply(const Vector_t &R,
return true;
}
return false;
}
......@@ -256,7 +254,6 @@ bool Bend::applyToReferenceParticle(const Vector_t &R,
return true;
}
return false;
}
......@@ -454,15 +451,15 @@ Vector_t Bend::calcCentralField(const Vector_t &R,
double deltaX) {
Vector_t B(0, 0, 0);
// double nOverRho = fieldIndex_m / designRadius_m;
// double expFactor = exp(-nOverRho * deltaX);
// double bxBzFactor = expFactor * nOverRho * R(1);
// Vector_t rotationCenter(-designRadius_m, R(1), 0.0);
// double cosangle = dot(R - rotationCenter, Vector_t(1, 0, 0)) / euclidian_norm(R - rotationCenter);
//double nOverRho = fieldIndex_m / designRadius_m;
//double expFactor = exp(-nOverRho * deltaX);
//double bxBzFactor = expFactor * nOverRho * R(1);
//Vector_t rotationCenter(-designRadius_m, R(1), 0.0);
//double cosangle = dot(R - rotationCenter, Vector_t(1, 0, 0)) / euclidian_norm(R - rotationCenter);
// B(0) = -bxBzFactor * cosangle;
// B(1) = expFactor * (1.0 - pow(nOverRho * R(1), 2.0) / 2.0);
// B(2) = -bxBzFactor * sqrt(1 - std::pow(cosangle, 2));
//B(0) = -bxBzFactor * cosangle;
//B(1) = expFactor * (1.0 - pow(nOverRho * R(1), 2.0) / 2.0);
//B(2) = -bxBzFactor * sqrt(1 - std::pow(cosangle, 2));
B(1) = 1.0;
......@@ -484,19 +481,23 @@ Vector_t Bend::calcEntranceFringeField(const Vector_t &R,
double engeFunc = gsl_spline_eval(entryFieldValues_m[0], Rprime(2), entryFieldAccel_m);
double engeFuncDeriv = gsl_spline_eval(entryFieldValues_m[1], Rprime(2), entryFieldAccel_m);
double engeFuncSecDeriv = gsl_spline_eval(entryFieldValues_m[2], Rprime(2), entryFieldAccel_m);
// double nOverRho = fieldIndex_m / designRadius_m;
// double expFactor = exp(-nOverRho * deltaX);
// double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
// double bXEntrance = -nOverRho * expFactor * Rprime(1) * engeFunc;
// double bYEntrance = (expFactor * engeFunc *
// (1.0 - 0.5 * trigFactor * pow(Rprime(1), 2.0)));
// double bZEntrance = -expFactor * Rprime(1) * engeFuncDeriv;
//double nOverRho = fieldIndex_m / designRadius_m;
//double expFactor = exp(-nOverRho * deltaX);
//double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
//double bXEntrance = -nOverRho * expFactor * Rprime(1) * engeFunc;
//double bYEntrance = (expFactor * engeFunc *
// (1.0 - 0.5 * trigFactor * pow(Rprime(1), 2.0)));
//double bZEntrance = expFactor * Rprime(1) * engeFuncDeriv;
// B(1) = (engeFunc *
// (1.0 - 0.5 * engeFuncSecDerivNorm * pow(Rprime(1), 2.0)));
B(1) = (engeFunc - 0.5 * engeFuncSecDeriv * pow(Rprime(1), 2.0));
B(2) = engeFuncDeriv * Rprime(1);
// (1.0 - 0.5 * engeFuncSecDerivNorm * pow(Rprime(1), 2.0)));
B(1) = (engeFunc - 0.5 * engeFuncSecDeriv * pow(Rprime(1), 2.0));
B(2) = engeFuncDeriv * Rprime(1);
}
return toEntranceRegion.rotateFrom(B);
......@@ -519,18 +520,20 @@ Vector_t Bend::calcExitFringeField(const Vector_t &R,
double engeFunc = gsl_spline_eval(exitFieldValues_m[0], Rprime(2), exitFieldAccel_m);
double engeFuncDeriv = gsl_spline_eval(exitFieldValues_m[1], Rprime(2), exitFieldAccel_m);
double engeFuncSecDeriv = gsl_spline_eval(exitFieldValues_m[2], Rprime(2), exitFieldAccel_m);
// double nOverRho = fieldIndex_m / designRadius_m;
// double expFactor = exp(-nOverRho * deltaX);
// double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
// double bXExit = -nOverRho * expFactor * Rprime(1) * engeFunc;
// double bYExit = (expFactor * engeFunc *
// (1.0 - 0.5 * trigFactor * pow(Rprime(1), 2.0)));
// double bZExit = expFactor * Rprime(1) * engeFuncDeriv;
//double nOverRho = fieldIndex_m / designRadius_m;
//double expFactor = exp(-nOverRho * deltaX);
//double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
//double bXExit = -nOverRho * expFactor * Rprime(1) * engeFunc;
//double bYExit = (expFactor * engeFunc *
// (1.0 - 0.5 * trigFactor * pow(Rprime(1), 2.0)));
//double bZExit = expFactor * Rprime(1) * engeFuncDeriv;
//B(1) = (engeFunc *
// (1.0 - 0.5 * engeFuncSecDerivNorm * pow(Rprime(1), 2.0)));
B(1) = (engeFunc - 0.5 * engeFuncSecDeriv * pow(Rprime(1), 2.0));
// B(1) = (engeFunc *
// (1.0 - 0.5 * engeFuncSecDerivNorm * pow(Rprime(1), 2.0)));
B(1) = (engeFunc - 0.5 * engeFuncSecDeriv * pow(Rprime(1), 2.0));
B(2) = engeFuncDeriv * Rprime(1);
}
return toExitRegion.rotateFrom(B);
......@@ -539,11 +542,10 @@ Vector_t Bend::calcExitFringeField(const Vector_t &R,
bool Bend::calculateMapField(const Vector_t &R, Vector_t &B) {
B = Vector_t(0.0);
bool verticallyInside = (std::abs(R(1)) < 0.5 * gap_m);
bool horizontallyInside = false;
Vector_t rotationCenter(-designRadius_m * cosEntranceAngle_m, R(1), designRadius_m * sinEntranceAngle_m);
if(inMagnetCentralRegion(R)) {
if (verticallyInside) {
double deltaX = 0.0;//euclidian_norm(R - rotationCenter) - designRadius_m;
......@@ -562,14 +564,14 @@ bool Bend::calculateMapField(const Vector_t &R, Vector_t &B) {
Vector_t BEntrance(0.0), BExit(0.0);
verticallyInside = (std::abs(R(1)) < gap_m);
if (inMagnetEntranceRegion(R)) {
horizontallyInside = true;
if (verticallyInside) {
BEntrance = calcEntranceFringeField(R, R(0));
}
}
if (inMagnetExitRegion(R)) {
horizontallyInside = true;
if (verticallyInside) {
......@@ -859,7 +861,8 @@ bool Bend::findIdealBendParameters(double chordLength) {
rotationZAxis_m += Physics::pi;
}
designRadius_m = chordLength / (2.0 * std::sin(angle_m / 2.0));
fieldAmplitude_m = ((refCharge / std::abs(refCharge)) *
fieldAmplitude_m = ((refCharge / std::abs(refCharge)) *
refBetaGamma * refMass /
(Physics::c * designRadius_m));
reinitialize = true;
......@@ -879,7 +882,6 @@ bool Bend::findIdealBendParameters(double chordLength) {
reinitialize = false;
}
return reinitialize;
}
......@@ -1186,6 +1188,7 @@ void Bend::setBendStrength() {
fieldAmplitude_m = ((charge / std::abs(charge)) * betaGamma * mass /
(Physics::c * designRadius_m));
// Search for angle if initial guess is not good enough.
findBendStrength(mass, gamma, betaGamma, charge);
}
......@@ -1776,4 +1779,4 @@ bool Bend::isInside(const Vector_t &R) const {
}
return (std::abs(R(1)) < 0.5 * gap_m && inMagnetCentralRegion(R));
}
\ No newline at end of file
}
......@@ -28,6 +28,8 @@
#include "gsl/gsl_spline.h"
#include "gsl/gsl_interp.h"
#include <gtest/gtest_prod.h>
#include <vector>
class Fieldmap;
......@@ -135,6 +137,8 @@ protected:
private:
FRIEND_TEST(Maxwell, Zeros);
// Not implemented.
void operator=(const Bend &);
......@@ -379,4 +383,4 @@ CoordinateSystemTrafo Bend::getBeginToEnd_local() const
return beginToEnd_lcs_m;
}
#endif // CLASSIC_BEND_H
\ No newline at end of file
#endif // CLASSIC_BEND_H
......@@ -22,6 +22,7 @@ set (_SRCS
Marker.cpp
Monitor.cpp
Multipole.cpp
MultipoleT.cpp
Offset.cpp
ParallelPlate.cpp
Patch.cpp
......@@ -71,6 +72,7 @@ set (HDRS
Marker.h
Monitor.h
Multipole.h
MultipoleT.h
Offset.h
ParallelPlate.h
Patch.h
......
......@@ -167,6 +167,7 @@ public:
, MONITOR
, MPSPLITINTEGRATOR
, MULTIPOLE
, MULTIPOLET
, MULTIPOLEWRAPPER
, OFFSET
, PARALLELPLATE
......@@ -705,4 +706,4 @@ bool ElementBase::isElementPositionSet() const
{ return elemedgeSet_m; }
#endif // CLASSIC_ElementBase_HH
\ No newline at end of file
#endif // CLASSIC_ElementBase_HH
This diff is collapsed.
/*
* Copyright (c) 2017, Titus Dascalu
* 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
* $Author: titus
*
* ---------------------------------------------------------------------
*
* 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"
#include "gsl/gsl_sf.h"
#include <vector>
class MultipoleT: public Component {
public:
/** Constructor: &name -> user-defined name */
explicit MultipoleT(const std::string &name);
/** Copy constructor */
MultipoleT(const MultipoleT &right);
/** Destructor */
~MultipoleT();
/** Inheritable copy constructor */
ElementBase* clone() const;
/** Return a dummy (0.) field value */
EMField &getField();
/** Return a dummy (0.) field value */
const EMField &getField() const;
/** Not implemented */
void getDimensions(double &zBegin, double &zEnd) const;
/** Calculate the field at some arbitrary position
*
* \param R -> position in the local 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
* \returns true if particle is outside the field map, else false
*/
bool apply(const Vector_t &R, const Vector_t &P, const double &t,
Vector_t &E, Vector_t &B);
/** 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
*/
bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
/** Initialise the MultipoleT
*
* \param -> bunch the global bunch object
* \param -> startField not used
* \param -> endField not used
*/
void initialise(PartBunchBase<double, 3>*, double &startField, double &endField);
/** Finalise the MultipoleT - sets bunch to NULL */
void finalise();
/** Return true if dipole component not zero */
bool bends() const;
/** Return the cell geometry */
PlanarArcGeometry& getGeometry();
/** Return the cell geometry */
const PlanarArcGeometry& getGeometry() const;
/** Accept a beamline visitor */
void accept(BeamlineVisitor& visitor) const;
/** Get the dipole constant B_0 */
double getDipoleConstant() const;
/** Set the dipole constant B_0 */
void setDipoleConstant(double B0);
/** Get the number of terms used in calculation of field components
* max power of z in Bz is 2 * maxOrder
*/
unsigned int getMaxOrder() const;
/** Set the number of terms used in calculation of field components
* max power of z in Bz is 2 * maxOrder
*/
void setMaxOrder(unsigned int maxOrder);
/** Get the maximum order in the given transverse profile
* -> highest power in given mid-plane field expansion
*/
unsigned int getTransMaxOrder() const;
/** Set the maximum order in the given transverse profile
* -> highest power in given mid-plane field expansion
*/
void setTransMaxOrder(unsigned int transMaxOrder);
/** Set transverse profile T(x)
*
* T(x) = B_0 + B1 x + B2 x^2 + B3 x^3 + ...
* n -> the order of the term (d^n/dx^n) to be set
* dTn -> value of n-th derivative
*/
void setTransProfile(unsigned int n, double Bn);
/** Get transverse profile n-th term */
double getTransProfile(int n) const;
/** Get all term of transverse profile */
std::vector<double> getTransProfile() const;
/** Set fringe field model
*
* Tanh model used here
* @f[ 1/2 * \left [tanh \left( \frac{s + s_0}{\lambda_{left}} \right)
* - tanh \left( \frac{s - s_0}{\lambda_{right}} \right) \right] @f]
* s0 -> centre field length and
* @f$ \lambda_{left} @f$ -> left end field length
* @f$ \lambda_{right} @f$ -> right end field length
* max_index -> (default 10) used to set up for differentiation - cannot
* calculate higher differentials than exist in max_index (this is always
* set to be 1 + twice the maxOrder_m when setting the fringe field)
* return true after fringe field is set
*/
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;
/** Set the entrance angle */
void setEntranceAngle(double entranceAngle);
/** Get the entrance angle */
double getEntranceAngle() const;
/** Get the bending radius
* not used
* when needed radius is found from length_m / angle_m
*/
double getBendRadius() const;
/** Set the length of the magnet
* if straight-> actual length
* if curved -> arc length
*/
void setLength(double length);
/** Get the length of the magnet */
double getLength() const;
/** not used */
double getChordLength() const;
/** Set the aperture dimensions
* this element only supports a rectangular aperture
*/
void setAperture(double vertAp, double horizAp);
/** Get the aperture dimensions
* returns a vector of 2 doubles
*/
std::vector<double> getAperture() const;
/** Set the angle of rotation of the magnet around its axis
* -> enables to make skew components
*/
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;
/** Set the step used in changing coord systems for variable radius
* -> units used: mm
* -> default: 50
* -> has a considerable effect on tracking time
*/
void setVarStep(double step);
/** Get the step used in changing coord system for variable radius */
double getVarStep() const;
private:
MultipoleT operator=(const MultipoleT& rhs);
// End fields
endfieldmodel::Tanh fringeField_l; // left
endfieldmodel::Tanh fringeField_r; // right
/** Field expansion parameters */
// number of terms used in calculating field components
unsigned int maxOrder_m = 0;
// highest power in given mid-plane field
unsigned int transMaxOrder_m = 0;
std::vector<double> transProfile_m;
// Geometry
PlanarArcGeometry planarArcGeometry_m;
/** Rotate frame for skew elements
* consecutive rotations:
* 1st -> about central axis
* 2nd -> azymuthal rotation
*/
Vector_t rotateFrame(const Vector_t &R);
/** Inverse of the 1st rotation in rotateFrame() method
* -> used to rotate B field back to global coord system
*/
Vector_t rotateFrameInverse(Vector_t &B);
/** Transform to Frenet-Serret coordinates for sector magnets
* -> if the radius is variable, the coord system is rotated
* step by step along the reference trajectory