Computer codes for beam dynamics analysis of cyclotronlike accelerators

Computer codes suitable for the study of beam dynamics in cyclotronlike (classical and isochronous cyclotrons, synchrocyclotrons, and fixed field alternating gradient) accelerators are reviewed. Computer modeling of cyclotron segments, such as the central zone, acceleration region, and extraction system is considered. The author does not claim to give a full and detailed description of the methods and algorithms used in the codes. Special attention is paid to the codes already proven and confirmed at the existing accelerating facilities. The description of the programs prepared in the worldwide known accelerator centers is provided. The basic features of the programs available to users and limitations of their applicability are described.


I. INTRODUCTION
Cyclotronlike [classical and isochronous cyclotrons, synchrocyclotrons, and fixed field alternating gradient (FFAG)] accelerators are broadly being used, e.g. in fundamental physics research as well as in applied areas such as medicine.With a steadily growing number of cyclotronlike accelerators in the design stage, the need for thorough accelerator modeling continues.The basic operational principles of cyclotrons such as time-constant magnetic field, spiral ion trajectory, injection and extraction systems lead to specific requirements on the computer codes intended for the beam dynamics analysis.In the 1950s and 1960s, analytical methods and programs in machine-language codes were commonly used.In the 1970s and 1980s, the first special computer codes for cyclotronlike accelerators were developed.Today there are many computer codes for the design of cyclotronlike accelerators with some developed specifically for cyclotronlike accelerators, while others evolved from more general codes.Some of the codes deal with a specific segment such as injection line, central region, accelerating zone, or extraction system, while others are appropriate for end-to-end cyclotronlike accelerator design.
A brief overview of the state-of-the-art codes for beam dynamics analysis of cyclotronlike accelerators is presented.The description of the main parameters of the codes and their scopes is given.

II. CYCLOPS
The CYCLOPS code [1] is perhaps the best-known tool for beam dynamics analysis in cyclotrons.It was initially developed at Oak Ridge National Laboratory [2] and Michigan State University (MSU) [3] and further elaborated at the Tri-University Meson Facility (TRIUMF).The program calculates all the important properties of closed orbits and the linear radial and vertical oscillations about these orbits.A companion code GOBLIN [4] calculates the accelerated orbit.
The CYCLOPS input data is a 2D magnetic field map with field components in the polar coordinate system in the cyclotron median plane providing the ability to track particles through measured magnetic fields.Numerical integration of the canonical equations of motion is provided.An independent variable is the azimuth.The code finds transfer matrices and tunes by integrating first-order ordinary differential equations.
The program provides at each energy, the tunes, eigenellipse parameters, and form factors for the linear oscillations (maximum value of modulation of betatron oscillation amplitude [5]), as well as data on the frequency error and phase slip.
TRIUMF has a special version of CYCLOPS that accommodates situations where the equilibrium orbit is slightly displaced from the nominal cyclotron median plane as a result of, for example, magnet imperfections.In addition to the usual map of axial component of magnetic induction B z , this version of the code requires values of radial and azimuthal field components B r , B θ and ∂Bz=∂z at z ¼ 0. It then uses modified equations of motion, which include zero-and first-order terms in z and p z to calculate the displaced equilibrium orbit and the resultant focusing properties.CYCLOPS is TRIUMF property [6], though it has been distributed to many labs.
CYCLOPS is not appropriate to perform precise ion tracking in the cyclotron central region or to match the particle trajectory with the elements of the extraction system.However, it does support the design of a magnetic system with stable rf phase slip and effective ion focusing, and is utilized worldwide to determine the necessary magnetic structure especially where high magnetic field accuracy is required, e.g. the TRIUMF cyclotron.As an illustration, the tune diagram of the TRIUMF [7] cyclotron shown in Fig. 1 was published in 1975, but CYCLOPS remains a useful and fast code.In addition, the program can calculate unstable fixed point orbits associated with nonlinear resonances, which allows studying resonance crossing and modifying the magnetic field to avoid dangerous crossings.For example, this analysis has allowed TRIUMF to successfully operate at the 2Q r ¼ 3 resonance.
CYCLOPS can also be used for simulation of FFAG accelerators [8].In this case, it is necessary to use a magnetic field map that has a very fine grid to get good equilibrium orbit solutions for the magnet with hard edges.An example is the Johnstone-Koscielniak medical FFAG for cancer therapy for 18-400 MeV=u carbon ions [9].In this case, the CYCLOPS results were satisfactory and the fields were obtained by smoothing the hard edges by fitting a Fourier series combined with the Lanczos σ-factors.

