1. OPALcycl
1.1. Introduction
OPALcycl, as one of the flavors of the OPAL framework, is a fully threedimensional parallel beam dynamics simulation program dedicated to future high intensity cyclotrons and FFAG. it tracks multiple particles which takes into account the space charge effects. For the first time in the cyclotron community, OPALcycl has the capability of tracking multiple bunches simultaneously and take into account the beambeam effects of the radially neighboring bunches (we call it neighboring bunch effects for short) by using a selfconsistent numerical simulation model.
Apart from the multiparticle simulation mode, OPALcycl also has two other serial tracking modes for conventional cyclotron machine design. One mode is the single particle tracking mode, which is a useful tool for the preliminary design of a new cyclotron. It allows one to compute basic parameters, such as reference orbit, phase shift history, stable region, and matching phase ellipse. The other one is the tune calculation mode, which can be used to compute the betatron oscillation frequency. This is useful for evaluating the focusing characteristics of a given magnetic field map.
In addition, the widely used plugin elements, including collimator, radial profile probe, septum, trimcoil field and charge stripper, are currently implemented in OPALcycl. These functionalities are very useful for designing, commissioning and upgrading of cyclotrons and FFAGs.
1.2. Tracking modes
According to the number of particles defined by the argument npart
in
beam
see Chapter [beam] , OPALcycl works in one of the following
three operation modes automatically.
1.2.1. Single Particle Tracking mode
In this mode, only one particle is tracked, either with acceleration or not. Working in this mode, OPALcycl can be used as a tool during the preliminary design phase of a cyclotron.
The 6D parameters of a single particle in the initial local frame must
be read from a file. To do this, in the OPAL input file, the command
line DISTRIBUTION
see Chapter [distribution] should be defined like
this:
Dist1: DISTRIBUTION, TYPE=fromfile, FNAME="PartDatabase.dat";
where the file PartDatabase.dat should have two lines:
1 0.001 0.001 0.001 0.001 0.001 0.001
The number in the first line is the total number of particles. In the
second line the data represents
x, p_x, y,
p_y, z, p_z
in the local
reference frame. Their units are described in
Section [variablesopalcycl].
Please don’t try to run this mode in parallel environment. You should
believe that a single processor of the 21^{st}
century is capable of doing
the single particle tracking.
1.2.2. Tune Calculation mode
In this mode, two particles are tracked, one with all data is set to
zero is the reference particle and another one is an offcentering
particle which is offcentered in both r
and
z
directions. Working in this mode, OPALcycl can
calculate the betatron oscillation frequency \nu_r
and
\nu_z
for different energies to evaluate the focusing
characteristics for a given magnetic field.
Like the single particle tracking mode, the initial 6D parameters of the two particles in the local reference frame must be read from a file, too. In the file should has three lines:
2 0.0 0.0 0.0 0.0 0.0 0.0 0.001 0.0 0.0 0.0 0.001 0.0
When the total particle number equals 2, this mode is triggered
automatically. Only the element CYCLOTRON
in the beam line is used and
other elements are omitted if any exists.
Please don’t try to run this mode in parallel environment, either.
1.2.3. Multiparticle tracking mode
In this mode, large scale particles can be tracked simultaneously, either with space charge or not, either single bunch or multibunch, either serial or parallel environment, either reading the initial distribution from a file or generating a typical distribution, either running from the beginning or restarting from the last step of a former simulation.
Because this is the main mode as well as the key part of OPALcycl, we will describe this in detail in Section [opalcycl:MultiBunch].
1.3. Variables in OPALcycl
OPALcycl uses the following canonical variables to describe the motion of particles:
 X

Horizontal position
x
of a particle in given global Cartesian coordinates [m].  PX

Horizontal canonical momentum [eV/c].
 Y

Longitudinal position
y
of a particle in global Cartesian coordinates [m].  PY

Longitudinal canonical momentum [eV/c].
 Z

Vertical position
z
of a particle in global Cartesian coordinates [m].  PZ

