OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
VerticalFFAMagnet.cpp
Go to the documentation of this file.
1 //
2 // Class VerticalFFAMagnet
3 // Defines the abstract interface for a vertical FFA magnet
4 // with vertical scaling fringe fields.
5 //
6 // Copyright (c) 2019 - 2023, Chris Rogers, STFC Rutherford Appleton Laboratory, Didcot, UK
7 // All rights reserved
8 //
9 // This file is part of OPAL.
10 //
11 // OPAL is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
18 //
20 
23 
24 #include <cmath>
25 
27  : Component(name), straightGeometry_m(1.) {
28 }
29 
31  Component(right),
32  straightGeometry_m(right.straightGeometry_m),
33  dummy(right.dummy),
34  maxOrder_m(right.maxOrder_m),
35  k_m(right.k_m),
36  Bz_m(right.Bz_m),
37  zNegExtent_m(right.zNegExtent_m),
38  zPosExtent_m(right.zPosExtent_m),
39  halfWidth_m(right.halfWidth_m),
40  bbLength_m(right.bbLength_m),
41  endField_m(right.endField_m->clone()),
42  dfCoefficients_m(right.dfCoefficients_m)
43 {
45 }
46 
47 
49 }
50 
53  magnet->initialise();
54  return magnet;
55 }
56 
58  return dummy;
59 }
60 
62  return dummy;
63 }
64 
68 }
69 
71  double& /*startField*/, double& /*endField*/) {
72  RefPartBunch_m = bunch;
73  initialise();
74 }
75 
77  RefPartBunch_m = nullptr;
78 }
79 
81  return straightGeometry_m;
82 }
83 
85  return straightGeometry_m;
86 }
87 
89  visitor.visitVerticalFFAMagnet(*this);
90 }
91 
93  if (std::abs(R[0]) > halfWidth_m ||
94  R[2] < 0. || R[2] > bbLength_m ||
95  R[1] < -zNegExtent_m || R[1] > zPosExtent_m) {
96  return true;
97  }
98  std::vector<double> fringeDerivatives(maxOrder_m+2, 0.);
99  double zRel = R[2]-bbLength_m/2.; // z relative to centre of magnet
100  for (size_t i = 0; i < fringeDerivatives.size(); ++i) {
101  fringeDerivatives[i] = endField_m->function(zRel, i); // d^i_phi f
102  }
103 
104  std::vector<double> x_n(maxOrder_m+1); // x^n
105  x_n[0] = 1.; // x^0
106  for (size_t i = 1; i < x_n.size(); ++i) {
107  x_n[i] = x_n[i-1] * R[0];
108  }
109 
110  // note that the last element is always 0, because dfCoefficients_m is
111  // of size maxOrder_m+1. This leads to better Maxwellianness in testing.
112  std::vector<double> f_n(maxOrder_m + 2, 0.);
113  std::vector<double> dz_f_n(maxOrder_m + 1, 0.);
114  for (size_t n = 0; n < dfCoefficients_m.size(); ++n) {
115  const std::vector<double>& coefficients = dfCoefficients_m[n];
116  for (size_t i = 0; i < coefficients.size(); ++i) {
117  f_n[n] += coefficients[i] * fringeDerivatives[i];
118  dz_f_n[n] += coefficients[i] * fringeDerivatives[i+1];
119  }
120  }
121  double bref = Bz_m * std::exp(k_m * R[1]);
122  B[0] = 0.;
123  B[1] = 0.;
124  B[2] = 0.;
125  for (size_t n = 0; n < x_n.size(); ++n) {
126  B[0] += bref * f_n[n+1] * (n+1) / k_m * x_n[n];
127  B[1] += bref * f_n[n] * x_n[n];
128  B[2] += bref * dz_f_n[n] / k_m * x_n[n];
129  }
130  return false;
131 }
132 
134  dfCoefficients_m = std::vector< std::vector<double> >(maxOrder_m+1);
135  dfCoefficients_m[0] = std::vector<double>(1, 1.);
136  if (maxOrder_m > 0) {
137  dfCoefficients_m[1] = std::vector<double>();
138  }
139  // n indexes like the polynomial order of the midplane expansion
140  // e.g. Bz = exp(mz) f_n y^n
141  // where y is distance from the midplane and z is height
142  for (size_t n = 2; n < dfCoefficients_m.size(); n+=2) {
143  const std::vector<double>& oldCoefficients = dfCoefficients_m[n-2];
144  std::vector<double> coefficients(oldCoefficients.size() + 2, 0);
145  // j indexes the derivative of f_0
146  for (size_t j = 0; j < oldCoefficients.size(); ++j) {
147  coefficients[j] += -1./(n)/(n-1) * k_m * k_m * oldCoefficients[j];
148  coefficients[j+2] += -1./(n)/(n-1) * oldCoefficients[j];
149  }
150  dfCoefficients_m[n] = coefficients;
151  }
152 }
153 
155  endField_m.reset(endField);
156  endField_m->setMaximumDerivative(maxOrder_m);
157 }
158 
159 void VerticalFFAMagnet::setMaxOrder(size_t maxOrder) {
160  if (endField_m.get()) {
161  endField_m->setMaximumDerivative(maxOrder);
162  }
163  maxOrder_m = maxOrder;
164 }
ElementBase * clone() const
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Abstract base class for electromagnetic fields.
Definition: EMField.h:188
void setEndField(endfieldmodel::EndFieldModel *endField)
void accept(BeamlineVisitor &visitor) const
VerticalFFAMagnet(const std::string &name)
StraightGeometry straightGeometry_m
Abstract base class for accelerator geometry classes.
Definition: Geometry.h:43
PartBunchBase< double, 3 > * RefPartBunch_m
Definition: Component.h:191
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &)=0
Apply the algorithm to a vertical FFA magnet.
Tps< T > exp(const Tps< T > &x)
Exponential.
Definition: TpsMath.h:165
void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField)
std::vector< std::vector< double > > dfCoefficients_m
const std::string name
BGeometryBase & getGeometry()
BMultipoleField dummy
std::unique_ptr< endfieldmodel::EndFieldModel > endField_m
Interface for a single beam element.
Definition: Component.h:50
void setMaxOrder(size_t maxOrder)
bool getFieldValue(const Vector_t &R, Vector_t &B) const
virtual void setElementLength(double length)
Set design length.
item[EANGLE] Entrance edge counterclockwise This enables to obtain skew at each point along the magnet
Definition: multipole_t.tex:32