III. COSY INFINITY
COSY INFINITY [10] uses advanced concepts of modern scientific computing such as graphics libraries and parallel computation.The code is broadly applicable including cyclotrons, synchrotrons, FFAGs, electrostatic accelerators, linacs, RFQ accelerators, spectrometers, and transport lines allowing injection-to-extraction simulations.COSY currently has more than 2500 registered users and has been extensively cross-checked and verified.The code is written in F77, F90, C þ þ, and runs on PC and UNIX platforms.The program is free, and can be downloaded from an MSU site [11].
There are various methods to perform tracking in COSY that preserve the simplistic symmetry inherent in Hamiltonian systems with minimal modifications allowing a reliable estimation of the dynamic aperture.The program performs the computation of all beam dynamics of the system to an arbitrary order, including out-of-plane expansions of fields and any nonlinear terms in the Hamiltonian.For efficient initial simulation and optimization, it is particularly useful to utilize very high order outof-plane expansions of suitable midplane models.There are various tools for analysis of nonlinear effects, including normal form-based computation of high-order amplitudedependent tune shifts and resonances.An algorithm for the evaluation of closed equilibrium orbits allows calculating betatron frequencies.
There are various algorithms to treat the fields of structural elements efficiently: the midplane field for a midplane symmetric element; the on-axis potential for straight elements like solenoids, quadrupoles, and higher multipoles; field representations from 3D data; measured fields.There is a possibility of importing TOSCA/OPERA [12] field maps.
The code implements rf and electrostatic structures intended for particle acceleration, bunching, and focusing.Field models for high accuracy evaluation of the Taylor maps of cylindrically symmetric electrostatic lenses (zero rf frequency) and time-varying cavities are used [13].Several procedures have been developed for COSY that make it easier to accurately model nonidealized electrostatic elements and their fringe fields [14].These procedures also allow an electrostatic element to be defined from measured or simulated potential data.TOSCA/OPERA can be used to provide the electrostatic potential map.
A method is implemented in COSY that allows the computation of space charge effects of arbitrary and large distributions of particles in an efficient and accurate way based on a variant of the fast multipole method (FMM) [15].This algorithm divides an arbitrary charge distribution into small boxes with a hierarchical structure.It then computes the multipole expansions and local expansions of charges far from the observer to achieve a computation efficiency that scales with the number of particles and computational errors that scales with a high power of the differential algebra order.The FMM algorithm is especially suited for beam dynamics simulations in accelerators with a large number of turns because of the efficiency and low computational error compared to other space charge algorithms.
Significant enhancements of the code COSY INFINITY for the particularly challenging case of FFAG accelerators have been implemented.There are several examples of applying COSY to the design of FFAG accelerators [16,17].A midplane field map can be generated using COSY from a special complex supported in the program.The code allows accurate studying of dynamic aperture in the accelerators, in which ions make a large number of turns during the acceleration process.The user can see effects of using various orders of out-of-plane expansion and the algorithm of symplectification (see Fig. 2).
COSY INFINITY can be used for compact and separated sector cyclotron design.The code was used to study the FIG. 1. Square of an axial frequency as a function of radius at TRIUMF cyclotron: 1-calculation with CYCLOPS (solid line), 2-measurements (data points) [7].
betatron frequencies of a 250 MeV superconducting proton cyclotron [18].Recently a graphical user interface interface developed by the Particle Accelerator Corporation (PAC) and COSY INFINITY and called FACT (FFAG and cyclotron tools) provides a graphical front end to the COSY.It fully develops and expands the complex field and edge profiles to accurately predict and optimize machine performance.An exact, 3D field expansion in polar coordinates is one of the output formats which can be used by other codes.

IV. OPAL
OPAL is an open-source tool for charged-particle optics in accelerator structures and beam lines [19].Using the MAD language with extensions, OPAL is derived from MAD9 and provides accurate tracking and comprehensive benchmarking against existing high-intensity accelerators.The framework provides parallel calculations, which become possible through the increasingly sophisticated mathematical models and evolving computer power available on the desktop as well as supercomputer centers.OPAL can run on a user's laptop as well as on the largest clusters.It may be the most effective tool for the study of space charge effects in particle accelerators.This goal can only be achieved with detailed 3D modeling capabilities and a sufficient number of simulation particles to obtain meaningful statistics on various quantities of the particle ensemble such as emittance, slice emittance, halo extension, etc.
The program has a special flavor OPAL-CYCL, which is also fully three dimensional, dedicated to future high intensity cyclotron and FFAG accelerator design.It tracks multiparticles, taking into account the space charge effects.OPAL-CYCL uses the fourth-order Runge-Kutta algorithm, second-order leapfrog scheme, and a time adaptive scheme [20].In addition to multiparticle simulation mode, OPAL-CYCL also has two other serial tracking modes for the 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 the reference orbit, phase shift history, stable region, and calculation of the eigenellipse.The other mode 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 plug-in elements, including the 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 cyclotrons and FFAGs.
In the code the magnetic field is read from an ASCII type file.The field data should be stored in the cylindrical coordinates frame (because the field map on the median plane of the cyclotron is usually measured in this frame).It is the real field map in the median plane of the existing cyclotron machine obtained using measurement equipment.Alternatively the field can be calculated numerically in the median plane or the 3D fields of the working gap for the region of interest by creating a 3D model using commercial software, such as TOSCA\OPERA [12] and ANSYS [21] during the design phase of a new machine.If the field on the median plane is calculated, the off median plane field can be obtained using the same expansion approach as the measured field map.
In the tune calculation mode, two particles are tracked; one, with all data set to zero, is the reference particle and the other is off-centered in both r and z directions.Working in this mode, the program can calculate the betatron oscillation frequencies Q r and Q z for different energies to evaluate the focusing characteristics for a given magnetic field.Theoretically, there is an eigenellipse for any given energy in the stable area of a fixed accelerator structure.When the initial phase space shape matches the eigenellipse, the beam envelop amplitude oscillation is minimal and the transmission efficiency maximal.The eigenellipse can be calculated by a single particle tracking of the off-centered particle, and its coordinates and momenta can be recorded at the same azimuthal position at each turn.In OPAL-CYCL the Fourier analysis method is used to plot the orbit difference between the reference particle and an initially displaced particle.The frequency point with the largest amplitude is the betatron tune value at the given energy.
In the current version of OPAL-CYCL, rf cavities are treated as straight lines with infinitely narrow gaps, and the electric field is a δ function plus a transit time correction.A two-gap cavity is treated as two independent single-gap cavities.The voltage profile of a cavity gap is read from the ASCII file.
OPAL-CYCL provides capability to include the central region of a compact cyclotron [22] where the user can generalize the reference coordinate system from 2D to 3D.The intermediate phase space data is stored in an ASCII file which can be used to plot the orbit.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 [23] and can be analyzed using the H5Root [24].
OPAL-CYCL was used to do numerical simulation of the PSI injector II and the PSI ring cyclotron.Nonlinear space charge effects, including mutual interaction of neighboring bunches, were studied [25].For the PSI ring cyclotron the total sum of controlled and uncontrolled losses is about 10 −4 .Beam dynamics simulations with this precision must include the effects of radial neighboring turns.It was found that neighboring bunch effects caused the beam bunch in the PSI ring to become more compact in the transverse direction with a slightly reduced energy spread contributing to a reduction in beam loss in high intensity operation (Fig. 3).
OPAL was also used for modeling the compact cyclotron CYCIAE-100 (CIAE, China) [26] and the DAEδALUS cyclotrons [27].Lattice design and tracking for the 1 GeV FFAG was also performed with this code [28].

