Commit 0d8ee91b authored by ext-rogers_c's avatar ext-rogers_c
Browse files

Add MultipoleT routines

parent 4768386e
\input{header}
\vinput{header}
\chapter{Elements}
\label{chp:element}
......@@ -1126,7 +1126,7 @@ label:MULTIPOLE, TYPE=string, APERTURE=real-vector,
\end{example}
\begin{kdescription}
\item[KN]
A real vector \seesec{anarray},
A real vector \seesec{anarray},
containing the normal multipole coefficients,
$K_n=\diffp[n]{B_y}{x}$
(default is $\SI{0}{\tesla\meter\tothe{-n}}$).
......@@ -1154,6 +1154,71 @@ A multipole has no effect on the reference orbit, i.e. the reference system at i
%it has the same effect on the reference orbit as a \keyword{SBEND}
%with the same length and deflection angle \texttt{KN*L}.
\clearpage
\section{General Multipole (extension)}
\label{sec:multipoleT}
\index{MULTIPOLET}
A \keyword{MULTIPOLET} is in \opalt a general multipole with extended features. It can represent a straight or curved magnet. In the curved case, the user may choose between constant or variable radius. This model includes fringe fields.
\begin{example}
label:MULTIPOLET, L=real, ANGLE=real, VAPERT=real, HAPERT=real,
LFRINGE=real, RFRINGE=real, TP=real-vector, VARRADIUS=bool;
\end{example}
\begin{kdescription}
\item[L]
Physical length of the magnet (meters), without end fields. (Default: 1 m)
\item[ANGLE]
Physical angle of the magnet (radians). If not specified, the magnet is considered to be straight (ANGLE=0.0). This is not the total bending angle since the end fields cause additional bending. The radius of the multipole is set from the LENGTH and ANGLE attributes.
\item[VAPERT]
Vertical (non-bend plane) aperture of the magnet (meters). (Default: 0.5 m)
\item[HAPERT]
Horizontal (bend plane) aperture of the magnet (meters). (Default: 0.5 m)
\item[LFRINGE]
Length of the left fringe field (meters). (Default: 0.0 m)
\item[RFRINGE]
Length of the right fringe field (meters). (Default: 0.0 m)
\item[TP]
A real vector \seesec{anarray}, containing the multipole coefficients of the field expansion on the mid-plane in the body of the magnet: the transverse profile $ T(x) = B_0 + B_1 x + B_2 x^2 + \dots $ is set by TP={$B_0$, $B_1$, $B_2$} (units: $ T \cdot m^{-n}$). The order of highest multipole component is arbitrary, but all components up to the maximum must be given, even if they are zero.
\item[MAXFORDER]
The order of the maximum function $f_n$ used in the field expansion (default: 5). See the scalar magnetic potential below. This sets for example the maximum power of $z$ in the field expansion of vertical component $B_z$ to $2 \cdot \text{MAXFORDER} $.
\item[EANGLE]
Entrance edge angle (radians).
\item[ROTATION]
Rotation of the magnet about its central axis (radians, counterclockwise). This enables to obtain skew fields. (Default 0.0 rad)
\item[VARRADIUS]
This is to be set TRUE if the magnet has variable radius. More precisely, at each point along the magnet, its radius is computed such that the reference trajectory always remains in the centre of the magnet. In the body of the magnet the radius is set from the LENGTH and ANGLE attributes. It is then continuously changed to be proportional to the dipole field on the reference trajectory while entering the end fields. This attribute is only to be set TRUE for a non-zero dipole component. (Default: FALSE)
\item[VARSTEP]
The step size (meters) used in calculating the reference trajectory for VARRARDIUS = TRUE. It specifies how often the radius of curvature is re-calculated. This has a considerable effect on tracking time. (Default: 0.1 m)
\end{kdescription}
Superposition of many multipole components is permitted.
The reference system for a multipole is a Cartesian coordinate system for straight
geometry and a $(x,s,z)$ Frenet-Serret coordinate system for curved geometry. In the latter case, the axis $\hat{s}$ is the central axis of the magnet.
\ifthenelse{\boolean{ShowMap}}{\seefig{straight}}{}.
\noindent The following example shows a combined function magnet with a dipole component
of 2 Tesla and a quadrupole gradient of 0.1 Tesla/m.
\begin{example}
M30:MULTIPOLET, L=1, RFRINGE=0.3, LFRINGE=0.2, ANGLE=PI/6, TP=$\left\{ 2.0, 0.1 \right\}$, VARRADIUS=TRUE;
\end{example}
The field expansion used in this model is based on the following scalar potential:
\begin{equation}
V = z f_0(x,s) + \frac{z^3}{3!} f_1(x,s) + \frac{z^5}{5!} f_2(x,s) + \dots
\end{equation}
Mid-plane symmetry is assumed and the vertical component of the field on the mid-plane is given by the user under the form of the transverse profile $T(x)$. The full expression for the vertical component is then
\begin{equation}
B_z = f_0 = T(x) \cdot S(s)
\end{equation}
where $S(s)$ is the fringe field. This element uses the Tanh model for the end fields, having only three parameters (the centre length $s_0$ and the fringe field lengths $\lambda_{left}$, $\lambda_{right}$):
\begin{equation}
S(s) = \frac{1}{2} \left[ tanh \left( \frac{s + s_0}{\lambda_{left}} \right) -
tanh \left( \frac{s - s_0}{\lambda_{right}} \right) \right]
\end{equation}
Starting from Maxwell's laws, the functions $f_n$ are computed recursively and finally each component of the magnetic field is obtained from $V$ using the corresponding geometries.
%\clearpage
\clearpage
\section{Solenoid}
......@@ -2018,4 +2083,4 @@ KHV:KICKER, HKICK=0.001, VKICK=0.0005;
The reference system for an orbit corrector is a Cartesian coordinate system.
\index{Elements|)}
\input{footer}
\ No newline at end of file
\input{footer}
......@@ -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"
......@@ -60,6 +61,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"
......@@ -689,6 +691,23 @@ 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()));
}
/**
*
*
......@@ -2127,6 +2146,7 @@ bool ParallelCyclotronTracker::getFieldsAtPoint(const double &t, const size_t &P
bool outOfBound = (((*sindex)->second).second)->apply(Pindex, t, Efield, Bfield);
Bfield *= 0.10; // kGauss --> T
Efield *= 1.0e6; // kV/mm --> V/m
......
......@@ -108,6 +108,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;
......@@ -144,6 +145,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;
......@@ -240,4 +244,4 @@ void BeamlineVisitor::visitRBend3D(const RBend3D &) {
}
#endif // CLASSIC_BeamlineVisitor_HH
\ No newline at end of file
#endif // CLASSIC_BeamlineVisitor_HH
......@@ -193,7 +193,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.
......@@ -222,7 +222,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)) {
......@@ -233,7 +232,6 @@ bool Bend::apply(const Vector_t &R,
return true;
}
return false;
}
......@@ -257,7 +255,6 @@ bool Bend::applyToReferenceParticle(const Vector_t &R,
return true;
}
return false;
}
......@@ -455,15 +452,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;
......@@ -485,19 +482,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);
......@@ -520,18 +521,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);
......@@ -540,11 +543,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;
......@@ -563,14 +565,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) {
......@@ -860,7 +862,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;
......@@ -880,7 +883,6 @@ bool Bend::findIdealBendParameters(double chordLength) {
reinitialize = false;
}
return reinitialize;
}
......@@ -1187,6 +1189,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);
}
......@@ -1777,4 +1780,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
......@@ -20,6 +20,7 @@ set (_SRCS
Marker.cpp
Monitor.cpp
Multipole.cpp
MultipoleT.cpp
Offset.cpp
ParallelPlate.cpp
Patch.cpp
......@@ -44,6 +45,8 @@ include_directories (
${CMAKE_CURRENT_SOURCE_DIR}
)
add_subdirectory(EndFieldModel)
ADD_OPAL_SOURCES(${_SRCS})
set (HDRS
......@@ -68,6 +71,7 @@ set (HDRS
Marker.h
Monitor.h
Multipole.h
MultipoleT.h
Offset.h
ParallelPlate.h
Patch.h
......@@ -90,4 +94,4 @@ set (HDRS
VariableRFCavity.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
\ No newline at end of file
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
......@@ -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
/*
* 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.
*/
#include "MultipoleT.h"
using namespace endfieldmodel;
MultipoleT::MultipoleT(const std::string &name):
Component(name),
length_m(1.0),
entranceAngle_m(0.0),
planarArcGeometry_m(1., 1.),
angle_m(0.0),
designRadius_m(0.0),
maxOrder_m(5),
rotation_m(0.0),
verticalApert_m(0.5),
horizApert_m(0.5),
variableRadius_m(false),
fringeField_l(endfieldmodel::Tanh()),
fringeField_r(endfieldmodel::Tanh()),
varStep_m(0.1),
transMaxOrder_m(0) {
}
MultipoleT::MultipoleT(const MultipoleT &right):
Component(right),
maxOrder_m(right.maxOrder_m),
transProfile_m(right.transProfile_m),
transMaxOrder_m(right.transMaxOrder_m),
entranceAngle_m(right.entranceAngle_m),
angle_m(right.angle_m),
designRadius_m(right.designRadius_m),
length_m(right.length_m),
planarArcGeometry_m(right.planarArcGeometry_m),
fringeField_l(right.fringeField_l),
fringeField_r(right.fringeField_r),
verticalApert_m(right.verticalApert_m),
horizApert_m(right.horizApert_m),
rotation_m(right.rotation_m),
variableRadius_m(right.variableRadius_m),
varStep_m(right.varStep_m),
dummy() {
RefPartBunch_m = right.RefPartBunch_m;
}
MultipoleT::~MultipoleT() {
}
ElementBase* MultipoleT::clone() const {
return new MultipoleT(*this);
}
void MultipoleT::finalise() {
RefPartBunch_m = NULL;
}
bool MultipoleT::apply(const Vector_t &R, const Vector_t &P,
const double &t,Vector_t &E, Vector_t &B) {
/** Rotate coordinates around the central axis of the magnet */
Vector_t R_prime = rotateFrame(R);
/** If magnet is not straight go to coords along the magnet */
if(angle_m != 0.0) {
Vector_t X = R_prime;
R_prime = transformCoords(X);
}
if (insideAperture(R_prime)) {
B[0] = getBx(R_prime);
B[1] = getBz(R_prime);
B[2] = getBs(R_prime);
return false;
}
else {
for(int i = 0; i < 3; i++) {
B[i] = 0.0;
}
return true;
}
}
bool MultipoleT::apply(const size_t &i, const double &t,
Vector_t &E, Vector_t &B) {
return apply(RefPartBunch_m->R[i], RefPartBunch_m->P[i], t, E, B);
}
Vector_t MultipoleT::rotateFrame(const Vector_t &R) {
Vector_t R_prime(3), R_pprime(3);
/** Apply two 2D rotation matrices to coords vector
* -> rotate around central axis => skew fields
* -> rotate azymuthaly => entrance angle
*/
// 1st rotation
R_prime[0] = R[0] * cos(rotation_m) + R[1] * sin(rotation_m);
R_prime[1] = - R[0] * sin(rotation_m) + R[1] * cos(rotation_m);
R_prime[2] = R[2];
// 2nd rotation
R_pprime[0] = R_prime[2] * sin(entranceAngle_m) +
R_prime[0] * cos(entranceAngle_m);
R_pprime[1] = R_prime[1];
R_pprime[2] = R_prime[2] * cos(entranceAngle_m) -
R_prime[0] * sin(entranceAngle_m);
return R_pprime;
}
Vector_t MultipoleT::rotateFrameInverse(Vector_t &B) {
/** This function represents the inverse of the rotation
* around the central axis performed by rotateFrame() method
* -> used to rotate B field back to global coord system
*/
Vector_t B_prime(3);
B_prime[0] = B[0] * cos(rotation_m) + B[1] * sin(rotation_m);
B_prime[1] = - B[0] * sin(rotation_m) + B[1] * cos(rotation_m);
B_prime[2] = B[2];
return B_prime;
}
bool MultipoleT::insideAperture(const Vector_t &R) {
if (abs(R[1]) <= verticalApert_m / 2. && abs(R[0]) <= horizApert_m / 2.) {
return true;
}
else {
return false;
}
}
Vector_t MultipoleT::transformCoords(const Vector_t &R) {
Vector_t X(3), Y(3);
// if radius not variable
if (not variableRadius_m) {
double radius = length_m / angle_m;
// Transform from Cartesian to Frenet-Seret along the magnet
double alpha = atan(R[2] / (R[0] + radius ));
if (alpha != 0.0) {
X[0] = R[2] / sin(alpha) - radius;
X[1] = R[1];
X[2] = radius * alpha;
}
else {
X = R;
}
}
// if radius is variable need to transform coordinates at each
// point along the trajectory
else {
double deltaAlpha, S = 0.0, localRadius;
double stepSize = varStep_m; // mm -> has a big effect on tracking time
if (abs(R[2]) <= stepSize) {
return R; // no transformation around origin
}
Y = R;
Vector_t temp;
while (abs(Y[2]) > stepSize and getRadius(S) != -1) {
localRadius = getRadius(S);
deltaAlpha = stepSize / localRadius;
if (R[2] < 0) {
deltaAlpha *= - 1.; // rotate in the other direction
}
temp = Y;
// translation
temp[0] += localRadius * (1 - cos(deltaAlpha));
temp[2] -= localRadius * sin(deltaAlpha);