SAD RMS Envelope Simulation Manual
Christopher K. Allen and Masanori Ikegami
KEK, Tsukuba, Ibaraki 305-0801 Japan
December, 2005

The capability for RMS envelope simulation, including space charge, has been
implemented in the SAD environment. This note is a summary of this feature. Included
is a discussion on running the simulation, along with a description of the simulation
parameters. Also included are a listing of the most used functions and their arguments, alisting of all the new functions implemented, along with a brief description and theircalling hierarchy, and a code excerpt demonstrating how to use the simulation in
SADScript.

1

Table of Contents
SAD RMS ENVELOPE SIMULATION PRIMER ................................................ 1
1

NOTES ON RUNNING THE SIMULATION ................................................... 3
1.1
1.2
1.3
1.4

Running the Simulation.................................................................................... 3
Arguments of ScheffSimulate .......................................................................... 3
Simulation Data................................................................................................. 5
Modifying Beamline Parameters ..................................................................... 6

2

BEAM PERVEANCE AND MULTI-CHARGED PARTICLE SIMULATION ..... 6

3

THE SECOND-ORDER MOMENT MATRIX ............................................... 7

4 IMPORTANT SAD FUNCTIONS AND THEIR LOCATIONS ............... ERROR!
BOOKMARK NOT DEFINED.
5FUNCTION LIST AND CALLING HIERARCHY ........................................... 9
5.1
5.2
5.3
5.4
5.5

Module Scheff.n................................................................................................. 9
Module Trace3dToSad.n................................................................................ 10
Module TwissUtility.n..................................................................................... 10
Module MatrixFunctions.n ............................................................................ 10Internal Additions to the SAD Interpreter................................................... 11

6

EXAMPLE CODE ..................................................................................... 127

OUTPUT .................................................................................................. 152

1 Notes on Running the Simulation
1.1

Units
The simulation uses all MKS units, except for beam energies and particle charges.
Particle charges are normalized by the unit charge q so that they appear as integral
quantities, positive or negative. Beam energies are in electron-Volts (including massenergies). All statistical quantities are given in the RMS values, for example, allemittances are specified as RMS emittances rather than effective emittance or normalized
emittance.
1.2

Running the Simulation
The general procedure for running the simulation is to load the input deck, specify the
initial beam parameters, convert them from Trace3D units to SAD units if necessary, thenrun the simulation with a call to ScheffSimulate[]. The results of the simulation
are returned by the function and may be processed as desired.
Before the call to ScheffSimulate[] can be made, there are several SAD
environment parameters that must be set. We list them below with at short explanation.
TRPT;
INS;
NOCOD;RFSW;! Simulate a transport section! Treat lattice as an insertion device (not part of a ring)! No closed orbit dispersion! Radio frequency standing wave – otherwise traveling wave structures

Also, to use the function PlotBeamBeta[] you should set the following parameter$DisplayFunction = CanvasDrawer;
1.3

Arguments of ScheffSimulate[]
ScheffSimulate[] uses an adaptive stepping procedure. Consequently, there
are several numerical parameters, as well as the beam parameters, that can be passed to
ScheffSimulate[] to fine tune this integration process. All of the numericalparameters are optional parameters having default values. These default values havebeen set for a reasonable performance/precision tradeoff. However, for some cases you
may wish to experiment to see which gives you the better performance and/or solution.
Consequently, we list here all the arguments, their description, and their default values ifthey are optional.
{{sn},{γn},{σn}} = ScheffSimulate[K0, σ0, h0:0.01, εsoln:10-5, ∆h:0.05, hmax:0.0]
o K0 = Initial generalized beam perveance of the beam, that is, the perveance at theentrance of the beamline. This parameter determines the space charge force and is3