V. ZGOUBI
The ZGOUBI computer code [29] calculates trajectories of charged particles in magnetic and electric fields.It uses a sixth-order Taylor series integration method with fields from analytic magnet models or from simulated or measured field maps.The initial version of the code, dedicated to ray tracing in magnetic spectrometers, was developed by Garreta and Faivre at CEN-Saclay in the early 1970s.It was perfected for the purpose of studying the four spectrometers SPES I, II, III, IV at the Laboratoire National Saturne (CEA-Saclay, France), and SPEG at Ganil (Caen, France), and large spectrometers at the PS experimental areas (CERN) and continues to be used in laboratories worldwide.Present developments of ZGOUBI have strongly benefited from the environment of the Groupe Théorie, Laboratoire National SATURNE, CEA/DSM-Saclay, in the period 1986-1998.The code has had ∼3500 downloads from ∼60 countries over the past ten years [30] and is available from the source Forge development site [31].
ZGOUBI is accompanied by documentation and a detailed user's manual.Practical examples are provided including input/output data files.It has evolved to allow the study of systems including complex sequences of optical elements such as dipoles, quadrupoles, arbitrary multipoles, FFAG magnets, and other magnetic or electric devices, and is also able to handle periodic structures.The code has several unique attributes including: (i) computation of optical functions, including periodic optics, from stepwise raytracing; (ii) multiturn tracking in circular accelerators, using field maps as well as analytical models of magnetic and electrostatic fields; (iii) spin tracking using the same numerical method as for the Lorentz equation; (iv) inclusion of stochastic photon emission and its effects on particle dynamics; (v) calculation of the synchrotron radiation electric field and spectra in arbitrary magnetic fields, using ray tracing; (vi) a possibility of using a mesh, which allows ray tracing from simulated or measured (1D, 2D or 3D) electric and magnetic field maps; (vii) a number of Monte Carlo procedures, including an unlimited number of trajectories, in-flight decay, stochastic radiation, etc.; (viii) a built-in fitting procedures allowing arbitrary variables and a large variety of constraints that are easily expandable; and (ix) simulation of time-varying power supplies.
The ZGOUBI numerical method for integrating the Lorentz equation is based on the Taylor series, which optimizes computing time and provides high accuracy.The integration variable is the arclength along the particle trajectory.
Magnetic field maps from measurements and calculated with TOSCA\OPERA and POISSON SUPERFISH can be based either on Cartesian or polar coordinates.Magnetic induction extrapolation off an arbitrary 2D plane, without any hypothesis of symmetries, is also available (e.g., it is used to handle Siberian snake field maps at RHIC).ZGOUBI provides three types of polynomial interpolation from the mesh, namely, a second-order interpolation, with either a 9-or a 25-point grid, or a fourth-order interpolation with a 25-point grid.If the field distribution result from measurements, the grid also smoothes possible field fluctuations.In 3D field maps magnetic induction and its derivatives up to the second order with respect to axes of Cartesian coordinate system are calculated by means of a secondorder polynomial interpolation.
The code calculates electric field and its derivatives in various ways.Several optical elements, defined analytically, have their field defined from the expression of the field and derivatives in the median plane.3D extrapolation of these off the median plane is drawn from Taylor expansions.Field maps treated in the actual version include the 1D axial maps with cylindrical symmetry and 2D field maps.The program supports accelerating elements such as a synchrotron cavity that can include synchrotron motion.A cavity with a fixed rf frequency can be used for cyclotron design.The user can specify a number of such elements that are used in the simulation.Voltage dependence on time, synchronous phase on time, or other parameters of the accelerating system can be imported from a separate external file.
Very effectively ZGOUBI can be used for calculation of particle trajectories, betatron tunes, beam envelopes, and space charge forces [32,33].Transverse space charge models were added to ZGOUBI using the interface PYZGOUBI.PYZGOUBI has been developed in Python, which allows scripting ZGOUBI and makes it easy to interface to other programs or functions.PYZGOUBI splits accelerator lattice elements into slices.Particles are tracked through a slice using ZGOUBI, and then their coordinates are returned to PYZGOUBI.Then a space charge kick is applied to the particles before tracking them through the next slice.Communication between PYZGOUBI and ZGOUBI is done via binary files; these can be in shared memory for increased performance.Multiple instances of the code can be run in parallel to reduce computational time for tracking a large number of particles.Three algorithms of calculations of space charge can be implemented [33]: a linear method, equivalent to the Kapchinskij-Vladimirskij (KV) envelope equations on a beam of uniform charge density; an N-body method sums up the Coulomb force between each pair of particles; a method that decomposes the beam into concentric elliptic rings, which is suitable for elliptically symmetric beams.
For example, ZGOUBI was used to study space charge dominated cases for the design of an accelerator-driven subcritical reactor nonscaling FFAG (Fig. 4).It is the strongest at low energies, whereas without space charge the tune would fall with increasing energy.While it is not possible to use this to cancel out tune shift, it does reduce the total shift during acceleration.
The code has tools for definition of particle losses on an accelerator structure by specifying boundary where the ion could be lost.
A series of auxiliary computing tools have been developed enabling the search for closed orbits in periodic machines, computation of optical functions and parameters, tune scans, dynamic aperture scans, spin dynamics data treatment, and graphic scripts based on GNUPLOT.An interactive graphic/analysis postprocessor that produces plots of ZGOUBI output data, Fourier analysis, ellipse matching, graphic interface to ZGOUBI (ZPOP) has undergone extensive developments.
The code is a widely used code for modeling FFAGs.PYZGOUBI adds a flexible layer to it, which can be used to add new features without making intrusive changes in the existing code.ZGOUBI was used as the on-line simulation code for the EMMA project after its use during the design period.It has been used in the analysis of data taken over the few EMMA run years [34].
Step-by-step particle tracking was done using TOSCA \OPERA3D field maps and the ZGOUBI code to study beam dynamics in the RACCAM FFAG [35].The code was used for tune analysis in the KEK 150 MeV FFAG [36] and for the study of the injection beam line and evaluation of the equilibrium orbit in the DAEδALUS cyclotron [37].It is at present extensively used for spin dynamics and other Siberian snake manipulations at the RHIC polarized proton complex and in the eRHIC EIC project.

