- Email: [email protected]
An interpolation scheme for cylindrically symmetric magnetic fields
NUCLEAR
INSTRUMENTS
AND METHODS
148 ( 1 9 7 8 )
163-166
; (~) N O R T H - H O L L A N D
PUBLISHING
CO.
AN INTERPOLATION SCHEME FOR CYLINDRICALLY SYMMETRIC MAGNETIC FIELDS D. E. LOBB
TRIUMF, Physics Department, University Of Victoria, Victoria, B.C., Canada Received 6 June 1977 A simple method based on Maxwell's equations and cylindrical symmetry is used to interpolate between values of a cylindrically symmetric magnetic field tabulated over a rectangular grid.
constants are calculated. For reasons to be explained below, the interpolation constants for Br and B, are calculated separately. c) Calculates the desired B, and B~ values using the appropriate interpolation constants. d) Converts the Br and B, values into the required Cartesian components Bx, By, B=. e) Returns these values to the Runge-Kutta routine.
I. Introduction
The magnet design program GFUN [Armstrong et al.L2)] provides for the tracking of particle trajectories through the magnetic field produced by a configuration of iron and coils. This tracking is quite slow since the fourth-order Runge-Kutta-Merson procedure requires the calculation of field values at five locations for each segment of Lhe ray. In this paper we present a faster method of doing these ray-tracing calculations for the particular case of axi-symmetric systems. Firstly, field values are calculated by GFUN over a rectangular mesh in the x_> 0, y = 0 half plane. (The z axis is Lhe axis of symmetry.) This yields values of Bx and .B=, with B,.= 0 by cylindrical symmetry. For the :subsequent calculations using cylindrical symmetry, the calculated Bx is interpreted as Br, with B 0 = 0 by cylindrical symmetry. These B r and B, field values are stored and are available to the interpolation scheme. The GFUN trajectory tracking routine uses ,Cartesian coordinates. Ray tracing in Cartesian coordinates avoids a potential source of numerical difficulty that can arise when cylindrical coordinates are used: there is a term in r-3 in the differential equation for the r-component of the trajectory which can get very large for trajectories which pass close to the symmetry axis. For further information on this point the interested reader is referred to the papers by Lindgren et al. 3) and Lindgren and SchneideP). In our calculation procedure the Runge-Kutta routine calls the interpolation routine which performs the following operations: a) Changes the Cartesian coordinates (x, y, z) to cylindrical coordinates (r, ~, z). b) Checks if either or both of the previous sets of interpolation constants are appropriate at this new point. If not, new interpolation
2. The interpolation method The combination of Maxwell's equations and cylindrical symmetry yields, for a static magnetic field in a source-free region of uniform magnetic permeability, the following equation for the vector potential [Olsen et al.S)]. A4' = ~=o ( - 1 ) V v ! ( v + l ) !
--
B(°2~)(z)'
(1)
which gives the magnetic field components 8r
=
--
v=o (--1)~ (l+v)-------~.' v! --
B,~ = 0, B~ = ~=o ( - 1)~ ~
Bg2v+l)(Z),
(2)
(3) B(°2~)(z)'
(4)
where Bo(z) is the axial magnetic field on the z-axis and the bracketed superscripts on B0 (z) represent derivatives. These equations have been used frequently and Olsen et al. 5) have commented on the inaccuracy that results in practical cases when derivatives are taken of a global polynomial fit to axial field data. Since higher order derivatives obtained from numerical data have large error, the coupling of these terms with large powers of r in eqs. (1)--(4) can introduce errors into calculated non-paraxial
164
D.E.
magnetic field values. Recently, Sk611ermo 6) has shown that the field components may be expressed in terms of complex integrals which are equivalent to the infinite series of eqs. (2) and (4) but which retain accuracy in regions not adjacent to the axis. However, in the method presented here, a global fit is not attempted but the form of the equations is used to obtain a local relationship in the vicinity of the point of interest. We are interested in & and B. at point P. The tabular entry which lies closest to P is that at 0 (r0, z0). We change variables: p = r/r o, (5) (6)
= (Z-Zo)/ro.
So point 0 now has coordinates p = 1, ff = 0. We assume that the function B0 (z), which is appropriate to this particular point P may be written as a fourth-order polynomial in ~: Bo(O = ~ + tiff + ~ 2 + 6~3 + e ~ .
(7)
[Note that B0 (~) need not have the interpretation of the axial magnetic field on the z-axis; here, it is just a general interpolation function.] This polynomial expansion yields B~4~(~)= constant, with higher derivatives zero. Thus, eq. (4) has terms in r
0 xp
z
r
i- . . . . . . . . . . . . . . . . . . . . . . .
xP 0
Z
Fig. 1. For point P lying within the dashed rectangle, the tabulated B z values used to determine B z at point P are located at the circled mesh intersections.
LOBB
p0, p2 and p4; eq. (2) has terms in p and results are: p
Br = - 7 ~ - ( p O
The
~ + 6(-3P~ 2 +~P3)+
+ e(-Zp~3 + 3p3~), B: = ~ + / ~
p3.
(8)
+ ~,(~2 _ ½p2) + 6~(~2 _ ~p2) +
+ e(~4 _ 3p2~2 + 3 p 4 ) .
(9)
In m a n y practical cases the variation of B: with is fairly slow, so eq. (9) is dominated by ~z. The fitted values of/3, y and 6 to be obtained from eq. (9) are likely to be inaccurate for the purposes of eq. (8). Thus, the B,. and Bz interpolations are carried out completely independently using tabulated B,. and B: values respectively. For the B,. interpolations we need appropriate values of/3, y, 6 and e. These are obtained using the tabulated B,. values at the four mesh points that form the corners of the mesh rectangle in which point P lies. If two of these mesh points lie on the symmetry axis, the values at the other two corners of the mesh rectangle are used to determine/3 and y with 6 and e set equal to zero. This lower order of accuracy is acceptable since Br = 0 on the symmetry axis and ]Brl is quite small near the symmetry axis. The tabular entries used to determine ~z,/3, y, 6 and e for the B~ calculation are illustrated in fig. 1. The two sets of linear equations obtained may be solved by any standard procedure to yield the separate interpolation constants for Br and Bz appropriate to point P. 3. A n e x a m p l e
In fig. 2 we present a four coaxial coil geometry and in figs 3 and 4 we present some representative field values calculated using GFUN. In fig. 5 we present the trajectories used to test the interpolation procedure. The trajectories were calculated in all cases using the same R u n g e - K u t t a method. The quantity of interest is the location where the trajectory intersects the z-axis. This is a severe test of the interpolation scheme since trajectory errors will accumulate. The results are presented in table 1. We see that a radial and axial spacing of 10 cm between tabular entries provides an accuracy of 0.04 m m in locating the cross-over location. To cover the region of space traversed by these trajec-
INTERPOLATION
165
SCHEME
r60"]
X
X --+
+ -400
-500
-i-300
+ 200
X I
H
O0
60 + 0
+ 10C
+
+ 300
+ 400
+ ~ 500 Z(cm)
X
I
Z =-450
+ 200
I
I
Z = -150
Z = 150
Z=450
Fig. 2. A system of four coaxial coils, current density 721 A / c m 2 in each coil.
tories requires 7 x 3 2 = 224 tabular entries. The Runge-Kutta procedure requires field values at five locations for each calculation interval: for a calculation interval of 10cm, field values at 150 24
I
2:2
I
r~5Ocm
20
f•
18 16
G
14 12
N I0 8
2 0
I
-150
I
-tO0
I
-50
0 Z icm)
I
I
50
100
150
locations are required for each trajectory spanning - 150 cm~~0.1 min on an IBM 370/158 for typical problems). In this situation a 30 segment trajectory would require > 12rain of computation time. With the interpolation method, once the table of field values has been written, all trajectories require the same calculation time regardless of the complexity of the magnetic geometry. The interpolation method allows a complete beam-optical examination of axi-symmetric systems, no matter how complicated.
Fig. 3. The longitudinal magnetic field at three different radii produced by the coil system of fig. 2. 60
i
0
I
I
-80
-40
I
I
I
40
80
,0
r~5Ocrn
G
I
i
~0
f
E :O
-2 -3
0
-4 -5 -6 -150
O -160 -100
-
0
0 Z (cm)
50
100
150
Fig. 4. The radial magnetic field at three different radii produced by the coil system of fig. 2.
-120
0 Z (cm)
Fig. 5. T h e radial position as a f u n c t i o n
120
160
o f distance for t h r e e
trajectories in the magnetic field produced by the coil system of f i g . 2 .
166
D.E.
LOBB
TABLE 1 The z value (cm/where the particular trajectory crosses the symmetry axis. The interval of calculation used for the Runge-Kutta procedure (3w) was a constant 2, 5 or 10 cm. Near the end of the trajectory a fine spacing of Aw = 0.1 cm was used; the location of the cross-over location was obtained by linear interpolation of these finely spaced results. The GFUN result was obtained by the Runge-Kutta procedure calling the GFUN field calculation procedure for each required value of the magnetic field. Trajectory A (P0 = 145 MeV/c, 00 = 35°) Radial spacing of tabular entries (cm)
Aw (cm)
Axial spacing of tabular entries (cm) 5 10 15
(Po
Trajectory B 140 MeV/c, 0o = 45°)
Axial spacing of tabular entries (cm) 5 10 15
Trajectory C (P0=146MeV/c, 0o=55 °) Axial spacing of tabular entries (cm) 5 10 15
5
2 5 10
149.83 149.84 149.84
149.82 149.82 149.82
149.79 149.79 149.75
150.02 150.02 150.02
150.03 150.03 150.02
150.03 150.03 150.02
149.97 149.97 149.97
149.96 149.96 149.96
149.93 149.96 149.94
10
2 5 10
149.75 t49.78 149.76
149.81 149.81 149.81
149.82 149.81 149.81
150.02 150.02 150.03
150.02 150.02 150.03
150.03 150.03 150.02
150.02 150.00 150.02
149.97 149.97 149.96
149.91 149.93 149.90
15
2 5 10
149.97 149.97 150.05
149.85 149.85 149.84
149.82 149.82 149.82
150.07 149.95 149.97
150.03 150.02 150.05
150.04 150.04 150.03
149.94 149.94 149.85
149.85 149.85 149.86
149.84 149.84 149.84
GFUN results
2 5 10
149.85 149.85 149.84
The author wishes to thank B. Henin for programming and testing the interpolation procedure, and C. W. Trowbridge for providing us with a copy of the program GFUN. References I) A. G. A. M. Armstrong, C. J. Collie, N. J. Diserens, M. J. Newman, J. Simkin and C. W. Trowbridge, Proc. Int. Conf. on Magnet technology (Rome, 1975) p. 168,
150.03 150.02 150.02
149.97 149.97 149.96
2) A. G. Armstrong, C. J. Collie, N. J. Diserens, M. J. Newman, J. Simkin and C. W. Trowbridge, Rutherford Laboratory Report RL-76-029/A (1976). 3) I. Lindgren, G. Petterson and W. Schneider, Nucl. Instr. and Meth. 22 (1963) 61. 4) I. Lindgren and W. Schneider, Nucl. Instr. and Meth. 22 (1963) 48. s) B. Olsen, G. Petterson and W. Schneider, Nucl. Instr. and Meth. 41 (1966) 325. 6) A. SkNlermo, Nucl. Instr. and Metb. 137 (1976) 339.
INSTRUMENTS
AND METHODS
148 ( 1 9 7 8 )
163-166
; (~) N O R T H - H O L L A N D
PUBLISHING
CO.
AN INTERPOLATION SCHEME FOR CYLINDRICALLY SYMMETRIC MAGNETIC FIELDS D. E. LOBB
TRIUMF, Physics Department, University Of Victoria, Victoria, B.C., Canada Received 6 June 1977 A simple method based on Maxwell's equations and cylindrical symmetry is used to interpolate between values of a cylindrically symmetric magnetic field tabulated over a rectangular grid.
constants are calculated. For reasons to be explained below, the interpolation constants for Br and B, are calculated separately. c) Calculates the desired B, and B~ values using the appropriate interpolation constants. d) Converts the Br and B, values into the required Cartesian components Bx, By, B=. e) Returns these values to the Runge-Kutta routine.
I. Introduction
The magnet design program GFUN [Armstrong et al.L2)] provides for the tracking of particle trajectories through the magnetic field produced by a configuration of iron and coils. This tracking is quite slow since the fourth-order Runge-Kutta-Merson procedure requires the calculation of field values at five locations for each segment of Lhe ray. In this paper we present a faster method of doing these ray-tracing calculations for the particular case of axi-symmetric systems. Firstly, field values are calculated by GFUN over a rectangular mesh in the x_> 0, y = 0 half plane. (The z axis is Lhe axis of symmetry.) This yields values of Bx and .B=, with B,.= 0 by cylindrical symmetry. For the :subsequent calculations using cylindrical symmetry, the calculated Bx is interpreted as Br, with B 0 = 0 by cylindrical symmetry. These B r and B, field values are stored and are available to the interpolation scheme. The GFUN trajectory tracking routine uses ,Cartesian coordinates. Ray tracing in Cartesian coordinates avoids a potential source of numerical difficulty that can arise when cylindrical coordinates are used: there is a term in r-3 in the differential equation for the r-component of the trajectory which can get very large for trajectories which pass close to the symmetry axis. For further information on this point the interested reader is referred to the papers by Lindgren et al. 3) and Lindgren and SchneideP). In our calculation procedure the Runge-Kutta routine calls the interpolation routine which performs the following operations: a) Changes the Cartesian coordinates (x, y, z) to cylindrical coordinates (r, ~, z). b) Checks if either or both of the previous sets of interpolation constants are appropriate at this new point. If not, new interpolation
2. The interpolation method The combination of Maxwell's equations and cylindrical symmetry yields, for a static magnetic field in a source-free region of uniform magnetic permeability, the following equation for the vector potential [Olsen et al.S)]. A4' = ~=o ( - 1 ) V v ! ( v + l ) !
--
B(°2~)(z)'
(1)
which gives the magnetic field components 8r
=
--
v=o (--1)~ (l+v)-------~.' v! --
B,~ = 0, B~ = ~=o ( - 1)~ ~
Bg2v+l)(Z),
(2)
(3) B(°2~)(z)'
(4)
where Bo(z) is the axial magnetic field on the z-axis and the bracketed superscripts on B0 (z) represent derivatives. These equations have been used frequently and Olsen et al. 5) have commented on the inaccuracy that results in practical cases when derivatives are taken of a global polynomial fit to axial field data. Since higher order derivatives obtained from numerical data have large error, the coupling of these terms with large powers of r in eqs. (1)--(4) can introduce errors into calculated non-paraxial
164
D.E.
magnetic field values. Recently, Sk611ermo 6) has shown that the field components may be expressed in terms of complex integrals which are equivalent to the infinite series of eqs. (2) and (4) but which retain accuracy in regions not adjacent to the axis. However, in the method presented here, a global fit is not attempted but the form of the equations is used to obtain a local relationship in the vicinity of the point of interest. We are interested in & and B. at point P. The tabular entry which lies closest to P is that at 0 (r0, z0). We change variables: p = r/r o, (5) (6)
= (Z-Zo)/ro.
So point 0 now has coordinates p = 1, ff = 0. We assume that the function B0 (z), which is appropriate to this particular point P may be written as a fourth-order polynomial in ~: Bo(O = ~ + tiff + ~ 2 + 6~3 + e ~ .
(7)
[Note that B0 (~) need not have the interpretation of the axial magnetic field on the z-axis; here, it is just a general interpolation function.] This polynomial expansion yields B~4~(~)= constant, with higher derivatives zero. Thus, eq. (4) has terms in r
0 xp
z
r
i- . . . . . . . . . . . . . . . . . . . . . . .
xP 0
Z
Fig. 1. For point P lying within the dashed rectangle, the tabulated B z values used to determine B z at point P are located at the circled mesh intersections.
LOBB
p0, p2 and p4; eq. (2) has terms in p and results are: p
Br = - 7 ~ - ( p O
The
~ + 6(-3P~ 2 +~P3)+
+ e(-Zp~3 + 3p3~), B: = ~ + / ~
p3.
(8)
+ ~,(~2 _ ½p2) + 6~(~2 _ ~p2) +
+ e(~4 _ 3p2~2 + 3 p 4 ) .
(9)
In m a n y practical cases the variation of B: with is fairly slow, so eq. (9) is dominated by ~z. The fitted values of/3, y and 6 to be obtained from eq. (9) are likely to be inaccurate for the purposes of eq. (8). Thus, the B,. and Bz interpolations are carried out completely independently using tabulated B,. and B: values respectively. For the B,. interpolations we need appropriate values of/3, y, 6 and e. These are obtained using the tabulated B,. values at the four mesh points that form the corners of the mesh rectangle in which point P lies. If two of these mesh points lie on the symmetry axis, the values at the other two corners of the mesh rectangle are used to determine/3 and y with 6 and e set equal to zero. This lower order of accuracy is acceptable since Br = 0 on the symmetry axis and ]Brl is quite small near the symmetry axis. The tabular entries used to determine ~z,/3, y, 6 and e for the B~ calculation are illustrated in fig. 1. The two sets of linear equations obtained may be solved by any standard procedure to yield the separate interpolation constants for Br and Bz appropriate to point P. 3. A n e x a m p l e
In fig. 2 we present a four coaxial coil geometry and in figs 3 and 4 we present some representative field values calculated using GFUN. In fig. 5 we present the trajectories used to test the interpolation procedure. The trajectories were calculated in all cases using the same R u n g e - K u t t a method. The quantity of interest is the location where the trajectory intersects the z-axis. This is a severe test of the interpolation scheme since trajectory errors will accumulate. The results are presented in table 1. We see that a radial and axial spacing of 10 cm between tabular entries provides an accuracy of 0.04 m m in locating the cross-over location. To cover the region of space traversed by these trajec-
INTERPOLATION
165
SCHEME
r60"]
X
X --+
+ -400
-500
-i-300
+ 200
X I
H
O0
60 + 0
+ 10C
+
+ 300
+ 400
+ ~ 500 Z(cm)
X
I
Z =-450
+ 200
I
I
Z = -150
Z = 150
Z=450
Fig. 2. A system of four coaxial coils, current density 721 A / c m 2 in each coil.
tories requires 7 x 3 2 = 224 tabular entries. The Runge-Kutta procedure requires field values at five locations for each calculation interval: for a calculation interval of 10cm, field values at 150 24
I
2:2
I
r~5Ocm
20
f•
18 16
G
14 12
N I0 8
2 0
I
-150
I
-tO0
I
-50
0 Z icm)
I
I
50
100
150
locations are required for each trajectory spanning - 150 cm~
Fig. 3. The longitudinal magnetic field at three different radii produced by the coil system of fig. 2. 60
i
0
I
I
-80
-40
I
I
I
40
80
,0
r~5Ocrn
G
I
i
~0
f
E :O
-2 -3
0
-4 -5 -6 -150
O -160 -100
-
0
0 Z (cm)
50
100
150
Fig. 4. The radial magnetic field at three different radii produced by the coil system of fig. 2.
-120
0 Z (cm)
Fig. 5. T h e radial position as a f u n c t i o n
120
160
o f distance for t h r e e
trajectories in the magnetic field produced by the coil system of f i g . 2 .
166
D.E.
LOBB
TABLE 1 The z value (cm/where the particular trajectory crosses the symmetry axis. The interval of calculation used for the Runge-Kutta procedure (3w) was a constant 2, 5 or 10 cm. Near the end of the trajectory a fine spacing of Aw = 0.1 cm was used; the location of the cross-over location was obtained by linear interpolation of these finely spaced results. The GFUN result was obtained by the Runge-Kutta procedure calling the GFUN field calculation procedure for each required value of the magnetic field. Trajectory A (P0 = 145 MeV/c, 00 = 35°) Radial spacing of tabular entries (cm)
Aw (cm)
Axial spacing of tabular entries (cm) 5 10 15
(Po
Trajectory B 140 MeV/c, 0o = 45°)
Axial spacing of tabular entries (cm) 5 10 15
Trajectory C (P0=146MeV/c, 0o=55 °) Axial spacing of tabular entries (cm) 5 10 15
5
2 5 10
149.83 149.84 149.84
149.82 149.82 149.82
149.79 149.79 149.75
150.02 150.02 150.02
150.03 150.03 150.02
150.03 150.03 150.02
149.97 149.97 149.97
149.96 149.96 149.96
149.93 149.96 149.94
10
2 5 10
149.75 t49.78 149.76
149.81 149.81 149.81
149.82 149.81 149.81
150.02 150.02 150.03
150.02 150.02 150.03
150.03 150.03 150.02
150.02 150.00 150.02
149.97 149.97 149.96
149.91 149.93 149.90
15
2 5 10
149.97 149.97 150.05
149.85 149.85 149.84
149.82 149.82 149.82
150.07 149.95 149.97
150.03 150.02 150.05
150.04 150.04 150.03
149.94 149.94 149.85
149.85 149.85 149.86
149.84 149.84 149.84
GFUN results
2 5 10
149.85 149.85 149.84
The author wishes to thank B. Henin for programming and testing the interpolation procedure, and C. W. Trowbridge for providing us with a copy of the program GFUN. References I) A. G. A. M. Armstrong, C. J. Collie, N. J. Diserens, M. J. Newman, J. Simkin and C. W. Trowbridge, Proc. Int. Conf. on Magnet technology (Rome, 1975) p. 168,
150.03 150.02 150.02
149.97 149.97 149.96
2) A. G. Armstrong, C. J. Collie, N. J. Diserens, M. J. Newman, J. Simkin and C. W. Trowbridge, Rutherford Laboratory Report RL-76-029/A (1976). 3) I. Lindgren, G. Petterson and W. Schneider, Nucl. Instr. and Meth. 22 (1963) 61. 4) I. Lindgren and W. Schneider, Nucl. Instr. and Meth. 22 (1963) 48. s) B. Olsen, G. Petterson and W. Schneider, Nucl. Instr. and Meth. 41 (1966) 325. 6) A. SkNlermo, Nucl. Instr. and Metb. 137 (1976) 339.