dependent upon particle energy. For a detailed description of generalized perveanceK see Section 2.
o σ0 = Initial second-order moment matrix of the beam. This is the 6×6 real, symmetricmatrix containing all the second-order moments of the particle beam (with respect tothe beam distribution). For a description of the initial moment matrix σ0, and themoment matrix σ in general see Section 3.
o h0 (optional) = The initial step size to used to start the adaptive integration algorithm.
The integration algorithm continually adjusts the actually step size h to maintain theerror tolerance specified by the argument εsoln (see below). Choosing an initial step
close too small will marginally slow the solution time, while choosing a value toolarge will not significantly affect the solution time. The current default value of h0 is
1 cm.
o εsoln (optional) = The integration algorithm is designed to maintain a specific accuracy
in the computed solution. Morever, it is designed to keep the step size h as large as
possible such that ||σsim − σsoln|| ≤ εsoln at each step, where σsim is the simulated
moment matrix and σsoln is the “exact” moment matrix. This technique provides thefastest solution time while maintaining the given error tolerance. The current default
value of εsoln is 10−5.
o ∆h (optional) = Slack tolerance parameter. For each integration step in the simulation,a new value of the step size h’ is computed. If h’ is greater than the previous step h,
the integration step must be rolled back and the solution recomputed with new step
size h’ in order to maintain the solution tolerance εsoln. However, if h’ is only
marginally smaller than h (perhaps due to noise or rounding errors), this is asignificant waste of CPU time. Thus, we only change the integration step size if|h’−h|/h > ∆h. (The reason for the absolute value is that we also must recompute theelement sub-transfer matrix exp(hA) if h is changed in either direction.) Thus, ∆hprovides slack for h to “jitter” before the adaptive stepping is triggered.
Consequently, ∆h can save significant CPU time, however, choosing a value to large
will compromise the accuracy of the solution and the solution tolerance εsoln. The
current default value of ∆h is 5%.
o hmax (optional) = Maximum allowable step size. For whatever reason, it is possible
for the user to specify a maximum step size for the adaptive stepping algorithm. Thisvalue will prevent the algorithm from taking any step sizes h greater than hmax. To
turn off this feature set hmax = 0. The current default value of hmax is 0 (no maximum
step size).
1.4

Results of ScheffSimulate[]
Here we list the returned values of ScheffSimulate[] along with a description.
No output of ScheffSimulate[] is stored in the SAD environment. It is all returnedby the function as a list of objects.

4{{sn},{γn},{σn}} = ScheffSimulate[K0, σ0, h0:0.01, εsoln:10-5, ∆h:0.05, hmax:0.0]o {sn} = Set of element entrance and exit positions along the beamline. The firstvalue in this list, s0, is the entrance position of the first beamline element, thatis, element n = 0. Typically this value is s0 = 0. The following values are theentrance positions of the beamline elements proceeding downstream. Notethat the exit position of element n is the entrance position of element n + 1.
Thus, the last value in {sn} is the exit position of the final beamline element.
o {γn} = Set of design relativistic parameters at each beamline element. Notethat since there is a one-to-one correspondence between this set and thenumber of beamline elements, there is one less value in {γn} as compared to{sn}.
o {σn} = Set of beam moment matrices at each location in {sn}. These values
are the main result of the RMS envelope simulation and contain all the beamstate data. For a complete description of the second-order moment matrix σ
see Section 3.
1.5

Simulation Data
The input deck is typically located in a separate file and is loaded into the SADScript
environment using the GetMAIN[InputDeckFile] function of SAD. Once the input deck
is loaded, you must extract the beamline that you wish to simulate with the call BL =ExtractBeamLine[BeamLineName]. Finally the beamline is specified by the FFS
function USE BL.
Within input deck file, along with the definition of the lattice under simulation, there
should be the following parameter statements
MASS= 0.939294 GEV;CHARGE = -1;
MOMENTUM = 0.610624 GEV;
which specify the mass, charge, and initial momentum of the beam particle, respectively(we have used the example of an 181 MeV proton). These parameters are used as design
parameters for the machine lattice.
The other beam parameters are specified as arguments to the main simulation
function ScheffSimulate[K0, σ0]. Specifically, these arguments are the initial
generalized beam perveance K0, the initial moment matrix σ0 = 〈zzT〉 (the other arguments
are optional numerical tuning parameters). For a complete description of theseparameters see Section 1.3.
The output data of ScheffSimulate[] is not stored in the SAD environment.
The returned values of element positions {sn}, gammas at each element {γn}, and momentmatrix at each element {σn} are all independent of SAD.