VI. CYCLONE, Z3CYCLONE
The CYCLONE code was developed in the mid-1960s at MSU as part of the K50 cyclotron project [38].The program has since been used to design central regions for the MSU K500, K1200, and K100 cyclotrons, and the Milan K800 cyclotron.In 1987 the program was brought to TRIUMF from MSU and was used in the design of the central region of the TR30 cyclotron.At TRIUMF the program has been modified [39] to make it more flexible and was used for simulation for TRIUMF cyclotron.The code is now also available at National Superconducting Cyclotron Laboratory MSU and is named Z3CYCLONE [40].CYCLONE and Z3CYCLONE are supported by TRIUMF and MSU accordingly and have similar characteristics.
CYCLONE is an orbit tracking code that can be used to study all regions of a cyclotron, but it is particularly suited for central region design work.Orbits are calculated using Runga-Kutta type integration in a variety of magnetic and electric field maps.The equations of motion used in CYCLONE are fully relativistic.The program is divided into three parts, each with a different integration approach.The first part integrates equations of motion in the given electric and magnetic field maps using time, as an independent variable.This is useful near the center of the cyclotron where orbits may deviate a lot from the circle.It is however slow because it needs to look up the magnetic field values at In part III the vertical motion is linear.In parts I and II the magnetic contribution is linear but the electric contribution can be either linear or nonlinear depending on the electric field mesh used.It is possible to use separate electric field maps in parts I and II.This allows having a more detailed mesh for calculation of electromagnetic fields in part I where the particle energy is still quite low.Provision is made to position these fields at arbitrary locations relative to the magnetic field.The field maps are stored in a format that is compatible with RELAX3D, so direct output from that code can be used as input.Both electric fields are stored in rectangular meshes.In part II the mesh may also contain multiple z-layers, which allow more accurate calculation of the vertical electric field.The magnetic field data is stored in a standard polar mesh.The CYCLONE code can handle three phase rf systems with spiral dees, operating in any harmonic mode.
Equilibrium orbits are generated internally and data can be stored in a file and read back in on future runs.There is a variety of configurations for the orbits, so x and p x can be calculated in any gap, at any angle, or along a fixed angle.Unfortunately, the program does not allow space charge calculations.CYCLONE uses a command driven input structure.The program manipulates data contained in many different files, most of which are optional.
The code is very popular in scientific organizations around the world and is very useful, especially for the design of the cyclotron central region.In addition to the MSU and TRIUMF cyclotrons design CYCLONE was used for modeling a lot of cyclotrons with internal ion source [41] and machines with axial injection [42].The code can be used not only for center structure optimization and betatron frequencies calculation but also for beam centering analysis (Fig. 5), as was done, for example, for the KIRAMS-430 superconducting cyclotron [43].

