CC BY
34
UDC 519.6, 537.8
B. Erdelyi, M. Berz, M. Lindemann
Вестник СПбГУ. Сер. 10. 2015. Вып. 4
DIFFERENTIAL ALGEBRA BASED MAGNETIC FIELD COMPUTATIONS AND ACCURATE FRINGE FIELD MAPS
Michigan State University, 567, Wilson Road, East Lansing, MI 48824, United States of America
For the purpose of precision studies of transfer maps of particle motion in complex magnetic fields, we develop a method for Differential Algebra based 3D field computation and multipole decomposition. It can be applied whenever a model of a magnet is given which consist of line wire currents, and the wires are utilized to represent both the coils and the iron parts via the so-called image current method. Such a model exists for most modern superconducting magnets and a large variety of others as well. It is stressed that it is the only practically possible way to extract the multipoles and its derivatives, and hence the transfer map of the particle motion, analytically to high order. We also study various related topics like aspects of computational complexity of the problem, Maxwellification of fields, importance of vanishing curl, etc., and its applications to very accurate computation of magnetic fields including fringe fields. Refs 19. Figs 4. Tables 6.
Keywords: Maxwell's equations, Biot—Savart law, multipoles, differential algebra.
Б. Эрдели, M. Берц, M. Линдеман
РАСЧЕТ МАГНИТНЫХ ПОЛЕЙ И ТОЧНЫХ МАТРИЦ ПЕРЕХОДА ДЛЯ КРАЕВЫХ ПОЛЕЙ С ПРИМЕНЕНИЕМ ДИФФЕРЕНЦИАЛЬНОЙ АЛГЕБРЫ
Мичиганский государственный университет, Соединенные Штаты Америки, 48824, Ист Ленсинг, Уилсон-роуд, 567
Для точного исследования матриц перехода для движения частиц в сложных магнитных полях разработан метод расчета трехмерных полей и декомпозиции мультиполей, основанный на дифференциальной алгебре. Он может применяться к любой модели магнита, состоящей из линейных проводов, и провода используются для представления катушек индуктивности и железных частей с помощью так называемого метода изображений. Такая модель существует для большинства как современных сверхпроводящих магнитов, так других. Следует подчеркнуть, что это практически единственный возможный путь извлечения мультиполей и их производных и, следовательно, аналитически матриц перехода высокого порядка для движения частиц. Также рассматриваются несколько смежных тем, таких как вычислительная сложность задачи, максвеллификация полей, важность исчезающих завитков и т. д., и их приложения к очень точным вычислениям магнитных полей, включая краевые поля. Библиогр. 19 назв. Ил. 4. Табл. 6.
Ключевые слова: уравнения Максвелла, закон Био—Савара, мультиполи, дифференциальная алгебра.
Introduction. The performance of modern high energy accelerators, such as the Large Hadron Collider at CERN depends critically on the field quality of superconducting or air coil magnets employed to guide and focus the circulating beam [1]. The nonlinearities
Erdelyi Bela — doctor of phylosophy; e-mail: berdelyi@niu.edu Berz Martin — doctor of phylosophy, professor; e-mail: berz@msu.edu
Lindemann Michael — doctor of phylosophy, professor; e-mail: lindemann@dhbw-loerrach.de
Эрдели Бела — доктор философии; e-mail: berdelyi@niu.edu
Берц Мартин — доктор философии, профессор; e-mail: berz@msu.edu
Линдеман Майкл — доктор философии, профессор; e-mail: lindemann@dhbw-loerrach.de
of these magnets drive resonances, potentially rendering the motion of particles at large amplitudes unstable [2]. The shrinkage of the useful region in space, called the dynamic aperture, due to magnet nonlinearities is detrimental to the stringent high luminosity requirements, so a careful design of these magnets is in order.
The design of air coil magnets is often performed by dedicated codes like ROXIE [3]. The accurate placement of the wires, followed by extensive optimization, produces an analytical model of the magnet. In a similar vein, many modern compact accelerators including certain FFAGs [4] as well as novel applications like those connected to planned precision measurement for electric dipole moments require very precise field arrangements [5].
Using the Biot-Savart law, the resulting magnetic field is computed on the surface of a cylinder coaxial with the optical axis, and subsequently Fourier analyzed numerically to reveal its multipole content. Several iterations are necessary to obtain a magnet model that satisfies the design specifications. However, once the magnet model is ready, the dynamical studies can use only the multipole data output of the magnet design codes, usually organized in tables. These contain the integrated values of the multipoles [6]. This is satisfactory for the interior of the multipoles where the fields are independent of the arc length s, which is used as the independent variable. However, it is not necessarily accurate enough for the end regions, the fringe fields, where the s dependence of the fields could result in unusual local dynamics not revealed by the integrated values.
In the following we develop a theory that solves this problem by computing the magnetic fields based on Differential Algebraic (DA) methods. This method allows not just the extraction of the multipole strengths, but also their full s dependence, allowing analytical computation of s derivatives, which are necessary for "exact" fringe field map computations. The first steps in this direction have been done by Caspi [7], which we here extend to a general theory. Subsequently we derive an improved, numerically stable version of the Biot-Savart law for straight line current wires in 3D and explain the principle of DA based field computation. This is followed by two methods for multipole extraction. The importance of enforcing Maxwell's equations is presented based on two methods to enforce Maxwell's equations, a local and a global approach. To conclude, various examples of the computation of multipoles and applications to map computation are given.
The Biot-Savart law and field computation. The magnetic field computation is based on the Biot-Savart law. To solve the equations of the motion, and hence to obtain the map of a magnetic element it is necessary to compute not just the value of the field at a certain point in space, but also its derivatives and hence its Taylor expansion. So, why do we want to use Differential Algebraic methods [8-10] to achieve this? In principle, it is possible to determine the derivatives analytically and implement it in some code to evaluate them. This was done using Mathematica, and the results are presented in Table 1. We mention that the calculations have been done in only one variable (x), for one field component (By), and one single line current. As it turns out, higher order computations are very resource intensive, and Mathematica usually ran out of memory at the computation of the order 10 derivative. Therefore, it is clear that for realistic magnet models, consisting of several 105 line currents, in at least 2 variables up to high orders this way is practically intractable. As a comparison, Table 2 shows the speed of the DA method on the same computer, which at the same time preserves the accuracy of the computed derivatives. Technically, to obtain the Taylor expansion of the field components around a specified point in space, all we need is the evaluation of the Biot-Savart law in DA.
Table 1. Results of transforming the analytic derivatives of the Biot—Savart law calculated with Mathematica to FORTRAN code
Order, n dnBy/dxn (0,0,0)
Generation Time Max. mem. Lines Evaluation time
0 0.02 891 4 0.3
1 0.09 912 20 0.8
2 0.33 1.026 69 2.4
3 0.94 1.306 188 7.0
4 2.50 1.932 457 16.0
5 6.26 3.213 1009 36.0
6 12.95 5.709 2078 72.9
7 27.23 10.275 4059 140.9
8 50.16 18.591 7567 267.0
9 96.73 33.132 13603 480.0
Note. Shown are the CPU time required to generate the code in sec, the required memory in kB, the length of the generated FORTRAN code, and the CPU time required for the evaluation of the derivative in msec.
Table 2. Evaluation time in sec of the Taylor expansion of the y component of the magnetic field in DA at various orders
Order, n dnBy/dxn (0,0,0)
Evaluation time
1 4.8-10~b
5 6.8-KT6
10 10.7-KT6
15 20.5-KT6
20 47.8-10-6
25 93.7-KT6
While the exact form of the formula is not critical for the evaluation of the magnetic field value at a certain point in the space as long it is mathematically accurate, it does have a significant influence when it is used to compute also the Taylor expansion of the field around a point. Some of the shortcomings of a naive implementation have been pointed out by Caspi [7]. Another numerical instability has been noticed by us when we utilized it at a point where some of the wire currents were exactly or almost collinear with the point of expansion. Therefore, several modifications of the standard form of the Biot-Savart law are needed to obtain a numerically stable version which has good behavior in all arising situations [11-13].
As a consequence of Ampere's law, the elementary magnetic flux density at a point r generated by a filamentary current wire dl situated at r' is given by the Biot-Savart formula
HoI dl x (r — r')
dB
4n
r — r'
To compute the magnetic field generated by an extended straight line current we parametrize the line by A G [0,1] and define r' (A) = rs + Al and re = rs + l, where rs, re represent the starting and end-point, respectively, of the line current to the point r. Integrating over the line
B = kl
dl x r'
~W
J |r.s + All 'J
+ Al|
dA
|rs + Al|3
3
1
1
with k = and introducing the shorthand notations a = |rs | , b = 2rs • l, c = |l| ,
the integral gives the result
dX
1
2b
2b
4c
(a + bX + cX2f/2 62 - 4acVv/« Va + b + c Va + b + ,
While mathematically accurate, this formula exhibits several severe numerical pitfalls that restrict its direct practical use, in particular when high-order derivatives are to be computed. Indeed, first the formula apparently exhibits a problem of cancellation of close-by numbers if b + c C a. Introduction of the quantity e = (b + c) /a yields
B
kl( lxr,)
y/a (52 - 4ac)
2b 1
1
4c
VlTeJ VI +£.
The first problem can be substantially alleviated now by observing that
1
1
Vl+e l+e + Vl+e
which yields the formula
B
kl( lxr,)
yfa (62 - 4ac)
2 be
4c
_l+e +Vl+e Vl+e.
(1)
However, there is a second numerical difficulty if the line current and the observation point are lying exactly or almost on the same line, because in this case b2 and 4ac assume similar values, which makes the evaluation of b2 — 4ac prone to numerical inaccuracy. To avoid this effect we rewrite the formula in terms of the angle 0 between l and rs. The relations among the angle and the products of vectors are
| sin d\ =
This implies the relationships
|1 x rs
HI • |r„
and cos 0 =
b2 - 4ac = -4 |rs|2 |l|2 sin2 6.
l-r„
|rs r
2be
4c
4 |rs||l| cos6 [2 |rs||l| cos6 + |l|2j 4 |i|2 |rs|
l+e+Vl+e Vl+e
|re| (|re| + |rs
|re
Finally, we obtain the magnetic field expressed in terms of rs and l as
B = —
kl (l X rs)
|rs|2 |rs + l| (|rs + l| + |rs|)
-|rs| +
|rs | cos2 6 + |l| cos 6 - |rs + l|
sin2 9
Denoting |rs| cos2 0 + |l| cos 0 = a and |rs + l| = ¡3, we manage to eliminate the sin2 0 term in the denominator with the help of the identity a — 3 = (a2 — 32) / (a + 3). Direct calculation shows that
a2 - f32 = - sin2 6 (jrs|2 cos2 6 + |rs + l|2) .
1
e
Altogether we obtain the final result kI (l X rs)
B
|rs|2 |rs + l| (|rs + l| + |rs|)
N +
! cos2 e + |rs + if
|rs | cos2 e + |i| cos e + |rs +1|
(2)
All expressions involve sums of positive numbers, so there is no cancellation problem. The only case where the expression is numerically unstable is when |rs | cos2 e +11| cos e + |rs + 1| approaches zero, that is e ^ n and |rs| < |1|. But this corresponds to a point in the close proximity of the wire where the field anyway diverges and is thus unavoidable.
Table 3. Taylor expansion of the magnetic field components Bx and By
r
s
Bx (x, s, y) By (x, s,y) Order Exp
0.9384473968456148E - 11 -.8339510041598585E - 11 0 0 0 0
0.1022945817408163E - 09 17.49010593690617 1 1 0 0
0.1430119891325347E - 09 -. 1304592905872117E - 09 1 0 1 0
17.49010593689972 0.9584448246993324E - 10 1 0 0 1
-.6067419448807509E - 09 0.5543819518033510E - 09 2 2 0 0
0.3602525017014357E - 08 -3.804751735973857 2 1 1 0
0.1205363555495387E - 08 —. 1106863867922359E - 08 2 0 2 0
-.2014572349218202E - 08 0.1941685518966341E - 08 2 1 0 1
-3.804751736131598 0.34442337840286 HE - 08 2 0 1 1
— . 5986214510200760E - 09 0.5524820548968856E - 09 2 0 0 2
— . 2282722315338770E - 07 -1.666257305421052 3 3 0 0
— . 9274588597119049E - 09 0.1244904457298190E - 09 3 2 1 0
0.7183921058029341E - 07 9.997543895896488 3 1 2 0
0.5557759530372408E - 09 — .2165498336204053E - 10 3 0 3 0
-4.998771979234908 — .3052480102017086E - 08 3 2 0 1
-.1133153901822226E - 06 0.1106802772765647E - 06 3 1 1 1
9.997543893641829 0.6954731157637895E - 07 3 0 2 1
— . 3357574107631933E - 08 -4.998771979536424 3 1 0 2
-.7398769930055948E - 09 — .5953132431457675E - 10 3 0 1 2
-1.666257304770050 — .2216494703055627E - 07 3 0 0 3
N o t e. The columns represent the expansion coefficients, the order in the expansion and the exponents of (x,s,y), respectively.
The necessary ingredients for the DA field calculation are the above formula and the magnet model, consisting of line wire currents. To this end, the entire field in space is calculated by summing up the fields created by wire currents. At each step, the evaluation of (2) in DA yields not only the value of the magnetic field at the respective point, but also its derivatives, that is the Taylor expansion with respect to coordinates.
The result for the return end of the Large Hadron Collider's High Gradient Quadrupole, up to order 3, is presented in Table 3. The correctness of our results have been checked against data obtained from Fermilab. The Fermilab data contains the values of the components of the field on the surface of a coaxial cylinder with the optical axis, and have been kindly supplied by G. Sabbi.
To show that the Biot-Savart law implementation based on (2) is much more stable then, for example, code based on (1), we use as an indicator the s component of the curl of the field, which, according to Maxwell's equations, should vanish. From Table 4 it is clear that the naive implementation goes wrong as low as second order in the curl. Due to space considerations we present the result only up to order 3, but our results show that the behavior of (2) is good up to very high orders. Also, we mention that probably this is the most straightforward and accurate way to compute the curl and hence verify whether Maxwell's equations are satisfied.
Table 4. Comparison of the s component of the curl, up to order 3, computed by two different implementation of the Biot—Savart law (eq. (2) and eq. (1) respectively)
(V X B)s from eq. (2) (V X B)s from eq. (1) Order Exp
0.6448175327022909E - 11 0.6451728040701710E - 11 0 0 0 0
0.3123336252824904E - 08 0.3123335699448115E - 08 1 1 0 0
0.157740931427952IE - 09 0.1577009633990656E - 09 1 0 1 0
0.3138928421006493E - 08 0.3138927988192986E - 08 1 0 0 1
0.6297175136893429E - 07 0.3368312531613524 2 2 0 0
0.1135643710736822E - 06 0.1135639355887008E - 06 2 1 1 0
0.2254658681977162E - 08 0.2237449336917052E - 08 2 0 2 0
0.6101880112296945E - 09 0.6101442684425235E - 09 2 1 0 1
0.1121600312625759E - 06 0.1121592163033647E - 06 2 0 1 1
-.6522627415961324E - 07 -.3368312554111546 2 0 0 2
-.1035163990081857E - 06 0.8482986704194671E - 06 3 3 0 0
0.3603849257016734E - 05 -1010523.362231316 3 2 1 0
0.2310457119847342E - 05 0.4430214737283222E - 06 3 1 2 0
0.236423638999561 IE - 07 — .3117654705420136E - 04 3 0 3 0
— .1977820552667708E - 05 — .3840225446083422E - 05 3 2 0 1
0.9434835135380232E - 08 0.9449170335074086E - 08 3 1 1 1
0.2256179271853398E - 05 0.6295506663533956E - 05 3 0 2 1
— .199990078186829 IE - 05 — .2090122526610116E - 05 3 1 0 2
— .3674740044701252E - 05 1010523.362276393 3 0 1 2
— .9278849333327333E - 07 0.5070636007076246E - 05 3 0 0 3
Multipole Extraction Algorithms. The Direct Method. Using the field computation of the preceding section, it is possible to extract the multipole content of magnetic fields directly, in a very elegant way that is arbitrary in order. In the following we assume straight optical axis. In the divergence-free, curl-free region of the magnets it is possible to derive the magnetic field components from a magnetic scalar potential that satisfies the Laplace equation. The general solution expanded in normal and skew components has the form
^
Vb (bkii(s)sml^ + ak,i(s) cos ¡4>) rk. (3)
k,l = 0
Defining 9k,i(s) and Mk,i(s) by
ak,i(s) V
we have an equivalent form
Vb =Y. Mk,i(s) cos (l^ + 0k,i(s)) rk,
k,i=o
which shows that any normal (skew) component can be obtained from the corresponding skew (normal) component by an s dependent rotation around the s axis. The link between the two forms in the other direction is
bk,i (s) = —Mk,i(s) sinQkj(s), ak,i (s) = Mk,i(s)cosQk,i (s).
Inserting (3) in the Laplace equation in cylindrical coordinates
Id (rdVB\ 1 d2VB d2VB
= ( —-— ] + — ^ ,o H—rr^- =
r dr \ dr J r2 d42 ds2
we obtain
Ë [{bk,i(s)(h2 - l2) + b'-2il(s))sm 14 + (ak,i(s)(k2 - l2) + a'^ (s)) cos l^rk-2 = 0, k,i=0
using the convention that the coefficients vanish for negative indices. Due to the fact that the above equation must hold for every r and 4, and the sin and cos are linearly independent
bKi(s)(k2 - l2) + b'k-2,i(s) = 0, ak,i(s)(h2 - l2) + a'k-2,i(s) = 0. Furthermore, it can be shown that the following recurrence relations hold
b(2in) (S) a(2in)(s) 6wW = n„=i(/2'_(/ + 2i/)2), al+2nJ(s) = n„=i(/2'_ {i + 2v)iy (5)
and the coefficients that cannot be obtained by these relations are zero. The bi,i(s), ai, i(s) are called normal and skew (true) multipoles respectively, while the terms that contain s derivatives are called pseudo-multipoles. It is worth mentioning that, as one can see from (4), the recurrence relations (5) do not hold in general for Mkil(s). However, when the 6kjl are s independent, (5) holds for Mkil(s). Inserting relations (5) in (3), we get for the potential
Vb = 22(fi(r,s)sml4 + gi(r, s) cos l4) rl
where
i=0
œ b(2n)(s)
2n __bi,i (s) 2n
Mr, ») - E - E n:=l{P-(i+2,)2)
r
n "
v=
= blAs)~Tn—TTr + ooTi—TT71—ô\r ~ ooAi,—Tv7—ôVn—o\r +•••'
4(l +1) 32(l +1)(l + 2) 384(l +1)(l + 2)(l + 3)
œ œ {2n)(„)
2n ai,i s 2
n=0 n=0
, , J'l,l (s) 2 al,l (S) 4__(S)_ 6,
ai'l[S) 4(/ + l)r 32(/ +!)(/ + 2)r 384(/+l)(/ + 2)(/ + 3) "' '
The functions fi(r, s) and gi(r, s) represent the "out of axis" expansion of the multipoles. The magnetic field components in cylindrical coordinates can be calculated using the well known formulas
dVB 1 dVB dVB
dr r dç ds
resulting in the expressions
œ
Br(r, 4, s) = go(r,s) + ^ fi(r,s)sinl4 + gi(r,s)cosl4ri-1, (6)
i=l
B0 (r,$,s) = ^[l (fi (r,s)cosl$ — gi(r,s)sinl$)] rl 1,
i=1 œ
Bs(r,$,s) = [fi (r,s)sinl$ + gl (r,s)cosl$
i=0
where prime denotes derivative with respect to s and
œ 2 œ (l + 2n)b( n (s)
fi(r, s) = V(/ + 2n)6;+2„ i(s)r = V — „„ ,';' ' r-" =
= (¿ + 2)6g(g) (/ + 4)bg(S) 4 (/ + 6)bg(S)
4(/+l) r 32(/+l)(/ + 2)r 384(/+l)(/ + 2)(/ + 3)
œ 2 œ (l + 2n)a(2,n'> (s) 2
Mr, s) = ¿2d + ^n)al+2n,l{s)r = £ ; ^^^ =
n=0 n=0 V=1 v v
g + 2)qg(g) (/ + 4)gff(S) 4 (/ + 6)gg(S)
4(/ + l) r 32(/+l)(/ + 2)r 384(/+l)(/ + 2)(/ + 3) "'
It can be seen that every multipole, except for l = 0, is multiplied by rl-1. For the special case l = 0, we get
œ j
k=1
B0(r,s) = 0,
œ
KM = -D-1)^1^^1^-)^.
k=0
In the DA picture, the field calculations are done locally, as a Taylor expansion of the field with respect to Cartesian coordinates x,y,s. Hence, we need the equations relating the cylindrical and Cartesian components of the magnetic field:
Bx(r,$,s) = Br (r,$,s)cos$ — B^(r,$, s)sm$, By (r,$,s) = Br (r,$,s)sin$ + B^(r,$,s)cos$,
and Bs (r, $, s) is unchanged. Clearly, if we evaluate the above equations in the midplane
(y = $ = 0), then
Br (r, $ = 0,s) \r^x = Bx (x, y = 0, s), Bj,(r, $ = 0,s) \r^x = By (x, y = 0, s).
So, all the information we need to extract the multipoles up to the order of calculation is in the Cartesian components of the fields in the midplane
Bx(x,y = 0,s)= g0(x,s)+^2 9i(x,s) ■ xl 1, By (x,y = 0,s) = ^ lfi (x,s) ■ xl 1.
=1 =1
This is possible due to the previously mentioned fact that any multipole strength of order l is multiplied by xl-1. Starting at l = 1, ai^s) is extracted as the x independent part
of Bx, and analogously 61,1(s) from By. Evaluating a1j1(s) and 61,1(s) at s = 0 yields the skew and normal dipole component. From a1j1(s) and 61,1(s) the functions f1(x,s), g1(x, s) are generated up to the order of calculation and subtracted from Bx(x, y = 0, s), and By (x,y = 0,s) respectively. This cancels the pseudo-multipoles generated by the s dependence of a1j1(s) and 61,1(s) (as in (5)), which otherwise would make the distinction between sextupole terms and pseudo-dipole terms impossible. The procedure can be iterated for the higher order multipoles, up to the order of calculation. After the k-th step, the remainder of the field components should contain just (k + l)-th and higher order multipoles.
However, there is an additional problem in the case of solenoidal fields (case l = 0). In this case we have an ao^s) in the potential, but its contribution vanishes from the field components Bx and By, so the function g0(x, s) cannot be generated from the information available in Bx and By. Fortunately, it can be generated from the Bs component, which evaluated at x = y = 0 yields a0 0(s). From this function we can calculate a0,0(s) up to a constant and generate the l = 0 contribution to Bx, g0(x, s). Once this is subtracted from Bx, the method works as previously described, starting with l = l.
Finally, two notes: the method relies on the fact that the magnetic field can be generated by a magnetic scalar potential that satisfies the Laplace equation. Therefore, it is really important that the curl of the field vanishes. If the fields are calculated from line currents by the Biot-Savart law, that means that the model should consist only of closed circuits to ensure vanishing curl. Maxwellification of the field ensures a better numerical stability of the algorithm. Secondly, only in the regions where the magnetic field is not s dependent the functions fl and gl are equal to the true multipoles, and lfl = fl, lgl = gl, an assumption that is sometimes made even for the s dependent region too.
Multipole Extraction by Analytical Fourier Transform. There is an alternate way to extract the multipoles. The field computation being performed in COSY in Cartesian coordinates (x,y,s), it is possible to perform an analytical pseudo-Fourier transform, meaning a series of coordinate transformation in DA, keeping the r and s dependence. This is solved by our S-Dependent Differential Algebraic Analytical Fourier Transform presented below.
Initially, the field components are in the form of (see eq. (10) on page 51). The first transformation is (x, y, s) ^ (r, cos $, sin $,s) by x = r cos $ and y = r sin $. At the same time, using (see eq. (11) on page 51) we switch to cylindrical coordinates, obtaining the field components in the form of (see eq. (12) on page 51). Note that we are going from a 3 variable representation to a 4 variable one, therefore some of the information in the new representation will turn out to be redundant for our purpose. The next transformation is (r, cos $, sin $,s) ^ (r,el$,e-l$,s), that is the exponential, and obviously a complex representation, using cos $ = (e+ e-l$) /2 and sin $ = (e— e-l$) / (2i).
It does not matter which field component is used for Fourier transformation, so assume that we are working with Br. Then, at this stage we have Br in the form of Br (r, ee-i$, s). Now, in principle, it is possible to recombine various products of powers of and to trigonometric functions involving multiple angles. However, the key point is to notice that we can obtain the true multipoles by setting e= 0. This is true due to the fact that all the terms of the form eiq$e-ip$, with q,p = 0, are responsible for the pseudo-multipoles. This becomes clear if one takes a closer look at (see eq. (13) on page 52). By setting eto zero we get rid of all the pseudo-multipole terms and we are going back to a 3 variable representation
n+1
Br(r,e-i*,s) = £ zl(s)e-il*rl-1, (7)
1=1
where zl(s) are complex functions of s. By comparison with (see eq. (14) on page 53) it is clear that one term of (7) can come only from [Al (s)cos ¡4 + Bl(s)sinl^j rl-1, which also can be expressed as
e-u, ( Bds) a4, (BM _
'22/ V 2 2
rl-1.
After setting eil* to zero and comparison with (7) we obtain the result Bi(s) = 2Re (zl(s)) and Al(s) = 2Im(zl(s)). As a final step, we take the true multipoles as given by
22
h,i(s) = yRe (zi(s)), aM(s) = ylm (z;(s)).
Enforcing Maxwell's Equations. Local Maxwellification. Given a magnet model consisting of line currents, it is possible to compute the magnetic field generated by the Biot-Savart law in the DA framework as a local Taylor expansion with respect to the Cartesian coordinates (x, y, s), as has been demonstrated in (2). The magnetic field should be divergence- and curl-free in the regions of interest, as implied by the Maxwell equations in a source-free region: V- B =0 and Vx B = 0. This is the case only when the magnet model consist of closed loops of current. Realistic magnet models, as for example the LHC HGQ end regions as modeled by the code ROXIE and supplied by G. Sabbi of FNAL, are not closed due to presence of image currents, "leads" and separate treatment of the two end regions (lead end and return end). One way to fix this problem is to input as much physical intuition as possible to close the magnet model, compute the field generated by this model, which should differ as little as possible from the original model. This is the first step of Maxwellification. Obviously, the solution is not unique due to infinitely many ways of closing the model. The closing is important to guarantee vanishing curl, as it is required by the Maxwell equations, and in this case the field is derivable from a scalar potential. From the DA computational point of view it is also important, because as seen above, it is enough to compute the field components only in the midplane. The computer time needed for computing a magnet model of several 105 line currents to high order in one end region of the LHC HGQ's scales much worse with the increase of the number of variables than with the increase of line currents, which should be linear.
Besides reduced computation time, a second step of the Maxwellification provides a way to correct for small computational errors or magnet model imperfections. One specific example is presented in the next section. Although the curl is already small, as shown in Table 4, the numerical stability of the multipole extraction algorithm is improved by local Maxwellification of the field, which is described below.
If we restrict ourselves to elements with straight optical axis for simplicity, the second step of Maxwellification proceeds as follows. Given Bx(x, 0, s), By(x, 0, s), Bs(x, 0, s) (y = 0 representing the midplane) we can compute the field components in the whole space in the following way: from a scalar potential V(x, y, s) that satisfies the Laplace equation
d2V(x,y, s) d2V{x,y,s) d2V(x,y, s) = Q dx2 dy2 ds2 '
the field results from the well-known relation B(x, y, s) = VV(x, y, s) (we neglect the sign
which is irrelevant in our discussion). We transform the Laplace equation to a fixed point problem by isolating the y derivative term
d2V(x, y , s) dy2
d2V(x, y, s) d2V(x, y, s)
dx2
ds2
and integrating twice with respect to y and considering the initial conditions in each of the steps, we obtain
y
V(x,y,s) = V(x, °,s) + J
dV(x,y", s)
dy"
dy' -
y''=0
y y
00
(d2V(x,y,s) d2V(x,y,s) \,n,,
{ dx2 + —JT2- )dy <*V,
so we obtain a fixed point problem for V(x,y,s). In the DA picture it converges to the exact solution in n steps, where n is the order of computation [14]. Since the Laplace equation is a second order PDE, we need two initial conditions. One is immediate from
dV (x, y" ,s)
dy"
dy' = yBy (x, 0, s),
because by definition dV/dy = By. This is already known, and the other initial condition to be calculated is the potential in the midplane V(x, 0, s).
In the ideal case the potential in the midplane is computed by a path integral, along an arbitrary path. This is the case when the curl of the initial field is exactly vanishing. Due to various causes previously mentioned this is almost never true. Nevertheless, the curl it is usually small. Then, one should use a path along which the field is deemed more accurately computed, yielding a potential, and subsequently field components that are close to the original, and curl that is almost vanishing. Hence the name Maxwellification.
Most of the time it is not obvious where the fields are computed more accurately. Then one could try different paths and choose the one giving the smallest curl. Convenient choices of paths are along the sides or diagonal of a rectangle in the midplane with opposite corners at (0, 0) and (x, s). In the midplane we have dV(r) = B(r) • dr, where r = (x, 0, s). Integrating, we get
V(r) = V(0) + J B (r") • dr'.
We can neglect the immaterial constant V(0), and integration along the sides in one direction gives
V (x-°s = ] B-{x",0•0)dx" + J B'(x 0s")ds".
00 Integration along the sides in the other direction gives
V (x-°- s) = I B-(°-°- '")ds" + i B-(x"-° s)"x"
y
y
s
x
For integration along the diagonal we set r' = Ar with A e [0,1]. Then, dr' = rdA and
V(x, 0, s) = x J Bx(Ax, 0,As)dA + s J Bs(Ax, 0,As)dA.
0 0
One can check by direct calculation that indeed Bx(x, 0,s) = dV(x, 0,s)/dx and Bs(x, 0, s) = dV(x, 0, s)/ds. This completes the Maxwellification procedure. Once we have V(x, y , s) we can compute the field components satisfying Maxwell's equations in the whole region of interest.
In DA the field components are computed as local Taylor expansions, so the method provides a local Maxwellification. That's why beside choosing the right path it might be useful to average over a certain region to decide which approach is the best. Finally, it should be obvious how to extend all the equations in the case of full 3D Maxwellification, if originally the field components are given in all 3 variables (x,y, s). For example, (8) is extended as
x s y
V xy,s) = J w^w + /«..O, + /By M,s)dy',
0 0 0
and in the same way in other cases. In this situation, of course, there are many more path choices and no fixed point transformation of the Laplace equation is needed.
Table 5. The s component of the curl, up to order n = 12, after Maxwellification of the field in Table 3
1
1
(V X B)s after Maxwellification n Exp
-.454747350886464 IE - 12 6 2 0 4
— .1455191522836685E - 10 7 2 1 4
0.3051757812500000E - 04 8 6 0 2
— .2328306436538696E - 09 8 4 2 2
— .3051757812500000E - 04 8 2 0 6
0.1862645149230957E - 08 9 4 3 2
— . 1862645149230957E - 08 10 6 2 2
— . 1490116119384766E - 07 10 4 4 2
— . 1490116119384766E - 07 10 2 4 4
— .1455191522836685E - 10 10 2 0 8
— . 284217094304040IE - 13 11 5 0 6
0.1490116119384766E - 07 11 4 1 6
— .355271367880050 IE - 14 11 2 0 9
0.5960464477539063E - 07 12 2 2 8
As an example, in Table 5 we present the s component of the curl of the Maxwellified field of Table 3. We mention that now we present the results up to order 12 in the curl and the first non-vanishing element occurs at order 6. Also, notice the improvement in comparison with the curl in Table 4.
Global Maxwellification. We saw that local Maxwellification is possible based on a closed magnet model consisting of line currents. One might imagine cases when the magnet model is not closed, and for some reason it is practically impossible to close it, or the actual closings change the original fields significantly. For such cases the S-Dependent Differential Algebraic Analytical Fourier Transform (SDDAAFT) provides a way for global Maxwellification and minimal modification of the original fields in a neighborhood of the optical axis. The only drawback compared to the previous method is that we need the
field computation of the unclosed model in all 3 variables (x, y, s), which implies increased computer time.
We start with the magnetic field vector B(x, y, s) representing the field of an unclosed magnet model computed using Biot-Savart law. Therefore, V- B = 0 and Vx B = 0. Then, exists another vector R(x,y, s), which stands for the field generated by fictitious line currents representing the closings of the model. Obviously, R is not unique, due to infinitely many possibilities of closing. It follows that V- R = 0 and Vx (B + R) = 0. Taking the cross product VxVx (B + R) = V(V - (B + R)) - A(B + R) = 0, we obtain A(Bi + Ri) = 0, Bi and Ri being the components in cylindrical coordinates of B and R. Now, we know from the direct method what is the structure of a function in cylindrical
coordinates that satisfies the Laplace equation. Hence, we get
^
Ri(r,4,s) = ^(fi,l (r, s) sin ¡4 + gi, l (r, s) cos ¡4) rl - Bi(r,4,s). (9)
l=0
Apparently, we get the smallest Ri in the vicinity of the optical axis if we choose the free parameters in fitl(r,s) and gill(r,s), the true multipoles, such that they cancel the corresponding terms in Bi. This way we fix uniquely the true multipoles, that are anyway the dominating part, and let Ri to contribute only for the pseudo-multipole parts. Here we define as being a true multipole of order l the s dependent function that is the coefficient of rl-1 cos ¡4 or rl-1 sin ¡4 respectively in the expression of Bi. This definition makes sense, as this is the case in general when the fields are derivable from a magnetic scalar potential.
Once the principle is understood, in practice we do not need to calculate explicitly Ri. It is enough to have Bi and extract the relevant terms, the true multipoles, then the out of axis expansion is performed, the potential is built up and the new fields are computed. The new fields will satisfy Maxwell's equation and corrections can be introduced from any part of the magnet, hence the name global Maxwellification.
Still, one thing remains to be proved. The solution is really unique if we prove that the true multipoles are invariant with respect to which component of the original field we choose, Br or B*. This result is easily obtained in case we impose the vanishing curl and divergence conditions, as has been shown in the direct method. It can be shown that this is the case without imposing any constraints on the coefficients. This is the subject of the appendix. Also, in the same appendix we stress the importance of the constraints imposed by the vanishing curl.
All the methods described have been implemented in the DA based code COSY INFINITY [15-17].
Examples and Computation of Maps. Using the methods developed in the previous sections, we computed the multipole strengths as a function of s for the LHC interaction region's High Gradient Quadrupoles. These quadrupoles have two end regions, the lead end and the return end, where the fields are s dependent. We computed the multipoles up to 28-poles for both ends. The field computation has been performed up to order 13, at 1 cm equally spaced points along the optical axis. In the following we restrict ourselves to present the results for the lead end only. Fig. 1 shows the extracted normal and skew quadrupole, duodecapole, and 20-pole components.
As previously mentioned, the map computation requires also the s derivatives of the multipoles. These are easily obtained as a by-product of the multipole extraction algorithms, because we always keep their s dependence. Derivative computation in DA in an elementary operation. It yields very accurate results without the need to resort to numerical differentiation. Some of the s derivatives of multipoles presented in Fig. 1 are
shown in Figures 2 and 3. The multipoles and their derivatives have been interpolated for plotting by a derivative preserving interpolation scheme. Also, the two multipole extraction algorithms were checked against each other and found to be in complete agreement.
Fig. 1. Multipole strengths of the lead end (shown are b2(s), b6 (s), bio(s), a2(s), ae(s) and aio(s))
210 10 MO10
-300 -200 -1Й0
-MO10 -2-1010 -3-1010
HI 200
10 1011 510»
ДО
-300 -200-100
-5-1011
-10-10 11
-1000
200"
10 10 23 5-10 23
№
-300-200-100 --510 23
0.3 0.2 0.1
-300 -200 -"100
-0.1 -0.2 -0.3
[1)6^200
Fig. 2. The s derivatives of normal multipole strengths at the lead end (shown are
b2(4) (s), b2(10) (s), b6(2) (s), b6(6) (s) and bw(2) (s))
The importance of vanishing curl has been stressed at several points throughout this paper. To show the influence on the extracted quadrupole strength and its derivatives of the effect of non-vanishing curl, we compare two cases: multipole extracted from a magnet model that generates field with non-vanishing curl, and multipole extracted from the same
aQ) 100 000 2
4-108 1 /1
. J IS 5
-300-200-100' 1 )) 200 -300
^t-108 1/
-8 10 8
№
41010 2-101° AA -1 D M LL s
-300-200^V00 1 -2 10l0\ -6-101« moo
10 1021 5-1021
Fig. 3. The s derivatives of skew multipole strengths at the lead end (shown are a2(2)(s), a2(4)(s), a2(10) (s), a6(2) (s), a6(6\s) and aio(2) (s))
Table 6. Opening aberrations (x\an), for the maps of the return and lead ends
Order, Opening aberrations (x\an)
n Return end Lead end
2 0 0
3 0.2497609620388847 0.2995410193258450
4 0.6293863605053251E - 13 — . 3054844395050998E - 13
5 0.1903231022089257 0.1962270889999381
6 0.6911276300083234E - 12 — .6125473066348172E - 12
7 0.1699451387437771 0.9019387229287767E - 01
8 —. 1364626459777126E - 09 0.2013346179046971E - 09
magnet model after all the open ended wires have been closed at "infinity" (meaning far away from the observation points). The result is contained in Fig. 4. It can be seen that although the agreement is still pretty good for the multipoles, the differences are amplified for the derivatives.
b2
10000_ -200 -100 "
8000
-50 000
6000 -75 000
4000 -100 000
2000 -125 000
s -150 000
№
-200 -100
100 200
200
4106 2-106
-200 -100 -2106
b<$
200
Fig. 4. Normal quadrupole strength, its first and second s derivatives extracted from Maxwellian and non-Maxwellian fields respectively
As an application we describe very accurate high-order map computations of s dependent fields. There are two ways to calculate maps. In the first case the following three steps are needed: using the magnet model, the field expansions at selected support points along the optical axis are computed. In case it is necessary, the Maxwellification is included in this step. Then follows the extraction of the multipoles. Finally, the multipoles are interpolated by Gaussian interpolation [18], and using the integration algorithm of COSY described in [17], the map is generated.
The alternate way's first step is the same. However, the scalar potential at support points is anyway computed in the process of Maxwellification. We use this potential to integrate the equations of the motion, with an interpolation scheme that preserves the derivatives at the support points, yielding the map.
Both methods have been implemented in COSY and they give essentially the same results. Especially at high orders the Gaussian method is faster, due to the smoothing properties of Gaussian interpolation [18]. Although to list the whole map would be too long, to get a feeling of the resulting fringe field maps for the above mentioned end regions we list the opening aberrations in both ends up to order 8 in Table 6. Once we have the maps, they can be employed for dynamics studies, which is actually the final purpose of the theory and methods developed in the present paper. We applied the methods to study the fringe field effects in the LHC. Some of the results can be found in [19].
Appendix: Structure of Br and B$ for Non-Maxwellian fields. We start with the Cartesian components
Bx(x,y,s)= E aij(s)x'V°,
ij=0 (10)
n
By(x,y,s) = E bij(s)xy0 , i,j=0
as given by the field computation, without imposing any relations among the coefficients aij and bj ; n is the order of computation. Transformation to cylindrical coordinates x = r cos y = r sin 4 gives
Br = Bx cos 4 + By sin 4, r x 4 y 4 (11) B$ = -Bx sin 4 + By cos 4-
Inserting (10) in (11) we obtain
n
Br (r, 4, s) = ^^ ri+j cos® 4 sinj 4 (aij (s) cos 4 + hj (s) sin 4),
i,j=0
=. (12) B$(r, 4, s) = ^^ ri+j cos® 4sinj 4 (-aij(s) sin 4 + b®j(s) cos 4) -
i,j=0
The next step is to transform the products of trigonometric functions into a sum of trigonometric functions involving multiple angles via
i+j + 1
cosi+1 4 sinj 4 = aj cos k4 + f3kj) sin k4,
k=0
i+j+1
cosi $ sinj+1 $ = Yk^ cos k$ + fi^ sin k$,
j oi,
k=0
(j) (j) (j) (j) where a£:), j3kj, y(), are real constants depending on j. Hence,
Br = £ r
i,j=0
i+j
i+j+1
j(s) (aj cos k$ + 3j) sin k$) +
k=0 i+j + 1
+ bij (s) (Ykj) cos k$ + 4j) sin k$
B, = £ r
i,j=0
i+j
k=0 i+j + 1
(13)
-aij (s) £ (Ykj cos k$ + 5j sin k$) +
k=0 i+j + 1
+ bij (s) Y^ (aj cos k$ + 3j) sin k$)
k=0
Now we can use the multipole definition (9) to retain in Br and B, only the true multipoles, i. e. we need to keep only the terms with k = i + j + 1, all others giving rise to pseudo-multipoles. We neglect the solenoidal terms, which are always treated best separately from Bs. Then, the components of the field containing just the true multipoles are
Br
B, =
ri+j
i,j=0
r
i,j=0
i+j
aij (sK(+j+1 + bij (s)Yi+j+i) cos (i + j + 1)$ +
+ (aij (s)3ffi+1 + bij (s^^ sin (i + j + 1) $] , (-aij (s)y j)j+1 + bij (s)a(+)j+^ cos (i + j + 1) $ +
+ (-aij (s)fi(+)j+1 + bij (s)j^+0 sin (i + j + 1) $ .
By expanding the trigonometric functions in terms of exponentials, it can be seen that: for j even
a(j) =2-(i+j) (_l)j/2 3(j) =0 Y(j) =0 S(j) =2-(i+j) (_l)j/2 ai+j+1 = 2 ( 1) , 3i+j + 1 = 0, Yi+j+1 = 0, °i+j+1 = 2 ( 1) ,
and for j odd
*i+J+1
= 0, 3^+1 = 2-(i+j) (-1)(j-1)/2 , Yjj+1 = -2-(i+j) (-1)(j-1)/2 , j =0.
Separation of the double sum into summation over i- and j-even, respectively j-odd leads to
Br Yl [aij(s)cos(i + j + 1)$ + bij(s)sin(i + j + 1)$]2-(i+j) (-1)j/2 +
i=0 I j=0 j-even
n
n
n
n
n
n
+ E ri+j -bij(s)cos(i + j + 1)4 + an(s)sin(i + j + 1)4] (-1)(j-1)/2
j=1 j-odd
n I n
= X) S E ri+j [bjj(s) c°s(* + J + ^ - ajj(s) sin( + ' + ^ 2~(i+j) i"1)''7' +
i=0 I j=0 j-even n
+ Y; ri+j K (s) cos(» + J + 1)^ + bij (s) sin(» + ' + 1)4] 2-(i+j) (-1)(j-1)72
j=1 j-odd
The symmetry of the above equations is clear. Br and B$ have the same number of terms, and if the symmetry holds term by term, then it also holds for their sum. By introducing a new index l = » + J we obtain relations of the form
Br = J2A (s)cos(l +1)4 + Bi (s)sin(l + 1)4] rl
l=0 n
Bf = Y,[Bi (s)cos(l + 1)4 - Ai (s)sin(l +1)4] rl
l=0
where Al (s), Bl (s) are sums over aij (s), bij (s) with i + j = l. Shifting the origin of summation to make comparison easier with the direct method, and using the convention that the l = 1 component corresponds to dipole gives the final form
n+1
Br(r,4,s) = J2[Al (s) cos 14 + Bl (s) sin 14] rl-1, (14)
l=i n+1
Bf(r,4,s) = [Bl(s)cos 14 - Al (s) sin 14] rl
l-1
l=1
By identification, it is apparent that (14) is of the same form as the field components derived from a scalar potential, containing only the true multipoles. For example, up to order 5 we can derive the following relations:
A1(s) A2(s) = aoo(s), = ^(aio(s)- bo1(s)),
A3(s) = ^ (a2°(s) ~ bn(s) - ao2(s)) ,
A4(s) = £ (a3o(s) - b21(s) - a 12(s) + bo3(s)),
As(s) = ^(a4o(s) - b31(s) - a22(s) + b13(s) + ao4(s)) ,
B1(s) B2(s) = boo(s), = ^ (6io(s) + ao1(s)) ,
Bs(s) = \(b2o(s) + an(s) - bo2(s)) ,
n
n
Ba{S) = i (630(5) + fl2i(s) - 612(5) - 003(5)),
B5(s) = -j- (b4o(s) + a31(s) - 622(5) - ai3(s) + &04(«)) • 16
The differences between the constrained case (by the Maxwell equations) and arbitrary coefficients now can be analyzed. Clearly, the dipole component will give the same result in every method. The differences start to show up beginning with the quadrupole component. For example the normal quadrupole is given in general by B2(s) = (b10(s) + a01(s))/2.
If we impose V-B = 0, as it is always the case for magnetic field computations, it gives just a10(s) + b01(s) + c00(s) = 0, that is, it does not impose any constraints between b10(s) and a01(s). On the other hand, if Vx B = 0, the s component imposes: b10(s) = a01(s). If the curl is not vanishing, i. e. b10(s) = a01(s), the method will take as the quadrupole component the average value. The same type of analysis can be performed on higher order multipoles to emphasize the importance of vanishing curl.
As a conclusion, we proved that the method can be used for global Maxwellification, with a unique solution, that alters the original fields by a minimal amount. However, we remind the reader that Br and B, do not contain all the terms, the whole field expressions for Br and B, have contributions from pseudo-multipoles that cannot be written in the form of (6) in case of non-vanishing curl.
Acknowledgments. We would like to thank G. Sabbi for many useful discussions and the considerable amount of data supplied for the Large Hadron Collider's interactions region High Gradient Quadrupoles. We would also like to thank Kyoko Makino for numerous discussions about implementational aspects in DA codes, Carol Johnstone for many discussions about FFAGs, Shashikant Manikonda for discussions of novel 3D coil arrangements to achieve very general field configurations, and various members of the Jiilich JEDI collaboration for many interesting discussions about dedicated elements for the manipulation of spin dynamics. The work was supported by the U.S. Department of Energy.
References
1. The LHC Study Group. The Large Hadron Collider — Conceptual Design. Technical Report AC/95-05 (LHC), CERN, 1995, 215 p.
2. Wei J., Ptitsin V., Pilat F., Tepikian S., Gelfand N., Wan W., Holt J. US-LHC IR magnet error analysis and compensation. Proceedings EPAC-1998. Stockholm, 1998, pp. 380—382.
3. Russenschuck S. A computer program for the design of superconducting accelerator magnets. Technical Report CERN AT/95-39, CERN, 1995, 366 p.
4. Johnstone C., Berz M., Snopok P. (eds.). Proceedings FFAG-2009 Intern. Conference. Singapore, World Scientific, 2011, 249 p. DOI:10.1142/S0217751X11053079.
5. Storage Ring Electric Dipole Moment Collaboration. URL: https://www.bnl.gov/edm/ (accessed 27.10.2015) and also JEDI Collaboration (Jiilich Electric Dipole Moment Investigations). URL: http://collaborations.fz-juelich.de/ikp/jedi/ (accessed 27.10.2015).
6. Sabbi G. Magnetic field analysis of HGQ coil ends. Technical Report TD-97-040. Batavia, Fermilab, 1997, 17 p.
7. Caspi S., Helm M., Laslett L. J., Brady V. O. An approach to 3D magnetic field calculation using numerical and differential algebra methods. Technical Report LBL-32624. Berkeley, Lawrence Berkeley Laboratory, 1992, 68 p.
8. Berz M. Modern Map Methods in Particle Beam Physics. San Diego, Academic Press, 1999. 318 p. Also available at http://bt.pa.msu.edu/pub (accessed 27.10.2015).
9. Berz M. Arbitrary order description of arbitrary particle optical systems. Nuclear Instruments and Methods, 1990, vol. A298, pp. 426-440.
10. Berz M. Differential algebraic description of beam dynamics to very high orders. Particle Accelerators, 1989, vol. 24, pp. 109-124.
11. Lindemann M. Rigorous Field Analysis of Superconducting Magnets and the Influence on
Nonlinear Dynamics in Particle Accelerators. Master's thesis. East Lansing, Michigan, USA, Michigan State University, 1998, 75 p.
12. Erdeiyi B. Symplectic Approximation of Hamiltonian Flows and Accurate ¡Simulation of Fringe Field Effects. PhD thesis. East Lansing, Michigan, USA, Michigan State University, 2001, 262 p.
13. Manikonda S. High Order Finite Element Methods to Compute Taylor Transfer Maps. PhD thesis. East Lansing, Michigan, USA, Michigan State University, 2006, 146 p.
14. Berz M. Differential Algebraic Techniques. Handbook of Accelerator Physics and Engineering. Eds. A. Chao, M. Tigner. Singapore, World Scientific, 1999, ch. 2.3.7, pp. 84-87.
15. Berz M. COSY INFINITY, a new arbitrary order optics Code. Computer Codes and the Linear Accelerator Community, Los Alamos LA-11857-C, 1990, pp. 137-156.
16. Berz M. Computational aspects of optics design and simulation: COSY INFINITY. Nuclear Instruments and Methods, 1990, vol. A298, pp. 473-479.
17. Makino K., Berz M. COSY INFINITY version 8. Nuclear Instruments and Methods, 1999, vol. A427, pp. 338-343.
18. Makino K., Berz M. Arbitrary order aberrations for elements characterized by measured fields. Optical Science, Engineering and Instrumentation '97, 1997, vol. 3155, pp. 221-227.
19. Erdelyi B., Berz M., Makino K. Detailed analysis of fringe field effects in the large hadron collider. Technical Report MSUCL-1129. National Superconducting Cyclotron Laboratory. East Lansing, Michigan, USA, Michigan State University, 1999, 37 p.
Статья рекомендована к печати проф. Д. А. Овсянниковым. Статья поступила в редакцию 10 сентября 2015 г.
