Commit c480c1bf authored by Martin Duy Tat's avatar Martin Duy Tat Committed by Martin Tat
Browse files

Added classes for handling differential operators used in MultipoleT

MultipoleT now converges for curved magnets
MultipoleT is split into three derived classed and a base class
parent 6ce2ae4d
......@@ -40,6 +40,10 @@
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/MultipoleT.h"
#include "AbsBeamline/MultipoleTBase.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
......@@ -745,7 +749,6 @@ void ParallelCyclotronTracker::visitMonitor(const Monitor &corr) {
// applyDrift(flip_s * corr.getElementLength());
}
/**
*
*
......@@ -772,6 +775,54 @@ void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}
/**
*
*
* @param multTstraight
*/
void ParallelCyclotronTracker::visitMultipoleTStraight(const MultipoleTStraight &multTstraight) {
*gmsg << "Adding MultipoleTStraight" << endl;
if (opalRing_m != NULL) {
opalRing_m->appendElement(multTstraight);
} else {
throw OpalException("ParallelCyclotronTracker::visitMultipoleTStraight",
"Need to define a RINGDEFINITION to use MultipoleTStraight element");
}
myElements.push_back(dynamic_cast<MultipoleTStraight *>(multTstraight.clone()));
}
/**
*
*
* @param multTccurv
*/
void ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &multTccurv) {
*gmsg << "Adding MultipoleTCurvedConstRadius" << endl;
if (opalRing_m != NULL) {
opalRing_m->appendElement(multTccurv);
} else {
throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius",
"Need to define a RINGDEFINITION to use MultipoleTCurvedConstRadius element");
}
myElements.push_back(dynamic_cast<MultipoleTCurvedConstRadius *>(multTccurv.clone()));
}
/**
*
*
* @param multTvcurv
*/
void ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &multTvcurv) {
*gmsg << "Adding MultipoleTCurvedVarRadius" << endl;
if (opalRing_m != NULL) {
opalRing_m->appendElement(multTvcurv);
} else {
throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius",
"Need to define a RINGDEFINITION to use MultipoleTCurvedVarRadius element");
}
myElements.push_back(dynamic_cast<MultipoleTCurvedVarRadius *>(multTvcurv.clone()));
}
/**
*
*
......@@ -3589,4 +3640,4 @@ void ParallelCyclotronTracker::injectBunch_m(bool& flagTransition) {
// After this, numBunch_m is wrong but not needed anymore...
numBunch_m--;
}
}
\ No newline at end of file
}
......@@ -134,6 +134,15 @@ public:
/// Apply the algorithm to a MultipoleT
virtual void visitMultipoleT (const MultipoleT &);
/// Apply the algorithm to a MultipoleTStraight
virtual void visitMultipoleTStraight (const MultipoleTStraight &);
/// Apply the algorithm to a MultipoleTCurvedConstRadius
virtual void visitMultipoleTCurvedConstRadius (const MultipoleTCurvedConstRadius &);
/// Apply the algorithm to a MultipoleTCurvedVarRadius
virtual void visitMultipoleTCurvedVarRadius (const MultipoleTCurvedVarRadius &);
/// Apply the algorithm to a Offset.
virtual void visitOffset(const Offset &);
......
......@@ -42,6 +42,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/RBend3D.h"
......@@ -136,6 +137,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 Probe.
virtual void visitProbe(const Probe &);
......@@ -381,6 +385,10 @@ inline void ParallelTTracker::visitMultipole(const Multipole &mult) {
itsOpalBeamline_m.visit(mult, *this, itsBunch_m);
}
inline void ParallelTTracker::visitMultipoleT(const MultipoleT &mult) {
itsOpalBeamline_m.visit(mult, *this, itsBunch_m);
}
inline void ParallelTTracker::visitProbe(const Probe &prob) {
itsOpalBeamline_m.visit(prob, *this, itsBunch_m);
}
......@@ -562,4 +570,4 @@ void ParallelTTracker::kickParticlesDKS() {
}
#endif
#endif // OPAL_ParallelTTracker_HH
\ No newline at end of file
#endif // OPAL_ParallelTTracker_HH
......@@ -193,3 +193,7 @@ install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/src")
# cmake-tab-width: 4
# indent-tabs-mode:nil
# End:
#Turn on gcov
#set(CMAKE_CXX_FLAGS "-g -O0 -Wall -fprofile-arcs -ftest-coverage")
#set(CMAKE_CXX_OUTPUT_EXTENSION_REPLACE 1)
......@@ -47,6 +47,9 @@ class Marker;
class Monitor;
class Multipole;
class MultipoleT;
class MultipoleTStraight;
class MultipoleTCurvedConstRadius;
class MultipoleTCurvedVarRadius;
class Offset;
class Patch;
class Probe;
......@@ -150,9 +153,18 @@ public:
/// Apply the algorithm to a multipole.
virtual void visitMultipole(const Multipole &) = 0;
/// Apply the algorithm to an arbitrary straight Multipole.
/// Apply the algorithm to an arbitrary Multipole.
virtual void visitMultipoleT(const MultipoleT &) = 0;
/// Apply the algorithm to an arbitrary straight Multipole.
virtual void visitMultipoleTStraight(const MultipoleTStraight &) = 0;
/// Apply the algorithm to an arbitrary curved Multipole of constant radius.
virtual void visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &) = 0;
/// Apply the algorithm to an arbitrary curved Multipole of variable radius.
virtual void visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &) = 0;
/// Apply the algorithm to a patch.
virtual void visitPatch(const Patch &) = 0;
......@@ -252,4 +264,4 @@ void BeamlineVisitor::visitRBend3D(const RBend3D &) {
}
#endif // CLASSIC_BeamlineVisitor_HH
\ No newline at end of file
#endif // CLASSIC_BeamlineVisitor_HH
add_subdirectory(EndFieldModel)
add_subdirectory(MultipoleTFunctions)
set (_SRCS
AlignWrapper.cpp
......@@ -24,6 +25,10 @@ set (_SRCS
Monitor.cpp
Multipole.cpp
MultipoleT.cpp
MultipoleTBase.cpp
MultipoleTStraight.cpp
MultipoleTCurvedConstRadius.cpp
MultipoleTCurvedVarRadius.cpp
Offset.cpp
ParallelPlate.cpp
Patch.cpp
......@@ -75,6 +80,10 @@ set (HDRS
Monitor.h
Multipole.h
MultipoleT.h
MultipoleTBase.h
MultipoleTStraight.h
MultipoleTCurvedConstRadius.h
MultipoleTCurvedVarRadius.h
Offset.h
ParallelPlate.h
Patch.h
......@@ -98,4 +107,8 @@ 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")
#Turn on gcov
#set(CMAKE_CXX_FLAGS "-g -O0 -Wall -fprofile-arcs -ftest-coverage")
#set(CMAKE_CXX_OUTPUT_EXTENSION_REPLACE 1)
This diff is collapsed.
This diff is collapsed.
/*
* Copyright (c) 2017, Titus Dascalu
* Copyright (c) 2018, Martin Duy Tat
* 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 "MultipoleTBase.h"
#include <gsl/gsl_sf_gamma.h>
#include "AbsBeamline/MultipoleTFunctions/tanhDeriv.h"
using namespace endfieldmodel;
MultipoleTBase::MultipoleTBase():
Component(),
fringeField_l(endfieldmodel::Tanh()),
fringeField_r(endfieldmodel::Tanh()),
maxOrder_m(3),
transMaxOrder_m(0),
length_m(1.0),
entranceAngle_m(0.0),
rotation_m(0.0),
boundingBoxLength_m(0.0),
verticalApert_m(0.5),
horizontalApert_m(0.5) {
}
MultipoleTBase::MultipoleTBase(const std::string &name):
Component(name),
fringeField_l(endfieldmodel::Tanh()),
fringeField_r(endfieldmodel::Tanh()),
maxOrder_m(3),
transMaxOrder_m(0),
length_m(1.0),
entranceAngle_m(0.0),
rotation_m(0.0),
boundingBoxLength_m(0.0),
verticalApert_m(0.5),
horizontalApert_m(0.5) {
}
MultipoleTBase::MultipoleTBase(const MultipoleTBase &right):
Component(right),
fringeField_l(right.fringeField_l),
fringeField_r(right.fringeField_r),
maxOrder_m(right.maxOrder_m),
transMaxOrder_m(right.transMaxOrder_m),
transProfile_m(right.transProfile_m),
length_m(right.length_m),
entranceAngle_m(right.entranceAngle_m),
rotation_m(right.rotation_m),
boundingBoxLength_m(right.boundingBoxLength_m),
verticalApert_m(right.verticalApert_m),
horizontalApert_m(right.horizontalApert_m),
dummy() {
RefPartBunch_m = right.RefPartBunch_m;
}
MultipoleTBase::~MultipoleTBase() {
}
bool MultipoleTBase::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);
/** Go to local Frenet-Serret coordinates */
R_prime[2] *= -1; //OPAL uses a different sign convention...
transformCoords(R_prime);
if (insideAperture(R_prime)) {
/** Calculate B-field in the local Frenet-Serret frame */
B[0] = getBx(R_prime);
B[1] = getBz(R_prime);
B[2] = getBs(R_prime);
/** Transform B-field from local to lab coordinates */
transformBField(B, R_prime);
B[2] *= -1; //OPAL uses a different sign convention...
return false;
} else {
for(int i = 0; i < 3; i++) {
B[i] = 0.0;
}
return true;
}
}
Vector_t MultipoleTBase::rotateFrame(const Vector_t &R) {
Vector_t R_prime(3), R_pprime(3);
/** Apply two 2D rotation matrices to coordinate 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 MultipoleTBase::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 coordinate 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 MultipoleTBase::setFringeField(const double &s0,
const double &lambda_l,
const double &lambda_r) {
fringeField_l.Tanh::setLambda(lambda_l);
fringeField_l.Tanh::setX0(s0);
fringeField_l.Tanh::setTanhDiffIndices(2 * maxOrder_m + 1);
fringeField_r.Tanh::setLambda(lambda_r);
fringeField_r.Tanh::setX0(s0);
fringeField_r.Tanh::setTanhDiffIndices(2 * maxOrder_m + 1);
return true;
}
double MultipoleTBase::getBz(const Vector_t &R) {
if (fabs(getFringeDeriv(0, R[2])) < 1.0e-12) {
return 0.0;
}
double Bz = 0.0;
std::size_t n = getMaxOrder() + 1;
while (n != 0) {
n--;
Bz = Bz * R[1] * R[1] + getFn(n, R[0], R[2])
/ gsl_sf_fact(2 * n);
}
return Bz;
}
double MultipoleTBase::getBx(const Vector_t &R) {
if (fabs(getFringeDeriv(0, R[2])) < 1.0e-12) {
return 0.0;
}
double Bx = 0.0;
std::size_t n = getMaxOrder() + 1;
while (n != 0) {
n--;
Bx = Bx * R[1] * R[1] + getFnDerivX(n, R[0], R[2])
/ gsl_sf_fact(2 * n + 1);
}
Bx *= R[1];
return Bx;
}
double MultipoleTBase::getBs(const Vector_t &R) {
if (fabs(getFringeDeriv(0, R[2])) < 1.0e-12) {
return 0.0;
}
double Bs = 0.0;
std::size_t n = getMaxOrder() + 1;
while (n != 0) {
n--;
Bs = Bs * R[1] * R[1] + getFnDerivS(n, R[0], R[2])
/ gsl_sf_fact(2 * n + 1);
}
Bs *= R[1] / getScaleFactor(R[0], R[2]);
return Bs;
}
double MultipoleTBase::getFringeDeriv(const std::size_t &n, const double &s) {
if (n <= 10) {
return (fringeField_l.getTanh(s, n) - fringeField_r.getNegTanh(s, n))
/ 2;
} else {
return tanhderiv::integrate(s,
fringeField_l.Tanh::getX0(),
fringeField_l.Tanh::getLambda(),
fringeField_r.Tanh::getLambda(),
n);
}
}
double MultipoleTBase::getTransDeriv(const std::size_t &n, const double &x) {
double func = 0.0;
std::size_t transMaxOrder = getTransMaxOrder();
if (n > transMaxOrder) {
return func;
}
std::vector<double> temp = getTransProfile();
for(std::size_t i = 1; i <= n; i++) {
for(std::size_t j = 0; j <= transMaxOrder; j++) {
if (j <= transMaxOrder_m - i) {
temp[j] = temp[j + 1] * (j + 1);
} else {
temp[j] = 0.0;
}
}
}
std::size_t k = transMaxOrder - n + 1;
while (k != 0) {
k--;
func = func * x + temp[k];
}
return func;
}
double MultipoleTBase::getFnDerivX(const std::size_t &n,
const double &x,
const double &s) {
if (n == 0) {
return getTransDeriv(1, x) * getFringeDeriv(0, s);
}
double deriv = 0.0;
double stepSize = 1e-3;
deriv += 1. * getFn(n, x - 2. * stepSize, s);
deriv += -8. * getFn(n, x - stepSize, s);
deriv += 8. * getFn(n, x + stepSize, s);
deriv += -1. * getFn(n, x + 2. * stepSize, s);
deriv /= 12 * stepSize;
return deriv;
}
double MultipoleTBase::getFnDerivS(const std::size_t &n,
const double &x,
const double &s) {
if (n == 0) {
return getTransDeriv(0, x) * getFringeDeriv(1, s);
}
double deriv = 0.0;
double stepSize = 1e-3;
deriv += 1. * getFn(n, x, s - 2. * stepSize);
deriv += -8. * getFn(n, x, s - stepSize);
deriv += 8. * getFn(n, x, s + stepSize);
deriv += -1. * getFn(n, x, s + 2. * stepSize);
deriv /= 12 * stepSize;
return deriv;
}
/*
* Copyright (c) 2017, Titus Dascalu
* Copyright (c) 2018, Martin Duy Tat
* 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_MULTIPOLETBASE_H
#define CLASSIC_MULTIPOLETBASE_H
/** ---------------------------------------------------------------------
*
* MultipoleTBase is a base class for a straight or curved combined \n
* function magnet (up to arbitrary multipole component) with fringe fields.\n
*
* ---------------------------------------------------------------------
*
* Class category: AbsBeamline \n
* $Author: Titus Dascalu, Martin Duy Tat, Chris Rogers
*
* ---------------------------------------------------------------------
*
* 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 "Algorithms/PartBunch.h"
#include "AbsBeamline/EndFieldModel/Tanh.h"
#include "Fields/BMultipoleField.h"
#include "Algorithms/Vektor.h"
#include "AbsBeamline/Component.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include <vector>
class MultipoleTBase: public Component {
public:
/** Default constructor */
MultipoleTBase();
/** Constructor
* \param name -> User-defined name
*/
explicit MultipoleTBase(const std::string &name);
/** Copy constructor */
MultipoleTBase(const MultipoleTBase &right);
/** Destructor */
~MultipoleTBase();
/** Return a dummy field value */
EMField &getField();
/** Return a dummy field value */
const EMField &getField() const;
/** Calculate the field at some arbitrary position \n
* 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
*/
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
* 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
*/
bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
/** Initialise the MultipoleT
* \param bunch -> Bunch the global bunch object