Jochem Snuverink1. OPAL-cycl
1.1. Introduction
OPAL-cycl, as one of the flavors of the OPAL framework, is a fully three-dimensional 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, OPAL-cycl has the capability of tracking multiple bunches simultaneously and take into account the beam-beam effects of the radially neighboring bunches (we call it neighboring bunch effects for short) by using a self-consistent numerical simulation model.
Apart from the multi-particle simulation mode, OPAL-cycl 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, trim-coil field and charge stripper, are currently implemented in OPAL-cycl. 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] , OPAL-cycl 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, OPAL-cycl 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 off-centering
particle which is off-centered in both r and
z directions. Working in this mode, OPAL-cycl 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. Multi-particle tracking mode
In this mode, large scale particles can be tracked simultaneously, either with space charge or not, either single bunch or multi-bunch, 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 OPAL-cycl, we will describe this in detail in Section [opalcycl:MultiBunch].
1.3. Variables in OPAL-cycl
OPAL-cycl uses the following canonical variables to describe the motion of particles:
- X
-
Horizontal position
xof a particle in given global Cartesian coordinates [m]. - PX
-
Horizontal canonical momentum [eV/c].
- Y
-
Longitudinal position
yof a particle in global Cartesian coordinates [m]. - PY
-
Longitudinal canonical momentum [eV/c].
- Z
-
Vertical position
zof 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 OPAL-t and OPAL-cycl
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)^2-1}.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 OPAL-cycl 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 OPAL-cycl, 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. OPAL-cycl uses Lagrange’s 5-point 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
field-read 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.414e-03 3.743e-03 8.517e-03 1.221e-02 2.296e-02
3.884e-02 5.999e-02 8.580e-02 1.150e-01 1.461e-01
1.779e-01 2.090e-01 2.392e-01 2.682e-01 2.964e-01
3.245e-01 3.534e-01 3.843e-01 4.184e-01 4.573e-01
......
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 post-processing 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 field-map
It is additionally possible to load 3D field-maps for tracking through
OPAL-cycl. 3D field-maps are loaded by sequentially adding new field
elements to a line, as is done in OPAL-t. 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="fdf-tosca-field-map.table", LENGTH_UNITS=10., FIELD_UNITS=-1e-4;
l1: Line = (ringdef, triplet, triplet);
The field-map with file name fdf-tosca-field-map.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.68275932346E-07 -5.3752492577 0.28280706805E-07 194.36351 0.0000000 79.516210 0.42525693524E-07 -5.3827955117 0.17681348191E-07 194.70861 0.0000000 78.667380 0.19766168358E-07 -5.4350026348 0.82540823165E-08 <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 field-map
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 two-gap cavity is treated as two
independent single-gap cavities. the spiral gap cavity is not
implemented yet. For more detail about the parameters of cyclotron
cavities, see Section [cavity-cycl].
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 field-map
The 3D RF field-map 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. OPAL-cycl uses a
fourth order Runge-Kutta algorithm and the second order Leap-Frog
scheme. The fourth order Runge-Kutta 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 sub-steps: first, the code tracks this
particle for t_1 = \tau - (t_c-t_{i-1}); 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
OPAL-cycl uses the same solvers as OPAL-t 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 Boris-Buneman time integrator (leapfrog-like scheme) which has per default only one force evaluation per step. The fourth order Runge-Kutta integrator keeps the space charge field constant for one step, although there are four external field evaluations. There is an experimental multiple-time-stepping (MTS) variant of the Boris-Buneman/leapfrog-method, which evaluates space charge only every N^th step, thus greatly reducing computation time while usually being still accurate enough.
1.8. Multi-bunch 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 OPAL-cycl, we developed a new fully consistent algorithm of
multi-bunch 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, multi-bunch
simulation starts from the injection points.
In the space charge calculation for multi-bunch, 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 OPAL-cycl 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 OPAL-cycl 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.
Multi-particle 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 file-format [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 off-centering 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="OPAL-cycl: 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.0e-06, EY = 1.0e-06, ET = 1.0e-06;
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.2E-3, CHARGE=1.0,
BFREQ= frequency;
Select, Line=l1;
TRACK,LINE=l1, BEAM=beam1, MAXSTEPS=nstep*turns, STEPSPERTURN=nstep,
TIMEINTEGRATOR="RK-4";
run, method = "CYCLOTRON-T", 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= 1e-06* EY= 1e-06* ET= 1e-06* 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> * Sigma-Matrix 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.344e-18 [C] OPAL> * ********************************************************************************** OPAL> * Selected Tracking Method == CYCLOTRON-T, NEW TRACK OPAL> * ********************************************************************************** OPAL> * Number of neighbour bunches= 1 OPAL> * DT = 1e-12 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.129e-05 [Beta Gamma] OPAL> * SIGMAPY = 3.639e-05 [Beta Gamma] OPAL> * SIGMAPZ = 4.49e-05 [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> * N-PROCESSORS 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.34353e-02 [nC] Qi = 4.34353e-09 [nC] OPAL> * Ekin = 5.80000e+02 [MeV] dEkin = 0.00000e+00 [MeV] OPAL> * rmax = ( 6.17284e-03 , 7.79039e-03 , 4.79664e-03 ) [m] OPAL> * rmin = ( -6.17191e-03 , -7.78887e-03 , -4.79646e-03 ) [m] OPAL> * rms beam size = ( 1.91510e-03 , 2.41332e-03 , 1.50548e-03 ) [m] OPAL> * rms momenta = ( 3.99780e-05 , 3.51484e-05 , 4.37486e-05 ) [beta gamma] OPAL> * mean position = ( -1.56300e-13 , 2.13395e-17 , -8.71013e-21 ) [m] OPAL> * mean momenta = ( 5.72459e-17 , 1.27218e+00 , 2.15345e-17 ) [beta gamma] OPAL> * rms emittance = ( 5.52530e-08 , 4.34281e-08 , 5.13673e-08 ) (not normalized) OPAL> * rms correlation = ( -3.96345e-01 , 7.58800e-01 , 1.24755e-01 ) OPAL> * hr = ( 3.98218e-04 , 5.02557e-04 , 6.39540e-04 ) [m] OPAL> * dh = 2.00000e+00 [%] OPAL> * t = 0.00000e+00 [s] dT = inf [s] OPAL> * spos = -8.71013e-21 [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.4758e-05 Deg, Z = -5.4519e-05 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.34353e-02 [nC] Qi = 4.34353e-09 [nC] OPAL> * Ekin = 5.80000e+02 [MeV] dEkin = 3.22716e-02 [MeV] OPAL> * rmax = ( 6.17284e-03 , 7.79039e-03 , 4.79664e-03 ) [m] OPAL> * rmin = ( -6.17191e-03 , -7.78887e-03 , -4.79646e-03 ) [m] OPAL> * rms beam size = ( 1.91510e-03 , 2.41332e-03 , 1.50548e-03 ) [m] OPAL> * rms momenta = ( 3.99780e-05 , 3.51500e-05 , 4.37486e-05 ) [beta gamma] OPAL> * mean position = ( 2.13000e+00 , 9.20405e-07 , -5.45188e-08 ) [m] OPAL> * mean momenta = ( -1.73993e-04 , 1.27218e+00 , -6.34947e-05 ) [beta gamma] OPAL> * rms emittance = ( 5.52530e-08 , 4.34328e-08 , 5.13673e-08 ) (not normalized) OPAL> * rms correlation = ( -3.96344e-01 , 7.58765e-01 , 1.24755e-01 ) OPAL> * hr = ( 3.98218e-04 , 5.02557e-04 , 6.39540e-04 ) [m] OPAL> * dh = 2.00000e+00 [%] OPAL> * t = 0.00000e+00 [s] dT = inf [s] OPAL> * spos = -5.45188e-08 [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.34353e-02 [nC] Qi = 4.34353e-09 [nC] OPAL> * Ekin = 5.80000e+02 [MeV] dEkin = 3.22716e-02 [MeV] OPAL> * rmax = ( 6.17284e-03 , 7.79039e-03 , 4.79664e-03 ) [m] OPAL> * rmin = ( -6.17191e-03 , -7.78887e-03 , -4.79646e-03 ) [m] OPAL> * rms beam size = ( 1.91510e-03 , 2.41332e-03 , 1.50548e-03 ) [m] OPAL> * rms momenta = ( 3.99780e-05 , 3.51500e-05 , 4.37486e-05 ) [beta gamma] OPAL> * mean position = ( 2.13000e+00 , 9.20405e-07 , -5.45188e-08 ) [m] OPAL> * mean momenta = ( -1.73993e-04 , 1.27218e+00 , -6.34947e-05 ) [beta gamma] OPAL> * rms emittance = ( 5.52530e-08 , 4.34328e-08 , 5.13673e-08 ) (not normalized) OPAL> * rms correlation = ( -3.96344e-01 , 7.58765e-01 , 1.24755e-01 ) OPAL> * hr = ( 3.98218e-04 , 5.02557e-04 , 6.39540e-04 ) [m] OPAL> * dh = 2.00000e+00 [%] OPAL> * t = 0.00000e+00 [s] dT = inf [s] OPAL> * spos = -5.45188e-08 [m] OPAL> * ********************************************************************************** OPAL> OPAL> * ----------------------------------------------- * OPAL> * Finalizing i.e. write data and close files : OPAL> * Finalize cyclotron OPAL> * ----------------------------------------------- *