Vertical canonical momentum [eV/c].
The independent variable is: t [s].
NOTE: unit conversion of momentum in OPALt and OPALcycl
Convert \beta_x \gamma
[dimensionless] to [mrad],
(\beta\gamma)_{\text{ref}}=\frac{P}{m_0c}=\frac{Pc}{m_0c^2}
P_x[{mrad}]=1000\times\frac{(\beta_x\gamma)}{(\beta\gamma)_{\text{ref}}}
Convert from [eV/c] to \beta_x\gamma
[dimensionless],
(\beta_x\gamma)=\sqrt{(\frac{P_x[{eV/c}]}{m_0c}+1)^21}.
This may be deduced by analogy for the y
and
z
directions.
1.3.2. The initial distribution in the local reference frame
The initial distribution of the bunch, either read from file or
generated by a distribution generator see Chapter [distribution], is
specified in the local reference frame of the OPALcycl Cartesian
coordinate system see Section [opalcycl:canon]. At the beginning of the
run, the 6 phase space variables (x, y, z, p_x, p_y, p_z)
are transformed to the global Cartesian coordinates using the starting
coordinates r_0
(RINIT
), \phi_0
(PHIINIT
), and z_0
(ZINIT
), and the starting momenta
p_{r0}
(PRINIT
), and p_{z0}
(PZINIT
) of
the reference particle, defined in the CYCLOTRON
element
see Section [cyclotron]. Note that p_{\phi 0}
is
calculated automatically from p_{total}
,
p_{r0}
, and p_{z0}
.
\begin{aligned}
X &= x\cos(\phi_0)  y\sin(\phi_0) + r_0\cos(\phi_0) \\
Y &= x\sin(\phi_0) + y\cos(\phi_0) + r_0\sin(\phi_0) \\
Z &= z + z_0\end{aligned}
\begin{aligned}
PX &= (p_x+p_{r0})\cos(\phi_0)  (p_y+p_{\phi 0})\sin(\phi_0) \\
PY &= (p_x+p_{r0})\sin(\phi_0) + (p_y+p_{\phi 0})\cos(\phi_0) \\
PZ &= p_z + p_{z0}\end{aligned}
1.4. Field Maps
In OPALcycl, the magnetic field on the median plane is read from an ASCII type file. The field data should be stored in the cylinder coordinates frame (because the field map on the median plane of the cyclotron is usually measured in this frame).
There are two possible situations. One is the real field map on median
plane of the exist cyclotron machine using measurement equipment.
Limited by the narrow gaps of magnets, in most cases with cyclotrons,
only vertical field B_z
on the median plane
(z=0
) is measured. Since the magnetic field data off the
median plane field components is necessary for those particles with
z \neq 0
, the field need to be expanded in Z
direction. According to the approach given by Gordon and Taivassalo, by
using a magnetic potential and measured B_z
on the median
plane, at the point (r,\theta, z)
in cylindrical polar
coordinates, the third order field can be written as
\begin{aligned}
B_r(r,\theta, z) & = & z\frac{\partial B_z}{\partial r}\frac{1}{6}z^3 C_r, \\
B_\theta(r,\theta, z) & = & \frac{z}{r}\frac{\partial B_z}{\partial \theta}\frac{1}{6}\frac{z^3}{r} C_{\theta}, \\
B_z(r,\theta, z) & = & B_z\frac{1}{2}z^2 C_z, \end{aligned}
where B_z\equiv B_z(r, \theta, 0)
and
\begin{aligned}
C_r & = & \frac{\partial^{3}{B_z}}{\partial r^{3}} + \frac{1}{r}\frac{\partial^{2}{B_z}}{\partial r^{2}}  \frac{1}{r^2}\frac{\partial{B_z}}{\partial r}
+ \frac{1}{r^2}\frac{\partial^{3}{B_z}}{{\partial r}{\partial \theta^2}}  2\frac{1}{r^3}\frac{\partial^{2}{B_z}}{\partial \theta^{2}}, \\
C_{\theta} & = & \frac{1}{r}\frac{\partial^{2}{B_z}}{{\partial r}{\partial \theta}} + \frac{\partial^{3}{B_z}}{{\partial r^2}{\partial \theta}}
+ \frac{1}{r^2}\frac{\partial^{3}{B_z}}{\partial \theta^{3}}, \\
C_z & = & \frac{1}{r}\frac{\partial{B_z}}{\partial r} + \frac{\partial^{2}{B_z}}{\partial r^{2}} + \frac{1}{r^2}\frac{\partial^{2}{B_z}}{\partial \theta^{2}}.
\end{aligned}
All the partial differential coefficients are on the median plane and can be calculated by interpolation. OPALcycl uses Lagrange’s 5point formula.
The other situation is to calculate the field on the median plane or the
3D fields of the working gap for interesting region numerically by
creating 3D model using commercial software, such as TOSCA, ANSOFT and
ANSYS during the design phase of a new machine. If the field on the
median plane is calculated, the field off the median plane can be
obtained using the same expansion approach as the measured field map as
described above. However, the 3D fields of the entire working gap should
be more accurate than the expansion approach especially at the region
not so close to the median plane in Z
direction.
In the current version, we implemented the three specific type
fieldread functions Cyclotron::getFieldFromFile() of the median plane
fields. That which function is used is controlled by the parameters
TYPE
of CYCLOTRON
see Section [cyclotron] in the input file.
1.4.1. CARBONCYCL type
If TYPE=CARBONCYCL
, the program requires the B_z
data
which is stored in a sequence shown in Figure [CYCLField].
We need to add 6 parameters at the header of a plain B_z
[kG] data file, namely, r_{min}
[mm],
\Delta r
[mm], \theta_{min}
[^{\circ}
],
\Delta \theta
[^{\circ}
], N_\theta
(total data
number in each arc path of azimuthal direction) and N_r
(total path number along radial direction). If \Delta r
or
\Delta \theta
is decimal, one can set its negative
opposite number. For instance, if
\Delta \theta = \frac{1}{3}{^{\circ}}
, the fourth line of
the header should be set to 3.0. Example showing the above explained
format:
3.0e+03 10.0 0.0 3.0 300 161 1.414e03 3.743e03 8.517e03 1.221e02 2.296e02 3.884e02 5.999e02 8.580e02 1.150e01 1.461e01 1.779e01 2.090e01 2.392e01 2.682e01 2.964e01 3.245e01 3.534e01 3.843e01 4.184e01 4.573e01 ......
1.4.2. CYCIAE type
If TYPE=CYCIAE
, the program requires data format given by ANSYS10.0.
This function is originally for the 100MeV cyclotron of CIAE, whose
isochronous fields is numerically computed by by ANSYS. The median plane
fields data is output by reading the APDL (ANSYS Parametric Design
Language) script during the postprocessing phase (you may need to do
minor changes to adapt your own cyclotron model):
/post1 resume,solu,db csys,1 nsel,s,loc,x,0 nsel,r,loc,y,0 nsel,r,loc,z,0 PRNSOL,B,COMP CSYS,1 rsys,1 *do,count,0,200 path,cyc100_Ansys,2,5,45 ppath,1,,0.01*count,0,,1 ppath,2,,0.01*count/sqrt(2),0.01*count/sqrt(2),,1 pdef,bz,b,z paget,data,table *if,count,eq,0,then /output,cyc100_ANSYS,dat *STATUS,data,,,5,5 /output *else /output,cyc100_ANSYS,dat,,append *STATUS,data,,,5,5 /output *endif *enddo finish
By running this in ANSYS, you can get a fields file with the name
cyc100_ANSYS.data. You need to put 6 parameters at the header of the
file, namely, r_{min}
[mm], \Delta r
[mm],
\theta_{min}
[^{\circ}
], \Delta \theta
[^{\circ}
],
N_\theta
(total data number in each arc path of azimuthal
direction) and N_r
(total path number along radial
direction). If \Delta r
or \Delta \theta
is
decimal,one can set its negative opposite number. This is useful is the
decimal is unlimited. For instance,if
\Delta \theta = \frac{1}{3} {^{\circ}}
, the fourth line of
the header should be 3.0. In a word, the fields are stored in the
following format:
0.0 10.0 0.0e+00 1.0e+00 90 201 PARAMETER STATUS DATA ( 336 PARAMETERS DEFINED) (INCLUDING 17 INTERNAL PARAMETERS) LOCATION VALUE 1 5 1 0.537657876 2 5 1 0.538079473 3 5 1 0.539086731 ...... 44 5 1 0.760777286 45 5 1 0.760918663 46 5 1 0.760969074 PARAMETER STATUS DATA ( 336 PARAMETERS DEFINED) (INCLUDING 17 INTERNAL PARAMETERS) LOCATION VALUE 1 5 1 0.704927299 2 5 1 0.705050993 3 5 1 0.705341275 ......
1.4.3. BANDRF type
If TYPE=BANDRF
, the program requires the B_z
data format
which is same as CARBONCYCL
. But with BANDRF
type, the program can
also read in the 3D electric field(s). For the detail about its usage,
please see Section [cyclotron].
1.4.4. Default PSI format
If the value of TYPE
is other string rather than above mentioned, the
program requires the data format like PSI format field file ZYKL9Z.NAR
and SO3AV.NAR, which was given by the measurement. We add 4 parameters
at the header of the file, namely, r_{min}
[mm],
\Delta r
[mm], \theta_{min}
[^{\circ}
],
\Delta \theta
[^{\circ}
], If \Delta r
or
\Delta \theta
is decimal,one can set its negative opposite
number. This is useful is the decimal is unlimited. For instance,if
\Delta \theta = \frac{1}{3}{^{\circ}}
, the fourth line of
the header should be 3.0.
1900.0 20.0 0.0 3.0 LABEL=S03AV CFELD=FIELD NREC= 141 NPAR= 3 LPAR= 7 IENT= 1 IPAR= 1 3 141 135 30 8 8 70 LPAR= 1089 IENT= 2 IPAR= 2 0.100000000E+01 0.190000000E+04 0.200000000E+02 0.000000000E+00 0.333333343E+00 0.506500015E+02 0.600000000E+01 0.938255981E+03 0.100000000E+01 0.240956593E+01 0.282477260E+01 0.340503168E+01 0.419502926E+01 0.505867147E+01 0.550443363E+01 0.570645094E+01 0.579413509E+01 0.583940887E+01 0.586580372E+01 0.588523722E+01 ......
1.4.5. 3D fieldmap
It is additionally possible to load 3D fieldmaps for tracking through
OPALcycl. 3D fieldmaps are loaded by sequentially adding new field
elements to a line, as is done in OPALt. It is not possible to add RF
cavities while operating in this mode. In order to define ring
parameters such as initial ring radius a RINGDEFINITION
type is loaded
into the line, followed by one or more SBEND3D
elements.
ringdef: RINGDEFINITION, HARMONIC_NUMBER=6, LAT_RINIT=2350.0, LAT_PHIINIT=0.0, LAT_THETAINIT=0.0, BEAM_PHIINIT=0.0, BEAM_PRINIT=0.0, BEAM_RINIT=2266.0, SYMMETRY=4.0, RFFREQ=0.2; triplet: SBEND3D, FMAPFN="fdftoscafieldmap.table", LENGTH_UNITS=10., FIELD_UNITS=1e4; l1: Line = (ringdef, triplet, triplet);
The fieldmap with file name fdftoscafieldmap.table is loaded, which is a file like
422280 422280 422280 1 1 X [LENGU] 2 Y [LENGU] 3 Z [LENGU] 4 BX [FLUXU] 5 BY [FLUXU] 6 BZ [FLUXU] 0 194.01470 0.0000000 80.363520 0.68275932346E07 5.3752492577 0.28280706805E07 194.36351 0.0000000 79.516210 0.42525693524E07 5.3827955117 0.17681348191E07 194.70861 0.0000000 78.667380 0.19766168358E07 5.4350026348 0.82540823165E08 <continues>
The header parameters are ignored  user supplied parameters
LENGTH_UNITS
and FIELD_UNITS
are used. Fields are supplied on points
in a grid in (r, y, \phi)
. Positions and field elements
are specified by Cartesian coordinates (x, y, z)
.
1.4.6. User’s own fieldmap
You should revise the function or write your own function according to
the instructions in the code to match your own field format if it is
different to above types. For more detail about the parameters of
CYCLOTRON
, please refer to Section [cyclotron].
1.5. RF field
1.5.1. Read RF voltage profile
The RF cavities are treated as straight lines with infinitely narrow
gaps and the electric field is a \delta
function plus a
transit time correction. the twogap cavity is treated as two
independent singlegap cavities. the spiral gap cavity is not
implemented yet. For more detail about the parameters of cyclotron
cavities, see Section [cavitycycl].
The voltage profile of a cavity gap is read from ASCII file. Here is an example:
6 0.00 0.15 2.40 0.20 0.65 2.41 0.40 0.98 0.66 0.60 0.88 1.59 0.80 0.43 2.65 1.00 0.05 1.71
The number in the first line means 6 sample points and in the following lines the three values represent the normalized distance to the inner edge of the cavity, the normalized voltage and its derivative respectively.
1.5.2. Read 3D RF fieldmap
The 3D RF fieldmap can be read from H5hut type file. This is useful for modeling the central region electric fields which usually has complicate shapes. For the detail about its usage, please see Section [cyclotron].
Please note that in this case, the E field is treated as a part of
CYCLOTRON
element, rather than a independent RFCAVITY
element.
1.6. Particle Tracking and Acceleration
The precision of the tracking methods is vital for the entire simulation
process, especially for long distance tracking jobs. OPALcycl uses a
fourth order RungeKutta algorithm and the second order LeapFrog
scheme. The fourth order RungeKutta algorithm needs four external
magnetic field evaluation in each time step \tau
. During
the field interpolation process, for an arbitrary given point the code
first interpolates Formula B_z
for its counterpart on the
median plane and then expands to this give point using
Equation [Bfield].
After each time step i
, the code detects whether the
particle crosses any one of the RF cavities during this step. If it
does, the time point t_c
of crossing is calculated and the
particle return back to the start point of step i
. Then
this step is divided into three substeps: first, the code tracks this
particle for t_1 = \tau  (t_ct_{i1})
; then it
calculates the voltage and adds momentum kick to the particle and
refreshes its relativistic parameters \beta
and
\gamma
; and then tracks it for
t_2 = \tau  t_1
.
1.7. Space Charge
OPALcycl uses the same solvers as OPALt to calculate the space charge effects see Chapter [fieldsolver].
Typically, the space charge field is calculated once per time step. This is no surprise for the second order BorisBuneman time integrator (leapfroglike scheme) which has per default only one force evaluation per step. The fourth order RungeKutta integrator keeps the space charge field constant for one step, although there are four external field evaluations. There is an experimental multipletimestepping (MTS) variant of the BorisBuneman/leapfrogmethod, which evaluates space charge only every N^th step, thus greatly reducing computation time while usually being still accurate enough.
1.8. Multibunch Mode
The neighboring bunches problem is motivated by the fact that for high intensity cyclotrons with small turn separation, single bunch space charge effects are not the only contribution. Along with the increment of beam current, the mutual interaction of neighboring bunches in radial direction becomes more and more important, especially at large radius where the distances between neighboring bunches get increasingly small and even they can overlap each other. One good example is PSI 590MeV Ring cyclotron with a current of about 2mA in CW operation and the beam power amounts to 1.2MW. An upgrade project for Ring is in process with the goal of 1.8MW CW on target by replacing four old aluminum resonators by four new copper cavities with peak voltage increasing from about 0.7MW to above 0.9MW. After upgrade, the total turn number is reduced from 200 turns to less than 170 turns. Turn separation is increased a little bit, but still are at the same order of magnitude as the radial size of the bunches. Hence once the beam current increases from 2mA to 3mA, the mutual space charge effects between radially neighboring bunches can have significant impact on beam dynamics. In consequence, it is important to cover neighboring bunch effects in the simulation to quantitatively study its impact on the beam dynamics.
In OPALcycl, we developed a new fully consistent algorithm of
multibunch simulation. We implemented two working modes, namely ,
AUTO
and FORCE
. In the first mode, only a single bunch is tracked
and accelerated at the beginning, until the radial neighboring bunch
effects become an important factor to the bunches’ behavior. Then the
new bunches will be injected automatically to take these effects into
account. In this way, we can save time and memory, and more important,
we can get higher precision for the simulation in the region where
neighboring bunch effects are important. In the other mode, multibunch
simulation starts from the injection points.
In the space charge calculation for multibunch, the computation region covers all bunches. Because the energy of the bunches is quite different, it is inappropriate to use only one particle rest frame and a single Lorentz transformation any more. So the particles are grouped into different energy bins and in each bin the energy spread is relatively small. We apply Lorentz transforming, calculate the space charge fields and apply the back Lorentz transforming for each bin separately. Then add the field data together. Each particle has a ID number to identify which energy bin it belongs to.
1.9. Input
All the three working modes of OPALcycl use an input file written in
mad
language which will be described in detail in the following
chapters.
For the Tune Calculation mode, one additional auxiliary file with the following format is needed.
72.000 2131.4 0.240 74.000 2155.1 0.296 76.000 2179.7 0.319 78.000 2204.7 0.309 80.000 2229.6 0.264 82.000 2254.0 0.166 84.000 2278.0 0.025
In each line the three values represent energy E
, radius
r
and P_r
for the SEO (Static Equilibrium
Orbit) at starting point respectively and their units are MeV, mm and
mrad.
A bash script tuning.sh is shown on the next page, to execute OPALcycl for tune calculations.
#!/bin/bash rm rf tempfile rm f plotdata rm f tuningresult exec 6<FIXPO_SEO N="260" j="1" while read u 6 E1 r pr # read in Energy, initil R and intial Pr of each # SEO from FIXPO output do rm rf tempfile echo n j = echo "$j" echo n energy= > tempfile echo n "$E1" >> tempfile echo ";" >> tempfile echo n r= >> tempfile echo n "$r" >> tempfile echo ";" >> tempfile echo n pr= >> tempfile echo n "$pr" >> tempfile echo ";" >> tempfile # execute OPAL to calculate tuning frquency and store opal testcycl.in commlib mpi info 0  \ grep "Max" >>tuningresult j=$[$j+1] done exec 6<& rm rf tempfile # post porcess exec 8<tuningresult rm f plotdata i="0" while [ $i lt $N ] do read u 8 a b ur1 d read u 8 aa bb uz1 dd echo "$ur1 $uz1" >>plotdata i=$[$i+1] done exec 8<& rm f tuningresult
To start execution, just run tuning.sh which uses the input file testcycl.in and the auxiliary file FIXPO SEO_. The output file is plotdata from which one can plot the tune diagram.
1.10. Output
Single Particle Tracking mode
The intermediate phase space data is stored in an ASCII file which can
be used to the plot the orbit. The file’s name is combined by input file
name (without extension) and trackOrbit.dat. The data are stored in
the global Cartesian coordinates. The frequency of the data output can
be set using the option SPTDUMPFREQ
of OPTION
statement
see Section [option]
The phase space data per STEPSPERTURN
(a parameter in the TRACK
command) steps is stored in an ASCII file. The file’s name is a
combination of input file name (without extension) and
afterEachTurn.dat. The data is stored in the global cylindrical
coordinate system. Please note that if the field map is ideally
isochronous, the reference particle of a given energy take exactly one
revolution in STEPPERTURN
steps; Otherwise, the particle may not go
through a full 360^ in STEPPERTURN
steps.
There are 3 ASCII files which store the phase space data around
0
, \pi/8
and \pi/4
azimuths.
Their names are the combinations of input file name (without extension)
and Angle0.dat, Angle1.dat and Angle2.dat respectively. The
data is stored in the global cylindrical coordinate system, which can be
used to check the property of the closed orbit.
Tune calculation mode
The tunes \nu_r
and \nu_z
of each energy are
stored in a ASCII file with name tuningresult.
Multiparticle tracking mode
The intermediate phase space data of all particles and some interesting
parameters, including RMS envelop size, RMS emittance, external field,
time, energy, length of path, number of bunches and tracking step, are
stored in the H5hut fileformat [bib:howison2010] and can be analyzed
using the H5root [bib:schietinger]. The frequency of the data output can
be set using the PSDUMPFREQ
option of OPTION
statement
see Section [option]. The file is named like the input file but the
extension is .h5.
The intermediate phase space data of central particle (with ID of 0) and
an offcentering particle (with ID of 1) are stored in an ASCII file.
The file’s name is combined by the input file name (without extension)
and trackOrbit.dat. The frequency of the data output can be set using
the SPTDUMPFREQ
option of OPTION
statement see Section [option].
1.11. Matched Distribution
In order to run matched distribution simulation one has to specify a
periodic accelerator. The function call also needs the symmetry of the
machine as well as a field map. The user then specifies the emittance
\pi{mmmrad}
.
/* * specify periodic accelerator */ l1 = ... /* * try finding a matched distribution */ Dist1:DISTRIBUTION, TYPE=GAUSSMATCHED, LINE=l1, FMAPFN=..., MAGSYM=..., EX = ..., EY = ..., ET = ...;
1.11.1. Example
Simulation of the PSI Ring Cyclotron at 580MeV and current 2.2mA. The
program finds a matched distribution followed by a bunch initialization
according to the matched covariance matrix.
The matched distribution algorithm works with normalized emittances,
i.e. normalized by the lowest energy of the machine. The printed
emittances, however, are the geometric emittances. In addition, it has
to be paid attention that the computation is based on
(x,x',y,y',z,\delta)
instead of
(x,p_{x},y,p_{y},z,p_{z})
. Since the particles are
represented in the latter coordinate system, the corresponding
transformation has to be applied to obtain the rms emittances that are
given in the output.
Input file
OPTION, ECHO = FALSE; OPTION, PSDUMPFREQ = 1; OPTION, SPTDUMPFREQ = 1; OPTION, PSDUMPEACHTURN = FALSE; OPTION, PSDUMPLOCALFRAME = FALSE; OPTION, REPARTFREQ = 20; OPTION, CZERO = FALSE; OPTION, TELL = TRUE; OPTION, VERSION = 10600; Title,string="OPALcycl: test matched distribution PSI Ring"; REAL Edes=.580; REAL gamma=(Edes+PMASS)/PMASS; REAL beta=sqrt(1(1/gamma^2)); REAL gambet=gamma*beta; REAL P0 = gamma*beta*PMASS; REAL brho = (PMASS*1.0e9*gambet) / CLIGHT; //value,{gamma,brho,Edes,beta,gambet}; REAL phi01=139.4281; REAL phi02=phi01+180.0; REAL phi04=phi01; REAL phi05=phi01+180.0; REAL phi03=phi01+10.0; REAL volt1st=0.9; REAL volt3rd=0.9*4.0*0.112; REAL turns = 0; REAL nstep = 0; REAL frequency=50.650; REAL frequency3=3.0*frequency; ring: Cyclotron, TYPE="RING", CYHARMON=6, PHIINIT=0.0, PRINIT=0.000174, RINIT=2130.0, SYMMETRY=8.0, RFFREQ=frequency, FMAPFN="s03av.nar", FMLOWE=72.0, FMHIGHE=590.0; rf0: RFCavity, VOLT=volt1st, FMAPFN="Cav1.dat", TYPE="SINGLEGAP", FREQ=frequency, RMIN = 1900.0, RMAX = 4500.0, ANGLE=35.0, PDIS = 416.0, GAPWIDTH = 220.0, PHI0=phi01; rf1: RFCavity, VOLT=volt1st, FMAPFN="Cav1.dat", TYPE="SINGLEGAP", FREQ=frequency, RMIN = 1900.0, RMAX = 4500.0, ANGLE=125.0, PDIS = 416.0, GAPWIDTH = 220.0, PHI0=phi02; rf2: RFCavity, VOLT=volt3rd, FMAPFN="Cav3.dat", TYPE="SINGLEGAP", FREQ=frequency3,RMIN = 1900.0, RMAX = 4500.0, ANGLE=170.0, PDIS = 452.0, GAPWIDTH = 250.0, PHI0=phi03; rf3: RFCavity, VOLT=volt1st, FMAPFN="Cav1.dat", TYPE="SINGLEGAP", FREQ=frequency, RMIN = 1900.0, RMAX = 4500.0, ANGLE=215.0, PDIS = 416.0, GAPWIDTH = 220.0, PHI0=phi04; rf4: RFCavity, VOLT=volt1st, FMAPFN="Cav1.dat", TYPE="SINGLEGAP", FREQ=frequency, RMIN = 1900.0, RMAX = 4500.0, ANGLE=305.0, PDIS = 416.0, GAPWIDTH = 220.0, PHI0=phi05; l1: LINE = (ring,rf0,rf1,rf2,rf3,rf4); Dist1:DISTRIBUTION, TYPE=GAUSSMATCHED, LINE=l1, FMAPFN="ring590_bfld.dat", MAGSYM= 8, EX = 1.0e06, EY = 1.0e06, ET = 1.0e06; Fs1:FIELDSOLVER, FSTYPE=FFT, MX=32, MY=32, MT=16, PARFFTX=true, PARFFTY=true, PARFFTT=true, BCFFTX=open, BCFFTY=open, BCFFTT=open; beam1: BEAM, PARTICLE=PROTON, pc=P0, NPART=1e7, BCURRENT=2.2E3, CHARGE=1.0, BFREQ= frequency; Select, Line=l1; TRACK,LINE=l1, BEAM=beam1, MAXSTEPS=nstep*turns, STEPSPERTURN=nstep, TIMEINTEGRATOR="RK4"; run, method = "CYCLOTRONT", beam=beam1, fieldsolver=Fs1, distribution=Dist1; endtrack; Stop;
Output
OPAL> Current settings of options: OPAL> OPTION,ECHO=FALSE,INFO=TRUE,TRACE=FALSE,VERIFY=FALSE,WARN=TRUE, OPAL> SEED=1.23457e+08,TELL=TRUE,PSDUMPFREQ=1,STATDUMPFREQ=10, OPAL> PSDUMPEACHTURN=FALSE,PSDUMPLOCALFRAME=FALSE,SPTDUMPFREQ=1, OPAL> REPARTFREQ=20,REBINFREQ=100,SCSOLVEFREQ=1,MTSSUBSTEPS=1, OPAL> REMOTEPARTDEL=1,SCAN=FALSE,RHODUMP=FALSE,EBDUMP=FALSE, OPAL> CSRDUMP=FALSE,AUTOPHASE=0,PPDEBUG=FALSE,SURFDUMPFREQ=1, OPAL> NUMBLOCKS=0,RECYCLEBLOCKS=0,NLHS=1,CZERO=FALSE,RNGTYPE="RANDOM", OPAL> SCHOTTKYCORR=FALSE,SCHOTTKYRENO=1,ENABLEHDF5=TRUE,ASCIIDUMP=FALSE, OPAL> BOUNDPDESTROYFQ=10,BEAMHALOBOUNDARY=0; OPAL> OPAL> *  OPAL> * About to find closed orbit and matched distribution OPAL> * I= 2.2 (mA) E= 580 (MeV) OPAL> * EX= 1e06* EY= 1e06* ET= 1e06* FMAPFN ring590_bfld.dat OPAL> * FMSYM= 8* HN= 6* PHIINIT= 0 OPAL> *  OPAL> * Converged (Ex, Ey, Ez) = (0.313794, 0.313794, 0.313794) pi mm mrad for E= 580 (MeV) OPAL> * SigmaMatrix OPAL> 3.914 & 0.6486 & 0 & 0 & 0.4542 & 1.362 \\ OPAL> 0.6486 & 0.6396 & 0 & 0 & 0.7265 & 0.2685 \\ OPAL> 0 & 0 & 6.239 & 1.198 & 0 & 0 \\ OPAL> 0 & 0 & 1.198 & 0.3859 & 0 & 0 \\ OPAL> 0.4542 & 0.7265 & 0 & 0 & 2.363 & 0.1752 \\ OPAL> 1.362 & 0.2685 & 0 & 0 & 0.1752 & 0.8944 \\ OPAL> * m before gsl_linalg_cholesky_decomp OPAL> * r(0, : ) = 1 & 0.4099 & 0 & 0 & 0.1493 & 0.7277 \\ OPAL> * r(1, : ) = 0.4099 & 1 & 0 & 0 & 0.5909 & 0.355 \\ OPAL> * r(2, : ) = 0 & 0 & 1 & 0.7722 & 0 & 0 \\ OPAL> * r(3, : ) = 0 & 0 & 0.7722 & 1 & 0 & 0 \\ OPAL> * r(4, : ) = 0.1493 & 0.5909 & 0 & 0 & 1 & 0.1205 \\ OPAL> * r(5, : ) = 0.7277 & 0.355 & 0 & 0 & 0.1205 & 1 \\ OPAL> * m after gsl_linalg_cholesky_decomp OPAL> * r(0, : ) = 1 & 0 & 0 & 0 & 0 & 0 \\ OPAL> * r(1, : ) = 0.4099 & 0.9121 & 0 & 0 & 0 & 0 \\ OPAL> * r(2, : ) = 0 & 0 & 1 & 0 & 0 & 0 \\ OPAL> * r(3, : ) = 0 & 0 & 0.7722 & 0.6354 & 0 & 0 \\ OPAL> * r(4, : ) = 0.1493 & 0.715 & 0 & 0 & 0.683 & 0 \\ OPAL> * r(5, : ) = 0.7277 & 0.06211 & 0 & 0 & 0.08229 & 0.6781 \\ OPAL> * Mass of simulation particle= 25.44 GeV/c^2 OPAL> * Charge of simulation particle= 4.344e18 [C] OPAL> * ********************************************************************************** OPAL> * Selected Tracking Method == CYCLOTRONT, NEW TRACK OPAL> * ********************************************************************************** OPAL> * Number of neighbour bunches= 1 OPAL> * DT = 1e12 OPAL> * MAXSTEPS = 0 OPAL> * Phase space dump frequency = 1 OPAL> * Statistics dump frequency = 10 w.r.t. the time step. OPAL> * ********************************************************************************** OPAL> * ************* D I S T R I B U T I O N ******************************************** OPAL> * OPAL> * Distribution type: MATCHEDGAUSS OPAL> * SIGMAX = 0.001978 [m] OPAL> * SIGMAY = 0.002498 [m] OPAL> * SIGMAZ = 0.001537 [m] OPAL> * SIGMAPX = 4.129e05 [Beta Gamma] OPAL> * SIGMAPY = 3.639e05 [Beta Gamma] OPAL> * SIGMAPZ = 4.49e05 [Beta Gamma] OPAL> * AVRGPZ = 0 [Beta Gamma] OPAL> * CORRX = 0.4099 OPAL> * CORRY = 0.7722 OPAL> * CORRZ = 0.1205 OPAL> * R61 = 0.7277 OPAL> * R62 = 0.355 OPAL> * R51 = 0.1493 OPAL> * R52 = 0.5909 OPAL> * CUTOFFX = 3 [units of SIGMAX] OPAL> * CUTOFFY = 3 [units of SIGMAY] OPAL> * CUTOFFZ = 3 [units of SIGMAZ] OPAL> * CUTOFFPX = 3 [units of SIGMAPX] OPAL> * CUTOFFPY = 3 [units of SIGMAPY] OPAL> * CUTOFFPZ = 3 [units of SIGMAPY] OPAL> * OPAL> * Distribution is injected. OPAL> * OPAL> * ********************************************************************************** OPAL> OPAL> * ************* B E A M ************************************************************ OPAL> * BEAM BEAM1 OPAL> * PARTICLE PROTON OPAL> * CURRENT 0.0022 A OPAL> * FREQUENCY 50.65 MHz OPAL> * CHARGE +e * 1 OPAL> * REST MASS 0.9383 GeV OPAL> * MOMENTUM 1.194 OPAL> * NPART 1e+07 OPAL> * ********************************************************************************** OPAL> OPAL> * ************* F I E L D S O L V E R ********************************************** OPAL> * FIELDSOLVER FS1 OPAL> * TYPE FFT OPAL> * NPROCESSORS 1 OPAL> * MX 32 OPAL> * MY 32 OPAL> * MT 16 OPAL> * BBOXINCR 2 OPAL> * GRRENSF INTEGRATED OPAL> * XDIM parallel OPAL> * YDIM parallel OPAL> * Z(T)DIM parallel OPAL> OPAL> * ********************************************************************************** /* * Adding Cyclotron skipped */ OPAL> * *********************** Bunch information in local frame: ************************ OPAL> * ************** B U N C H ********************************************************* OPAL> * NP = 10000000 OPAL> * Qtot = 4.34353e02 [nC] Qi = 4.34353e09 [nC] OPAL> * Ekin = 5.80000e+02 [MeV] dEkin = 0.00000e+00 [MeV] OPAL> * rmax = ( 6.17284e03 , 7.79039e03 , 4.79664e03 ) [m] OPAL> * rmin = ( 6.17191e03 , 7.78887e03 , 4.79646e03 ) [m] OPAL> * rms beam size = ( 1.91510e03 , 2.41332e03 , 1.50548e03 ) [m] OPAL> * rms momenta = ( 3.99780e05 , 3.51484e05 , 4.37486e05 ) [beta gamma] OPAL> * mean position = ( 1.56300e13 , 2.13395e17 , 8.71013e21 ) [m] OPAL> * mean momenta = ( 5.72459e17 , 1.27218e+00 , 2.15345e17 ) [beta gamma] OPAL> * rms emittance = ( 5.52530e08 , 4.34281e08 , 5.13673e08 ) (not normalized) OPAL> * rms correlation = ( 3.96345e01 , 7.58800e01 , 1.24755e01 ) OPAL> * hr = ( 3.98218e04 , 5.02557e04 , 6.39540e04 ) [m] OPAL> * dh = 2.00000e+00 [%] OPAL> * t = 0.00000e+00 [s] dT = inf [s] OPAL> * spos = 8.71013e21 [m] OPAL> * ********************************************************************************** OPAL> OPAL> OPAL> * Phase space dump 0 (global frame) at integration step 0 OPAL> * T = 0 ns, Live Particles: 10000000 OPAL> * E = 580 MeV, beta * gamma = 1.2722 OPAL> * Bunch position: R = 2130 mm, Theta = 2.4758e05 Deg, Z = 5.4519e05 mm OPAL> * Local Azimuth = 0.0078362 Deg, Local Elevation = 0.0028597 Deg OPAL> OPAL> * *********************** Bunch information in global frame: *********************** OPAL> * ************** B U N C H ********************************************************* OPAL> * NP = 10000000 OPAL> * Qtot = 4.34353e02 [nC] Qi = 4.34353e09 [nC] OPAL> * Ekin = 5.80000e+02 [MeV] dEkin = 3.22716e02 [MeV] OPAL> * rmax = ( 6.17284e03 , 7.79039e03 , 4.79664e03 ) [m] OPAL> * rmin = ( 6.17191e03 , 7.78887e03 , 4.79646e03 ) [m] OPAL> * rms beam size = ( 1.91510e03 , 2.41332e03 , 1.50548e03 ) [m] OPAL> * rms momenta = ( 3.99780e05 , 3.51500e05 , 4.37486e05 ) [beta gamma] OPAL> * mean position = ( 2.13000e+00 , 9.20405e07 , 5.45188e08 ) [m] OPAL> * mean momenta = ( 1.73993e04 , 1.27218e+00 , 6.34947e05 ) [beta gamma] OPAL> * rms emittance = ( 5.52530e08 , 4.34328e08 , 5.13673e08 ) (not normalized) OPAL> * rms correlation = ( 3.96344e01 , 7.58765e01 , 1.24755e01 ) OPAL> * hr = ( 3.98218e04 , 5.02557e04 , 6.39540e04 ) [m] OPAL> * dh = 2.00000e+00 [%] OPAL> * t = 0.00000e+00 [s] dT = inf [s] OPAL> * spos = 5.45188e08 [m] OPAL> * ********************************************************************************** OPAL> OPAL> * Beginning of this run is at t = 0 [ns] OPAL> * The time step is set to dt = inf [ns] OPAL> * Single particle trajectory dump frequency is set to 1 OPAL> * The frequency to solve space charge fields is set to 1 OPAL> * The repartition frequency is set to 20 OPAL> OPAL> *  Start tracking  * OPAL> OPAL> *  DONE TRACKING PARTICLES  * OPAL> * Final energy of reference particle = 580 [MeV] OPAL> * Total phase space dump number(includes the initial distribution) = 1 OPAL> * One can restart simulation from the last dump step (restart 0) OPAL> * OPAL> * Finished during turn 1 (0 turns completed) OPAL> * Cave: Turn number is not correct for restart mode OPAL> OPAL> * *********************** Bunch information in global frame: *********************** OPAL> * ************** B U N C H ********************************************************* OPAL> * NP = 10000000 OPAL> * Qtot = 4.34353e02 [nC] Qi = 4.34353e09 [nC] OPAL> * Ekin = 5.80000e+02 [MeV] dEkin = 3.22716e02 [MeV] OPAL> * rmax = ( 6.17284e03 , 7.79039e03 , 4.79664e03 ) [m] OPAL> * rmin = ( 6.17191e03 , 7.78887e03 , 4.79646e03 ) [m] OPAL> * rms beam size = ( 1.91510e03 , 2.41332e03 , 1.50548e03 ) [m] OPAL> * rms momenta = ( 3.99780e05 , 3.51500e05 , 4.37486e05 ) [beta gamma] OPAL> * mean position = ( 2.13000e+00 , 9.20405e07 , 5.45188e08 ) [m] OPAL> * mean momenta = ( 1.73993e04 , 1.27218e+00 , 6.34947e05 ) [beta gamma] OPAL> * rms emittance = ( 5.52530e08 , 4.34328e08 , 5.13673e08 ) (not normalized) OPAL> * rms correlation = ( 3.96344e01 , 7.58765e01 , 1.24755e01 ) OPAL> * hr = ( 3.98218e04 , 5.02557e04 , 6.39540e04 ) [m] OPAL> * dh = 2.00000e+00 [%] OPAL> * t = 0.00000e+00 [s] dT = inf [s] OPAL> * spos = 5.45188e08 [m] OPAL> * ********************************************************************************** OPAL> OPAL> *  * OPAL> * Finalizing i.e. write data and close files : OPAL> * Finalize cyclotron OPAL> *  *