VII. SNOP
The SNOP program [44] was developed at JINR (Dubna, Russia) to analyze beam dynamics in compact cyclotrons and is an extension of the cyclotron beam dynamics analysis code (CBDA) [45].SNOP provides full 3D particle tracing in cyclotrons from injection to extraction.The main features of SNOP are usage of 3D electric and magnetic field maps, 3D beam space charge calculation, and analysis of the beam losses on the structure elements of the facility.It is possible to simulate realistically the axial injection line, central zone, acceleration region, and extraction.SNOP calculates trajectories of macroparticles in electromagnetic fields, solving a complete system of equations of motion without any simplification.It uses the classical fourth-order Runge-Kutta method to solve the system of the equations with time as an independent variable.For space charge calculation the user can choose either the direct particle-toparticle or the particle-in-cell method.
Optical elements are available to construct a computer model of a cyclotron, including magnetic quadrupoles, solenoids, multiharmonic buncher, spiral electrostatic inflector, harmonic coils, electrostatic deflectors, magnetic channels, and stripping foils.The main magnetic field can be uploaded as both a 2D median-plane field in polar coordinates and a 3D map.In the first case, a full magnetic field can be restored.Practically all field formats are correlated with the TOSCA\OPERA3D program, which allows electromagnetic field calculations.
Since the 3D electric field can be used, any accelerating cavity type is supported including spiral dees.The program uses the static 3D field map of the dees and allows changing peak dee voltage separately for each dee.It is possible to introduce dependence of the dee voltage on the radius.The Cartesian coordinate system is used in the program as the most generally applicable to determining electromagnetic fields of the elements with a complicated structure, such as the inflector, extraction system elements, or harmonic coils.In addition, SNOP allows using analytical definition of the fields of some elements, for instance, in the extraction system.This feature is very suitable at the initial stage of the computer modeling.A distinctive feature of SNOP is the possibility to use 3D electric and magnetic field maps that are asymmetrical relative to the cyclotron median plane.The user can also use a built-in command to introduce the first harmonic of the magnetic field with varying parameters.Importantly, this allows evaluation of the effects of various system misalignments in support of developing misalignment tolerances.
SNOP is intended for the Windows platform and distributed freely by the authors as an executable file with the corresponding user's manual.The code is convenient to use due to its user-friendly interface (Fig. 6).SNOP is structured in such a way that there are dedicated blocks responsible for facility structure units in the beam transport line and accelerator itself.There are also separate menus to control parameters of particle acceleration and to assign parameter values for the equation of particle motion including beam space charge calculation.
The program permits controlling parameters of electromagnetic devices for beam acceleration and focusing in so that it is possible to modify available parameters without editing any files manually.SNOP allows tracing both a single trajectory and a 3D particle bunch.Also, the program allows calculating amplitudes of radial betatron oscillations.SNOP has no internal tool for direct calculation of betatron tunes, but the user can find the equilibrium orbit and calculate frequencies using its classical definition.In this way, two particles are tracked, one moves along the equilibrium orbit, and the other is off-centered in either the r or the z direction.Betatron tune is a number of particle oscillations around its equilibrium orbit divided by the number of turns that the ion makes in the magnetic field.
In SNOP there is a possibility of detecting loss of particles when they cross the surfaces of the accelerator structure elements.One of the options is description of the structure geometry by analytical surfaces.This permits fast estimation of the particle losses in the facility regions where this geometry input is valid.For more accurate calculation of losses the geometry can be imported from the CAD program with sufficiently detailed presentation of the mechanical model.It can be used in those parts of the facility where the geometry is more complex and cannot be described by analytical surfaces.Fast algorithms for determining the particle intersection with the surfaces of the bodies are used to decrease the required computer time for the calculations in SNOP.
Full information about motion of ions is contained in the ASCII files.The general output file contains coordinates, velocities, energy, rf phase of particles, etc. SNOP supports using such systems as MATHCAD and AUTOCAD, with which data exchange can be carried out.MATHCAD is used at the initial data generation and for analysis of the results.AUTOCAD permits one to specify geometry of the objects when calculating particle losses and to depict positions of the accelerated and lost particles against the background of the real facility structure (Fig. 7).
Since SNOP has general compact cyclotron blocks, input data can be presented in the form of 3D field maps, and setup geometry is imported from the CAD program, the code can be used for the design of any cyclotron.It can be applied not only to compact machines but also to separated sector cyclotrons and FFAG accelerators.SNOP is installed at several laboratories in the world.It was successfully used for designing of the unique ultracompact superconducting cyclotron ION-12SC [46], modeling the NIRS-930 cyclotron [47], and the separated sector cyclotron for carbon therapy K1600 [48], and for some other cyclotronlike accelerators.Future plans include extensions of the code for large ring cyclotrons and synchrocyclotrons.For example, an algorithm of using time-dependent rf frequency is under development.