5

1.6

Modifying Beamline Parameters
ScheffSimulate[] retrieves the beamline parameters directly from the SAD
environment. So it is possible to change a beamline element parameter directly then
immediately re-run the simulation. For example, we may use a call to the SADScript
environment such as
ELEMENT[“K1”, “QuadHor02”]= PI/3;

which will set the parameter “K1” of the element named “QuadHor02” to PI/3. The
simulation can then be re-run with the new element setting. However, except for MASS,CHARGE, and MOMENTUM, ScheffSimulate[] ignores all other beam input
parameters in the SAD environment. The other beam parameters must be set using thevalues of initial generalized beam perveance K0 and initial moment matrix σ0 (asarguments to ScheffSimulate[]).

2 Perveance and Simulating Particles with Multiple Charges
The generalized beam perveance K is given by
K=

qI
I
q
I
1
1 1
,
=
=
3 3
2
2 3
2πε 0 mc β γ2πε 0 E R β 2 γ 3
2πε 0 mγ v

(1)where q is the unit charge, ε0 is the electric permittivity of free-space, v is the particle
velocity, c is the speed of light, β = v/c is the normalize velocity, γ is the relativistic factor,
I = qN l v is the average current with N l being the number of particles per unit length, and
ER = mc2/q is the rest energy of the beam particle in electron-Volts.
For an RF system with frequency f, bunch charge Q is related to I as
Q=

I
,
f

(2)

where we assume the bunch spacing is βλ. Substituting Eq. (2) to (1), we find
K=

qQf
Q
q
f
Q 1
f
,
=
=
3 3
2
2 3
2 3
2πε 0 mc β γ2πε 0 E R β γ
2πε 0 mγ v

(3)which means that we need information on q, m, and Q for determining K. Physicallyspeaking, the space-charge field E is determined by Q, and the space-charge forces acted
on the individual particles are determined by the product of q and E. So, it is natural tohave qQ in K.
Interpreting “m” as “m/q”, you can simulate multi-charged particles (with q other than
1 or −1), because K depends on q/m, not on q and m independently. Thus, for a particleof charge nq and the same mass m, we would substitute ER → ER/n in all the simulation
functions (e.g., ComputePerveance[], TraceToSadLongTwiss[], etc.)

63 The Second-Order Moment Matrix
3.1

Definition of the Moment Matrix
We denote the phase space coordinates of a charged particle as z = (x x’ y y’ z dp/p0)T
where x, y, and z are the offsets of the particle from the synchronous particle position in
the three orthogonal directions, x ′ = dx / ds = x& / s& = p x / p 0 , y ′ = dy / ds = y& / s& = p y / p 0 ,s is the path length parameter, the over-dot represents differentiation with respect to time,
p0 is the design momentum of the synchronous particle, and dp is the difference inmomentum from the design momentum. Every beam particle can be identified by its
position z in phase space. Consider a distribution of particles in phase space, represented
by the beam density function f : ℜ6 → ℜ. Thus, f(z)d6z is the number of particles withina differential volume d6z at location z in phase space. For any function g : ℜ6 → ℜ on
phase space, the moment of g with respect to the distribution f, denoted 〈g〉, is defined by〈g〉 ≡ ∫g(z)f(z)d6z. The matrix of second-order moments of the beam, or moment matrix σ,
is defined by σ ≡ 〈zzT〉, where zzT is the outer product of the phase space coordinates.
Expanding out this definition we see that σ appears as

σ = zz T⎛ x2⎜⎜ xx ′
⎜
⎜ xy
=⎜
′
⎜ xy
⎜ xz
⎜
⎜ xδ
⎝

