Skip to content

GitLab

  • Projects
  • Groups
  • Snippets
  • Help
    • Loading...
  • Help
    • Help
    • Support
    • Community forum
    • Submit feedback
  • Sign in
O OPALManualWiki
  • Project overview
    • Project overview
    • Details
    • Activity
    • Releases
  • Repository
    • Repository
    • Files
    • Commits
    • Branches
    • Tags
    • Contributors
    • Graph
    • Compare
    • Locked Files
  • Issues 0
    • Issues 0
    • List
    • Boards
    • Labels
    • Service Desk
    • Milestones
    • Iterations
  • Merge requests 0
    • Merge requests 0
  • CI/CD
    • CI/CD
    • Pipelines
    • Jobs
    • Schedules
  • Operations
    • Operations
    • Incidents
    • Environments
  • Analytics
    • Analytics
    • CI/CD
    • Code Review
    • Issue
    • Repository
    • Value Stream
  • Wiki
    • Wiki
  • Snippets
    • Snippets
  • Members
    • Members
  • Activity
  • Graph
  • Create a new issue
  • Jobs
  • Commits
  • Issue Boards
Collapse sidebar
  • snuverink_j
  • OPALManualWiki
  • Wiki
  • opalcycl

Last edited by Jochem Snuverink Sep 11, 2017
Page history

opalcycl

Table of Contents
  • 1. OPAL-cycl
    • 1.1. Introduction
    • 1.2. Tracking modes
    • 1.3. Variables in OPAL-cycl
    • 1.4. Field Maps
    • 1.5. RF field
    • 1.6. Particle Tracking and Acceleration
    • 1.7. Space Charge
    • 1.8. Multi-bunch Mode
    • 1.9. Input
    • 1.10. Output
    • 1.11. Matched Distribution

1. 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 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 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].

Figure 1: 2D field map on the median plane with primary direction corresponding to the azimuthal direction, secondary direction to the radial direction

CarbonFieldFormat

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> * ----------------------------------------------- *
Clone repository
  • autophase
  • beam command
  • benchmarks
  • control
  • conventions
  • distribution
  • elements
  • fieldmaps
  • fieldsolvers
  • format
  • geometry
  • Home
  • introduction
  • lines
  • opal madx
View All Pages