VIII. TRACK
The computer code named TRACK was developed at Paul Scherrer Institut (PSI) (Zurich, Switzerland) as a general purpose program for particle tracking in electric and magnetic fields [49].The program was originally conceived to facilitate the design and analysis of beam line magnets, spectrometers, and separators.However, the built-in command language interpreter allows a broader use of the program for beam dynamics analysis in cyclotronlike accelerators.The code uses the method of tiny arcs for particle tracking.The program follows particle trajectories in three dimensions using an algorithm that contains the analytical solution of the equation of motion.The integration is done in small constant field steps of varying size, optimized for the local field variation.The input fields can be obtained from commercially available codes, mathematical models, or field measurements.Particle histories can be calculated as a function of position, momentum, mass and time in combined static and alternating magnetic and electric fields.The program provides 3D visualization of the particle tracks and evolution of the phase space, and marks the particle positions as a function of the phase when a time-dependent field is used.The program was successfully used to design the compact superconducting cyclotron ACCEL/VARIAN [50].Since the program was not written for a specific application, its command language interpreter has to be used to extract relevant data from the particle tracks at specific locations, moments in time or at other conditions specified by the user.For the ACCEL cyclotron precessional beam extraction was analyzed with TRACK (Fig. 8).

IX. PHASCOL
The PHASCOL code [51] was developed at Joint Institute for Nuclear Research (JINR) (Dubna, Russia) for particle dynamics computation allowing for space charge effects in cyclotrons.The Runge-Kutta method of the fourth order is used for numerical integration of equations of motion.All computations are done in the double precision.For numerical modeling of space charge effects, the beam is subdivided into a series of bunches.Each bunch contains a set of macroparticles.PHASCOL provides a possibility of adding a certain number of injected bunches following one another in the assigned time intervals.Two methods are used for calculating the self field of the beam: the method of pairwise interactions and the particle-in-cell method.In the first method, the electric field in the location of the macroparticle is calculated by summing up the Coulomb field, which acts from the side of the remaining particles.In the second method, with the aid of fast Fourier transform, the 3D Poisson equation on the grid that covers beam is solved numerically.
The magnetic field is entered using the radial dependencies of Fourier harmonics in the median plane of the accelerator.Outside the median plane the magnetic field components are computed in accordance with the Maxwell equations.
For visual representation of the input data and the results of computation, the graphical utilities of the Microsoft Fortran Power Station (FPS) 4.0 are applied.These graphs show the starting and final positions of the particles in the transverse and longitudinal phase planes, and the radial and axial particle trajectories in the phase planes versus either time or the azimuth angle.A signal of the differential probe located at an arbitrary azimuth and having definite radial size is also computed and plotted versus the radius.
The code was used for beam dynamics analysis in the JINR Phasotron, for design of a high-current separated sector cyclotron that was considered as an injector into the Phasotron [52], and for other accelerators.

X. MATLAB CODES
MATLAB [53] is a multiparadigm numerical computing environment and fourth-generation programming language.It can be used for accelerator modeling because numerous toolboxes of this system can provide the researcher with all necessary subprograms, including integration of equations of motion and space charge calculation.Due to its powerful graphical potential, MATLAB can be useful as a postprocessor for analysis of simulation results.So, it was logical that MATLAB was used as a platform for program complexes for cyclotron beam dynamics simulations, one of which was developed at JINR (Dubna, Russia) [54].
To accommodate any kind of acceleration setup, programs for three types of differential equation systems were written to solve systems of differential equations of particle motion written in the Cartesian coordinate system (to simulate particle motion in transport lines) and in the cylindrical coordinate system (to simulate particle motion inside cyclotron).To consider space charge effects, it is essential to use time as an independent variable.If space charge forces are not considered, it is more convenient to use the azimuth angle as an independent variable (in cylindrical coordinate system).Ion beam dynamics simulations can be made through the injection line, central region (Fig. 9), accelerating area, extraction system, and high energy beam transport.The program also allows analyzing ion beam losses due to charge exchange with the residual gas.
3D electric accelerating and magnetic field maps provided with external programs can be imported to the program.A measured magnetic field map also can be used.The program has an option of analytical estimation of the cavity voltage distribution over the radius.
The program complex was successfully used for designing the compact cyclotron CYTRACK that operates in Dubna (Russia) [55].One more application of the program is a design study of a new superconducting cyclotron for proton therapy SC202 [56] that is under development by the collaboration of JINR (Dubna, Russia) and ASIPP (Hefei, China).
Another program for designing the cyclotron central region using MATLAB was developed at INFN-LNS (Italy) [57] and was used for beam dynamics analysis in the superconducting cyclotron SCENT [58].The program produces different files describing the central trajectory and some input files for the mechanical drawing code that allows generating the mechanical model of the inflector and its collimators.This model can be imported in TOSCA \OPERA.The code allows easily printing the input file of the particle coordinates in the phase space.Many parameters of the beam distribution like the size and the position of the beam and the kind of the distribution, uniform or Gaussian, can be selected.The electric and magnetic maps are imported into a dedicated code to evaluate the real trajectory.This code allows evaluating the beam dynamics by integrating the Lorentz-Force equation using the Cartesian coordinates and time as an independent variable of integration.Development of an option for using the angular azimuth as independent variable of integration is in progress.From the analysis of the trajectory the drawing of the central region electrodes is modified and a new simulation is restarted.The code allows entering the coordinates, velocity, and phase for a number of particles.Figure 10 shows the graphical interface of the code.Its outputs are trajectories in the central region and evolution of the particle energy gain.