xx ′

xy

xy ′

xz

x′ 2

x ′y

x ′y ′

x ′z

x ′y

y2

yy ′

yz

x ′y ′

yy ′

y′2

y ′z

x ′z

yz

y ′z

z2

x ′δ

yδ

y ′δ

zδ

xδ ⎞
⎟
x ′δ ⎟⎟yδ ⎟,y ′δ ⎟⎟zδ ⎟⎟δ 2 ⎟⎠(4)

where we have abbreviated dp/p0 as δ. This matrix is the main object of concern in thesimulation containing all the particle beam information.
3.2

Initial Moment Matrix
To initialize the simulation, we need an initial value for σ, which we denote as σ0.
This quantity is computed from the initial Twiss parameters of the beam in all three phaseplanes. Given the initial Twiss parameters (α x , β x , ε~x ) , (α y , β y , ε~y ) , (α z , β z , ε~z ) , foreach phase plane x, y, and z, respectively, the value of σ0 should be
0
0
0
0 ⎞
⎛ β xε~x − α xε~x
⎜
⎟
~
~
0
0
0
0 ⎟
⎜ − α xε x γ xε x⎜ 0β yε~y − α yε~y
0
0
0 ⎟
⎟,
σ0 = ⎜− α yε~y γ yε~y
0
0
0 ⎟
⎜ 0
⎜ 0
β zε~z − α zε~z ⎟⎟
0
0
0
⎜
⎜ 0
− α z ε~z γ z ε~z ⎟⎠
0
0
0
⎝

(5)

where the tilde indicate RMS quantities. Note that there is not coupling between phaseplanes for the initial σ matrix. During the simulation this condition may change (e.g., in7

the presence of bending magnets, misalignments, etc.).
The function
CorrelationMatrix6D[] is provided for computing this matrix given the initial
Twiss parameters.
3.3

Propagation of the Moment Matrix
In linear beam optics the phase space coordinates of a particle transform as zn+1 = Φnzn where Φn is the transfer matrix of beamline element n, zn is the phase space coordinateof the particle at the entrance of n, and zn+1 is the phase space coordinate at the exit of n.
With no space charge, we may derive the propagation equations for the moment matrix
directly from its definition and the propagation equations for the single particle. For abeamline element n the result is

σn+1 = ΦnσnΦnT(6)where Φn is the transfer matrix of element n, σn is the moment matrix at the entrance of n,and σn+1 is the moment matrix at the exit of n. This is the tenet upon which thesimulation is built. However, to include space charge effects we must determine the self
forces (from σ) then augment the dynamics σn+1 = ΦnσnΦnT accordingly.