XI. ADVANCED ORBIT CODE
The advanced orbit code (AOC) [59] was developed by Ion Beam Applications (IBA) (Louvain-la-Neuve, Belgium) and is intended for beam dynamics studies not only for cyclotrons but also for synchrocyclotrons where the number of turns is very large.The code facilitates design studies of critical systems and processes in medical and industrial accelerators.Examples include: (i) injection into and extraction from cyclotrons, (ii) central region, beam-capture and longitudinal beam dynamics studies in synchrocyclotrons, (iii) studies of resonance crossings, (iv) stripping extraction, (v) beam simulation from the ion source to the extraction system, (vi) space charge effects, (vii) beam transmission studies in gantries, or (viii) calculation of Twiss functions.The AOC tracks accelerator orbits under the combined action of static magnetic and rf electric fields, using time as an independent variable.The equations of motion are integrated by a fifthorder Runge-Kutta method with adaptive step-size control.This results in high precision and facilitates the tracking of an accelerator from the ion source up to the beam extraction system.
For the first type, space charge effects can be included.The fully relativistic and fully 3D electric and magnetic self-fields are calculated for a bunch in free space, by accumulating for each particle the contributions of all others.For calculation of the space charge a single step or a more precise converging iterative solution can be chosen.The code uses multithreading to decrease calculation time.
There are two ways to import the static magnetic field map of a magnet into the AOC.The first method assumes that the field is specified in a 2D plane and the first or higher order field expansion is used to obtain the field outside of this plane.Such 2D maps can be obtained from measurements or from magnet calculations.For higher (third) order expansion, also the 2D Laplacian in the 2D plane needs to be given.Also median plane errors can be specified in the same plane.In the second method, the full 3D field components are specified on a 3D grid and the fields are calculated by the Lagrange polynomial interpolation.The degree of interpolation is a free parameter.The full 3D grid may consist of multiple 3D meshes each with different topology allowing a finer mesh where necessary.The same approach is used for the 3D electrical field maps.The rf electric field may be generated from different sources with different rf frequencies (including dc), rf phases, and electrode voltages.There is also a possibility of calculating electric field maps from a simplified dee geometry (using Gaussian profiles expanded around straight gaps) if 3D maps are not available.The 3D field maps can be obtained from 3D finite element packages such as TOSCA\OPERA3D.
It is very significant that the AOC code can be used for designing synchrocyclotrons with a large number of turns.
In this way, the code tracks only the relevant parameters during acceleration: the rf phase, the energy, and the orbit center coordinates.This code is referred to as the "phase motion code."In synchrocyclotrons, the rf and its voltage are modulated in time in order to ensure resonant acceleration up to final radius.There is longitudinal beam dynamics as in synchrotrons and the beam can only be captured during a short time window.This process was studied for the IBA superconducting synchrocyclotron S2C2 [60] and the result is shown in Fig. 11.

XII. OTHER CYCLOTRON CODES
This chapter deals with computer codes that can be used for designing cyclotrons, but are out of date or information about them is very scarce.
The computer code named CASINO (calculation of spiral inflector orbits) [61,62] was popular in designing cyclotron central regions and was used for modeling some accelerators [63].The program is a fourth-order Runge-Kutta tracking code.It tracks particles in a magnetic field which may be specified as homogeneous, B z ðzÞ (the field along the magnet axis), or B z ðz; rÞ (as provided by POISSON for example).The electric field may be computed in an analytic manner or computed from an electric potential map.
NAJO is a general multiparticle code for studying particle motion in cyclotrons [64].The program was developed for GANIL separated sector cyclotrons, but it could be adapted to other configurations.The code can accommodate any magnetic field configuration and its main limitation comes from the shape of the accelerating gaps which are presently restricted to radial ones.Accelerating field effects are expressed as kicks applied at the gap centers allowing for a complete decoupling in the treatment of the magnetic and electric fields.NAJO provides two methods for calculation of space charge effects: the particle-to-particle interaction method and the equivalent continuous distribution method.Both methods take into account neighboring bunches and image effects.The code was used for space charge analysis in some cyclotrons [65].
NORTICA (numerical orbit tracking in cyclotrons with analysis) is a package of computer codes for beam dynamics simulations developed at NSCL [66].The Fortran programming language was used for the numerical algorithms realization, Tcl/Tk for the user interface and C for the program interface.The codes are capable of simulating beam dynamics in both coupled cyclotron facility (CCF) accelerators, the K500 and K1200 superconducting cyclotrons.The codes can also be used for the stand-alone operation at any cyclotron.The package consists of three major parts: (i) The first part includes procedures for fitting main and trim coil currents in order to shape the main magnetic field.(ii) The next part provides the actual beam dynamics simulation in the previously generated field with imperfections and additional contributions to the total magnetic field (extraction elements and bumps).(iii) The last part involves the detailed simulation and setting of the extraction system parameters.
The user-friendly interface is provided by the program and allows the user to see clearly the dependence of the beam and field characteristics on different input parameters, which make the process of the cyclotron tuning quicker.
The General Particle Tracer (GPT) [67] is a software package developed to study 3D charged particle dynamics in electromagnetic fields.The purpose of the GPT is to aid in the design of accelerators and beam lines by using modern, particle tracking techniques.The GPT is capable of tracking any number of particles through complex electromagnetic fields taking all 3D effects and space charge forces into account.Space charge effects can be calculated using a 1D, 2D, or 3D model depending on the type of simulation.This produces highly accurate and reliable results.The GPT offers these capabilities in a very user-friendly and easy to customize package.Due to the general character of the code and its flexibility, the GPT is suited for a large number of purposes, including cyclotrons and FFAG accelerators.
The fields of physical structures are represented by the so-called elements.The GPT comes with a large set of standard elements for basic structures such as solenoids, quadrupoles, and accelerating structures.However, for specific cases custom elements can easily be added.
The equations of motion for the macroparticles are solved relativistically in the time domain using a fifth-order embedded Runge-Kutta integrator with adaptive step-size control.This allows the user to select the required accuracy, while the simulation time is always kept to a minimum.
Elements can be positioned anywhere and in any direction in 3D space.This gives the user complete freedom in constructing the experiment, and also components placed off center or rotated can be simulated.The effect of fringe fields can be taken into account and measured or externally calculated fields can be used in the simulation.The field map elements can read electromagnetic field configurations calculated by external programs like TOSCA\OPERA and the POISSON set of codes.
The GPT is accompanied by a number of pre processing and postprocessing programs and interfaces to other software packages to ease the design process.These tools are combined with the GPT and integrated into the Windowsbased graphical user interface, GPTWIN.The GPT on UNIX machines uses command-line versions of these programs.

XIII. CONCLUSION
The codes, which are created for beam dynamics analysis in cyclotronlike accelerators and available for third-party use, cover necessity for the design of a cyclotron or a FFAG accelerator including development of the central region, acceleration zone, and injection and extraction systems.Well-established methods for numerical integration of equations of motion using as independent integration variable both time and the arclength along the particle trajectory or azimuth are employed.Some of the codes utilize special symplectification algorithms that are useful for accelerators with a large number of turns.Measured and calculated magnetic and electric fields can be used.
CYCLOPS is the fastest and best-known code that provides equilibrium orbits calculation, orbital frequency analysis, and betatron frequencies estimation.However, it does not support central region design, extraction system development, and 3D beam dynamics analysis.For hard-edge magnets, preliminary preparation of the input data is necessary.
COSY INFINITY is a useful tool for studying beam optics, particularly at higher order.The code allows accurate analysis of the beam dynamics in FFAGs, cyclotrons, and synchrocyclotrons, in which ions make a large number of turns during the accelerating process.Though COSY provides high precision, better codes exist for the modeling of the cyclotron central region.
OPAL is very effective and the fastest tool for space charge effect studies in particle accelerators.The OPAL- CYCL version was created especially for cyclotrons and FFAGs and provides all necessary tools.The code can be used for the cyclotron central region design as well.
ZGOUBI is a popular tool for the design of FFAG accelerators as well as the beam dynamics simulations in cyclotrons.The code allows using many optical elements and is a flexible and reliable code, but is less user friendly for space charge calculations and employing 3D electric field maps.ZGOUBI is accompanied by documentation with practical examples.Authors of the code are very contactable and happy to provide any help.
The CYCLONE code is a useful, simple and accurate tool that is popular for designing the cyclotron central region.It provides basic possibilities for beam dynamics analysis and cyclotron structure design.But the program does not support space charge calculations and is not intended for extraction system development.
The SNOP code is a convenient tool for the beam dynamics analysis in cyclotronlike accelerators.The code was specially developed for a compact cyclotron, but can also be used for separated sector machines and FFAG accelerators.SNOP is a user-friendly code, which permits conducting calculations for a wide class of facilities apart from cyclotrons.The program supports calculations from injection line to extraction system.SNOP provides full 3D calculations with space charge effects included.However, the program does not have an internal mechanism of equilibrium orbit and betatron tunes calculation.
Codes that were developed on the basis of MATLAB offer opportunities to the user from the point of view of powerful graphical potential.MATLAB has toolboxes and the user can use the prepared tools for study.But using external framework leads to a substantial increase in the calculation time in comparison with the codes that use Fortran, C=C þ þ=C, etc.
The advanced orbit code (AOC), which is successfully used for cyclotrons, can also be applied to synchrocyclotrons with a large number of turns.The code was developed by well-known experts in the cyclotron field and tested on existing accelerators.

FIG. 3 .
FIG. 3. Top view of 1 mA bunch distributions at the turn 130 in the local frame in the PSI ring cyclotron.The results are obtained from single bunch (left) and seven bunches (right) simulations [25].

FIG. 4 .
FIG.4.Cell tune as a function of energy and peak current for 10 mm mrad beam[33].

FIG. 5 .
FIG. 5.The single particle trajectory in the first few turns with the central point of the orbit at each time step (shown in blue) in the KIRAMS-430 cyclotron [43].

FIG. 8 .
FIG. 8. Radial momentum component p r and energy as a function of radius near extraction.The data points represent all the turns at extraction and the bar at 81.6 cm represents the extraction septum [50].

FIG. 11 .
FIG. 11.Filling of the longitudinal separatrix at 3 MeV in S2C2 (top).The accepted rf phases at injection as a function of the relative injection time (bottom).The color codes indicate the link between capture time and position inside the separatrix [60].