4 SADScript Packages and Important Functions
The following is a listing of the new SADScript packages used in the RMS envelope
simulation. A description of each package is given, along with a listing of the functions
that you will typically use most in the SAD RMS envelope simulation. The functions areincluded with their arguments, shown with the standard notion found in the literature. A
complete listing of all the implemented functions, along with a brief description, is
provided in Section 5.
{ oldsad/Packages/Scheff.no {{sn},{γn},{σn}} = ScheffSimulate[K0, σ0, h0:0.01, ε:10-5, ∆h:0.05, hmax:0.0]o {{sn},{γn},{Φn}} = GetBeamLineElementData[]o SaveBeamMatrixData[file, {sn}, {γn}, {σn}]o PlotBeamBeta[{sn},{σn}]This package is the main simulation package and includes the simulation driver,
ScheffSimulate[K0, σ0, ∆s:0.01] whose arguments are the initial generalized beam
perveance K0, the initial moment matrix σ0 = 〈zzT〉 and the (optional) initial step size h0,
the (optional) solution error tolerance ε, the (optional) slack tolerance in computing new
stizes ∆h, and the (optional) maximum allowable step size hmax. (For a complete
description of ScheffSimulate[] see Section 1.3). Its function directly depends
upon one other module, MatrixFunctions.n. You will typically use the function
ComputePerveance[] from Trace3dToSad.n, to compute the beam perveance
from
the
beam
parameters,
andCorrelationMatrix6D[]
fromTwissUtility.n to compute the initial moment matrix σ0 from the initial Twiss
parameters of the beam.

8

{ oldsad/Packages/Trace3dToSad.no K = ComputePerveance[f, ER, W, Q]o {α,β,ε}SAD = TraceToSadTransTwiss[{α,β,ε}T3D]o {α,β,ε}SAD = TraceToSadLongTwiss[f, Er, W, {α,β,ε}T3D]
A convenience module for converting the parameters used in Trace3D (i.e., those foundin the Trace3D input file) to those used by SAD, and vice-versa. The two codes use
different units so it is helpful to have functions to convert between them. The most used
are listed above, where the argument f is the RF frequency (Hz), ER is the beam particle
rest energy (=mc2) in (eV), W is the beam kinetic energy (eV), Q is the bunch charge (C),and {α,β,ε} are the Twiss parameters for a single phase plane. This module has no
dependencies.

{ oldsad/Packages/TwissUtility.no σ = CorrelationMatrix6D[{α,β,ε}x, {α,β,ε}y, {α,β,ε}z]
A convenience module for working with the Twiss parameters of a beam. The most usedfunction will typically be CorrelationMatrix6D[], which builds a 6×6 momentmatrix σ from Twiss parameters in the three phase planes. The moment matrix is block
diagonal and composed of 2×2 sub-matrices from each of the uncoupled phase planes.
This module has no dependencies.
{ oldsad/Packages/MatrixFunctions.no F = MatrixLog[Φ]o Φ = MatrixExp[F]Many matrix functions are included in this module. The two most important are thematrix logarithm MatrixLog[] and matrix exponential MatrixExp[]. However,there are also functions for computing matrix inner products and norms. The matrix
logarithm function is designed for symplectic arguments and is not robust in general.
The function may no converge for arguments sufficiently far from the identity matrix,
even though the argument does, indeed, have a logarithm.

5 Function List and Calling Hierarchy
The follow is list of all the new functions added to the SAD Environment, both as
SADScript and compiled code. A brief description of the function is provided in themargin. For functions implemented in SADScript, a full description is given in thesource code, including the arguments and the returned values. All SADScript modules
are located under the oldsad/Packages/ directory in the repository.
5.1 Module Scheff.n
ScheffSimulate
ScheffZeroPropagate

ScheffPropagate
ScheffPropElem- Main RMS simulation function- envelope simulation w/out space charge (K=0)- envelope simulation with space charge

9StepSigmaMatrix – steps the moment matrix a specified length
CompStepSize - computes the new step size for integration
ScheffGenerator1- compute the space charge matrix in beam frame
ScheffGenerator2- compute the space charge matrix in the lab frame
ScheffDecoup- decouple the phase planes in configuration space
ScheffDecoupRotMatrix – get matrix in SO(6) zeroing largest 〈xixj〉
SaveBeamMatrixData- save solution data to disk

PlotBeamBeta- plot solution beta function in three phase planes

GetBeamlineElementData- return beamline data {{sn},{γn},{Φn}}5.2 Module Trace3dToSad.n
ComputePerveance- Compute generalized beam perveance
TraceToSadCoords- Convert Trace3D phase space coordinates to SAD coords
SadToTraceCoords- Convert SAD phase space coordinates to Trace3D coords
TraceToSadTransTwiss- Convert Trace3D transverse Twiss parameters to SAD
TraceToSadLongTwiss- Convert Trace3D longitudinal Twiss parameters to SAD
SadToTraceTransTwiss- Convert SAD transverse Twiss parameters to Trace3D
SadToTraceLongTwiss- Convert SAD longitudinal Twiss parameters to Trace3D
5.3 Module TwissUtility.n
CorrelationMatrix6D- Build full 6×6 second-order moment matrix from Twiss
CorrelationMatrixTrans- Build transverse 4×4 moment matrix from Twiss
CorrelationMatrix- Build 2×2 moment matrix for single phase plane

GetEnvelopeFromTwiss- Compute (√〈x2〉, √〈x2〉, ε) from Twiss parameters

GetTwissFromEnvelope- Compute (α, β, ε) from the beam envelope

5.4 Module MatrixFunctions.n
MatrixLog- Computes the log of a matrix close to identity
MatrixInnerProd2- Computes the l2 inner product of a matrix in ℜn×n

MatrixZassen
MatrixCommutator- Zassenhaus formula of two matrices to 2nd order- Computes the commutator of two matrices

MatrixLogTaylor- Taylor expansion of matrix logarithm function

MatrixExp
MatrixNorm2- Compute the exponential of a matrix- Compute the l2 norm (squared) of a matrix in ℜn×n

10

MatrixInnerProd2
MatrixExpTaylor- Induces the l2 norm

- Taylor expansion of the matrix exponent. function

5.5 Internal Additions to the SAD Interpreter
EllipticRd(x,y,z) - Computes the Carlson elliptic integral of the second kind given by theformulaR D ( x, y , z ) ≡

dt
3 ∞
1/ 2
∫
0
2 (t + x ) (t + y )1 / 2 (t + z )3 / 2

116

Example Code

Below is a code excerpt from a SADScript file simulating the beam transport sectionof the J-PARC machine. This excerpt lists the essential code for a complete RMS
envelope simulation. The complete file can be found on the SAD cluster at the location/usr/users/ckallen/code/testing/ScheffTest.sad.
!!========================================!!
!! Initialize SAD!!
!!========================================FFS;! Begin SADScript!
! Load Beamline!
GetMAIN["~ckallen/J-Parc/linac/simdb-LI_L3BT01-nopmq0000.sad"];
L3BT01 = ExtractBeamLine["L3BT01all"];
USE L3BT01;!
!
!

Initialize SAD Environment – These are important!

TRPT;
INS;
CAL;

! a transport line, not a ring
! mimic an insertion device! do SAD calculation of beamline

NOCOD;
RFSW;! no closed orbit correction! this is for RF gaps, use standing wave structure$DisplayFunction = CanvasDrawer;! if you want to use PlotBeamBeta[]12!
! INITIALIZE SIMULATION!

!
! TRACE3D Parameters!
f = 324.0e6;! RF frequency (Hz)Er = 939.29432e6;
W = 181.0338e6;
XI = 30.0e-3! particle rest energy (eV)! beam kinetic energy (eV)! beam current (A)vecTwissXt3d = {-0.44117, 5.774, 1.889};
vecTwissYt3d = {0.21808, 6.4229, 1.706};
vecTwissZt3d = {0.3095, 2.0888, 466.99};!
!
!

Numerical Parametersh0 = 0.01;
errSoln = 1.0e-5;
hslack = 0.05;
hmax = 0.0;

! initial step length! solution error tolerance! adaptive step backlash tolerance! maximum step length (=0 turned off)!
! Compute SAD Parameters (Convert from Trace3D)!
Q = XI/f;! beam bunch charge (C)K0 = ComputePerveance[f, Er, W, Q];

vecTwissX = TraceToSadTransTwiss[vecTwissXt3d];
vecTwissY = TraceToSadTransTwiss[vecTwissYt3d];
vecTwissZ = TraceToSadLongTwiss[f, Er, W, vecTwissZt3d];
sig0 = CorrelationMatrix6D[vecTwissX, vecTwissY, vecTwissZ];13

!
! RUN SIMULATION!

! Run simulation{lstPos, lstGamma, lstSig} = ScheffSimulate[K0, sig0, h0, errSoln, hslack, hmax];! Store results
SaveBeamMatrixData[“RmsEnvOutput.txt”, lstPos, lstGamma, lstSig];! Look at the Results
PlotBeamBeta[lstPos, lstSig];
Exit[];

14

7 Output

15