The OPAL Framework: Version 2024.2

LaserEmission-1.in

All supplementary files can be found in the laser emission regression test.

Injector II is a separated sector cyclotron specially designed for pre-acceleration (inject: 870 keV, extract: 72 MeV) of high intensity proton beam for Ring cyclotron. It has 4 sector magnets, two double-gap acceleration cavities (represented by 2 single-gap cavities here) and two single-gap flat-top cavities.

Following is an input file of Single Particle Tracking mode for PSI Injector II cyclotron.

Injector2.in

The supplementary files should be placed in the same directory.

To run OPAL on a single node, just use this command:

Here shows some pictures using the resulting data from single particle tracking using OPAL-cycl.

Left plot of Figure 5 shows the accelerating orbit of reference particle. After 106 turns, the energy increases from 870 keV at the injection point to 72.16 MeV at the deflection point.

From theoretic view, there should be an eigen ellipse for any given energy in stable area of a fixed accelerator structure. Only when the initial phase space shape matches its eigen ellipse, the oscillation of beam envelop amplitude will get minimal and the transmission efficiency get maximal. We can calculate the eigen ellipse by single particle tracking using betatron oscillation property of off-centered particle as following: track an off-centered particle and record its coordinates and momenta at the same azimuthal position for each revolution. Figure 6 shows the eigen ellipse at symmetric line of sector magnet for energy of 2 MeV in Injector II.

Right plot of Figure 5 shows very good agreement of the tune diagram by OPAL-cycl and FIXPO. The trivial discrepancy should come from the methods they used. In FIXPO, the tune values are obtained according to the crossing points of the initially displaced particle. Meanwhile, in OPAL-cycl, the Fourier analysis method is used to manipulate orbit difference between the reference particle and an initially displaced particle. The frequency point with the biggest amplitude is the betatron tune value at the given energy.

Following is the input file for single bunch tracking with space charge effects in Injector II.

Injector2-sc.in

For the supplementary files see above, PSI Injector II Cyclotron

To run OPAL on single node, just use this command:

To run OPAL on N nodes in parallel environment interactively, use this command instead:

If restart a job from the last step of an existing .h5 file, add a new argument like this:

Figure 7 and Figure 8 are simulation results, shown by H5root code.

From the view of numerical simulation, the difference between Injector II and Ring cyclotron comes from two aspects:

Figure 9 shows a single particle tracking result and tune calculation result in the PSI Ring cyclotron. Limited by size of the user guide, we don’t plan to show too much details as in Injector II.

The GUNGAUSSFLATTOPTH and ASTRAFLATTOPTH distribution types are two common types previously implemented to simulate electron beams emitted from photocathodes in an electron photoinjector. These are no longer explicitly supported and are instead now defined as specialized sub-types of the distribution type FLATTOP. That is, the emitted distributions represented by GUNGAUSSFLATTOPTH and ASTRAFLATTOPTH can now be easily reproduced by using the FLATTOP distribution type and we would encourage use of the new command structure.

Having said this, however, old input files that use the GUNGAUSSFLATTOPTH and ASTRAFLATTOPTH distribution types will still work as before, with the following exception. Previously, OPAL had a Boolean OPTION command FINEEMISSION (default value was TRUE). This OPTION is no longer supported. Instead you will need to set the distribution attribute Table 28 to a value that is 10 \(\times\) the value of the distribution attribute Table 26 in order for your simulation to behave the same as before.

The FROMFILE, GAUSS and BINOMIAL distribution types have changed from previous versions of OPAL. However, legacy distribution commands should work as before as long as the momentum units are converted to the current convention (see section on units).

Input momentum can be given without units i.e. as \(\beta\gamma\), or in eV/c. Up until OPAL 2.2, eV was used instead of eV/c, but this was changed since eV is a unit of energy rather than momentum. To adapt old files with momentum in eV, such that they can work for the newer OPAL versions, the following formula can be used:

and you will need to set the distribution attribute INPUTMOUNITS to:

We define the global cartesian coordinate system, also known as floor coordinate system with \(K\), a point in this coordinate system is denoted by \((X, Y, Z) \in K\). In Figure 10 of the accelerator is uniquely defined by the sequence of physical elements in \(K\). The beam elements are numbered \(e_0, \ldots , e_i, \ldots e_n\).

A local coordinate system \(K'_i\) is attached to each element \(e_i\). This is simply a frame in which \((0,0,0)\) is at the entrance of each element. For an illustration see Figure 10. The local reference system \((x, y, z) \in K'_n\) may thus be referred to a global Cartesian coordinate system \((X, Y, Z) \in K\).

The local coordinates \((x_i, y_i, z_i)\) at element \(e_i\) with respect to the global coordinates \((X, Y, Z)\) are defined by three displacements \((X_i, Y_i, Z_i)\) and three angles \((\Theta_i, \Phi_i, \Psi_i)\).

\(\Psi\) is the roll angle about the global \(Z\)-axis. \(\Phi\) is the yaw angle about the global \(Y\)-axis. Lastly, \(\Theta\) is the pitch angle about the global \(X\)-axis. All three angles form right-handed screws with their corresponding axes. The angles (\(\Theta,\Phi,\Psi\)) are the Tait-Bryan angles [12].

The displacement is described by a vector \(\mathbf{v}\) and the orientation by a unitary matrix \(\mathcal{W}\). The column vectors of \(\mathcal{W}\) are unit vectors spanning the local coordinate axes in the order \((x, y, z)\). \(\mathbf{v}\) and \(\mathcal{W}\) have the values:

where

We take the vector \(\mathbf{r}_i\) to be the displacement and the matrix \(\mathcal{S}_i\) to be the rotation of the local reference system at the exit of the element \(i\) with respect to the entrance of that element.

Denoting with \(i\) a beam line element, one can compute \(\mathbf{v}_i\) and \(\mathcal{W}_i\) by the recurrence relations

where \(\mathbf{v}_0\) corresponds to the origin of the LINE and \(\mathcal{W}_0\) to its orientation. In OPAL-t they can be defined using either X, Y, Z, THETA, PHI and PSI or ORIGIN and ORIENTATION, see Simple Beam Lines.

In order to calculate space charge in the electrostatic approximation, we introduce a co-moving coordinate system \(K_{\text{sc}}\), in which the origin coincides with the mean position of the particles and the mean momentum is parallel to the z-axis.

In order to compute statistics of the particle ensemble, \(K_s\) is introduced. The accompanying tripod (Dreibein) of the reference orbit spans a local curvilinear right handed system \((x,y,s)\). The local \(s\)-axis is the tangent to the reference orbit. The two other axes are perpendicular to the reference orbit and are labelled \(x\) (in the bend plane) and \(y\) (perpendicular to the bend plane).

Both coordinate systems are described in Figure 11.

The reference orbit consists of a series of straight sections and circular arcs and is computed by the Orbit Threader i.e. deduced from the element placement in the floor coordinate system.

To facilitate the change for users we will provide a compatibility mode. The idea is that the user does not have to change the input file. Instead OPAL-t will compute the positions of the elements. For this it uses the bend angle and chord length of the dipoles and the position of the elements along the trajectory. The user can choose whether effects due to fringe fields are considered when computing the path length of dipoles or not. The option to toggle OPAL-t’s behavior is called IDEALIZED. OPAL-t assumes per default that provided ELEMEDGE for elements downstream of a dipole are computed without any effects due to fringe fields.

Elements that overlap with the fields of a dipole have to be handled separately by the user to position them in 3D.

We split the positioning of the elements into two steps. In a first step we compute the positions of the dipoles. Here we assume that their fields don’t overlap. In a second step we can then compute the positions and orientations of all other elements.

The accuracy of this method is good for all elements except for those that overlap with the field of a dipole.

The OrbitThreader integrates a design particle through the lattice and setups up a multi map structure (IndexMap). Furthermore when the reference particle hits an rf-structure for the first time then it auto-phases the rf-structure, see Appendix Auto-phasing Algorithm. The multi map structure speeds up the search for elements that influence the particles at a given position in 3D space by minimizing the looping over elements when integrating an ensemble of particles. For each time step, IndexMap returns a set of elements \(\mathcal{S}_{\text{e}} \subset {e_0 \ldots e_n}\) in case of the example given in Figure 10. An implicit drift is modelled as an empty set \(\emptyset\).

This file is used to log the statistical properties of the bunch in the ASCII variant of the SDDS format [13]. It can be viewed with the SDDS Tools [14] or GNUPLOT. The frequency with which the statistics are computed and written to file can be controlled With the option STATDUMPFREQ (see Option Statement). The name of the output file is <input_file_name >.stat. The information that is stored are found in Table 5. Additionally, more statistical information about the beam can be saved. The option COMPUTEPERCENTILES allow whether the 68.27 (1 sigma for normal distribution), the 95.45 (2 sigmas), the 99.73 (3 sigmas) and the 99.994 (4 sigmas) percentiles for the beam size and the normalized emittance should be computed (see Table 6)). Whereas, DUMPBEAMMATRIX control whether to write the 6-dimensional beam matrix (upper triangle only) to statatistic file (see Table 7)).

OPAL-t computes the statistics of the bunch for every MONITOR that it passes. The name of the output file is data/<input_file_name >_Monitors.stat. The information that is written can be found in the following table.

OPAL-t copies the input file into this file and replaces all occurrences of ELEMEDGE with the corresponding position using X, Y, Z, THETA, PHI and PSI. The name of the output file is data/<input_file_name >_3D.opal.

The trajectory of the reference particle is stored in this ASCII file. The name of the output file is data/<input_file_name >_DesignPath.dat. The content of the columns are listed in the following table.

(for a straight multipole)

Most accelerator modeling codes use the hard-edge model for magnets - constant Hamiltonian. Real magnets always have a smooth transition at the edges - fringe fields. To obtain a multipole description of a field we can apply the theory of analytic functions.

Assuming that \(A\) has only a non-zero component \(A_s\) we get

These equations are just the Cauchy-Riemann conditions for an analytic function \(\tilde{A} (z) = A_s (x, y) + i V(x,y)\). So the complex potential is an analytic function and can be expanded as a power series

with \(\lambda_n, \mu_n\) being real constants. It is practical to express the field in cylindrical coordinates \((r, \varphi, s)\)

From the real and imaginary parts of equation () we obtain

Taking the gradient of \(-V(r, \varphi)\) we obtain the multipole expansion of the azimuthal and radial field components, respectively

Furthermore, we introduce the normal multipole coefficient \(b_n\) and skew coefficient \(a_n\) defined with the reference radius \(r_0\) and the magnitude of the field at this radius \(B_0\) (these coefficients can be a function of s in a more general case as it is presented further on).

To obtain a model for the fringe field of a straight multipole, a proposed starting solution for a non-skew magnetic field is

It is straightforward to derive a relation between coefficients

By identifying the term in front of the same powers of \(r\) we obtain the recurrence relation

The solution of the recursion relation becomes

Therefore

The transverse components of the field are

where the following gradients determine the entire potential and can be deduced from the function \(C_{n,0}(z)\) once the harmonic \(n\) is fixed.

The z-directed component of the filed can be expressed in a similar form

The gradient functions \(g_{rn}, g_{\varphi n}, g_{zn}\) are obtained from

One preferred model to approximate the gradient profile on the central axis is the k-parameter Enge function

where \(d(z)\) is the distance to the field boundary and \(\lambda\) characterizes the fringe field length.

(fixed radius)

We consider the Frenet-Serret coordinate system \( ( \hat{\mathbf{x}}, \hat{\mathbf{s}}, \hat{\mathbf{z}} )\) with the radius of curvature \( \rho \) constant and the scale factor \( h_s = 1 + x/ \rho\). A conversion to these coordinates implies that

To simplify the problem, consider multipoles with mid-plane symmetry, i.e.

The most general form of the expansion is

Maxwell’s equations \( \nabla \cdot \mathbf{B} = 0 \) and \( \nabla \times \mathbf{B} = 0 \) in the above coordinates yield

The substitution of General form into Maxwell’s equations allows for the derivation of recursion relations. Maxwell equations gives

Equating the powers in \(x^i z^{2k}\)

A similar result is obtained from Maxwell equations

The last equation from \(\nabla \times \mathbf{B} = 0\) should be consistent with the two recursion relations obtained. Maxwell equations implies

This results follows directly from a factors and c factors; therefore the relations are consistent. Furthermore, the last required relations is obtained from the divergence of B

Using the relation (a factors) to replace the \(a\) coefficients with \(b\)’s we arrive at

All the coefficients above can be determined recursively provided the field \(B_z\) can be measured at the mid-plane in the form

where \(B_{i,0}\) are functions of \(s\) and they can model the fringe field for each multipole term \(x^n\). As an example, for a dipole magnet, the \(B_{1,0}\) function can be model as an Enge function or \(tanh\).

(variable radius of curvature)

The difference between this case and the above is that \(\rho\) is now a function of \(s\), \(\rho(s)\). We can obtain the same result starting with the same functional forms for the field (General form). The result of the previous section also holds in this case since no derivative \(\frac{\partial}{\partial s}\) is applied to the scale factor \(h_s\). If the radius of curvature is set to be proportional to the dipole filed observed by some reference particle that stays in the centre of the dipole

This is a different, more compact treatment The derivation is more clear if we gather the variables together in functions. We assume again mid-plane symmetry and that the z-component of the field in the mid-plane has the form

where \(T(s)\) is the transverse field profile and \(S(s)\) is the fringe field. One of the requirements of the symmetry is that \(B_z(z) = B_z(-z)\) which using a scalar potential \(\psi\) requires \(\frac{\partial \psi}{\partial z}\) to be an even function in z. Therefore, \(\psi\) is an odd function in z and can be written as

The given transverse profile requires that \(f_0(x,s) = T(x)S(s)\), while \(\nabla^2 \psi = 0\) follows from Maxwell’s equations as usual, more explicitly

For a straight multipole \(h_s = 1\). Laplace’s equation becomes

By equating powers of \(z\) we obtain the recursion relation

The general expression for any \(f_n(x,s)\) is then obtained from the mid-plane field by

For a curved multipole of constant radius \(h_s = 1 + \frac{x}{\rho} \quad \text{with} \quad \rho = const.\) The corresponding Laplace’s equation is

Again we substitute with the functional form of the potential to get the recursion

Finally consider what changes for \(\rho = \rho (s)\). Laplace’s equation is

The last step is again the substitution to get

If the radius of curvature is proportional to the dipole field in the centre of the multipole (the dipole component of the transverse field is a constant \(T_{dipole}(x) = B_0\) and

This expression can be replaced in ([eq.40]) to obtain a more explicit equation.

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:

where the file PartDatabase.dat should have two lines:

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 Units.

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.

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:

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.

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 Multi-bunch Mode.

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 Variables in OPAL-cycl). 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 \((X, Y, Z, PX, PY, PZ)\) using the starting coordinates \(r_0\) (RINIT), \(\phi_0\) (PHIINIT), and \(z_0\) (ZINIT), and the starting momenta \(p_{total}\) (defined by the beam energy specified in the Beam Command), \(p_{r0}\) (PRINIT), and \(p_{z0}\) (PZINIT) of the reference particle, defined in the CYCLOTRON element. Note that \(p_{\phi 0}\) is calculated automatically from \(p_{total}\), \(p_{r0}\), and \(p_{z0}\).

If TYPE=CARBONCYCL, the program requires the \(B_z\) data which is stored in a sequence shown in Figure 14.

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:

The number of values per column is arbitrary.

If TYPE=CYCIAE, the program requires data format given by ANSYS10.0. This function is originally for the 100 MeV cyclotron of CIAE, whose isochronous fields is numerically computed 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):

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 if the decimal is unlimited. For instance, if \(\Delta \theta = \frac{1}{3} {^{\circ}}\), the fourth line of the header should be -3.0. In other words, the fields are stored in the following format:

BANDRF fieldmap TYPE allows to read both the magnetic field and the electric field from the CYCLOTRON element. The magnetic field data (\(B_z\)) are stored in the same format as CARBONCYCL. Regarding the electric field, the 3D RF field-map can be read from H5Hut type file (.h5part file extension) making use of a conversion tool integrated into OPAL: ascii2h5block. For the details about its usage, please see Read 3D RF field-map. The electric field map can be imported from some field files, which can be optionally superposed, providing different meshing specifications for the field map. This is useful for modeling the central region electric fields which usually has complicated shapes. Please note that in this case, the E field is treated as a part of CYCLOTRON element, rather than an independent RFCAVITY element.

If TYPE=RING, 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 if the decimal is unlimited. For instance, if \(\Delta \theta = \frac{1}{3}{^{\circ}}\), the fourth line of the header should be -3.0.

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.

The field-map with file name fdf-tosca-field-map.table is loaded, which is a file like

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)\).

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 Cyclotron.

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 Cyclotron.

The voltage profile of a cavity gap is read from ASCII file. Here is an example:

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.

The 3D RF field-map can be read from H5Hut type file for TYPE=BANDRF. To generate the h5part field files, it is necessary to use the ascii2h5block conversion tool. This tool is available after compiling OPAL with ENABLE_BANDRF=ON. This program builds h5part field files from ASCII file output of ANSYS or OPERA 3D field data. The command to use it should be like:

The ASCII input files must be adapted for an adequate conversion. The first row of each input field map that you wish to combine should specify the grid size in each direction, that is integers with the amount of steps (or different values) in \(x\), \(y\) and \(z\). The first three columns contains the spatial coordinates specified in meters, and the last three columns store the field data. Positions and field are supplied on points in a Cartesian coordinates grid:[(x, y, z)]. Example showing the requiered format for ASCII file:

See OPAL-t Statistics output. The only difference is in the cyclotron coordinate convention. The z (or s) coordinate refers to the vertical direction and y corresponds to the longitudinal or beam travel direction. In addition, OPAL-cycl includes some extra data, as is described in the following table.

The intermediate phase space data is stored in an ASCII file which can be used to 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.

The phase space data per STEPSPERTURN (a parameter in the TRACK command) steps is stored in an ASCII file. The file’s is named <input_file_name >-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 STEPSPERTURN steps; Otherwise, the particle may not go through a full 360\(^{\circ}\) in STEPSPERTURN 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.

The tunes \(\nu_r\) and \(\nu_z\) of each energy are stored in a ASCII file named tuningresult.

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 [19] and can be analyzed using H5root [20], [21]. The frequency of the data output can be set using the PSDUMPFREQ option of OPTION statement. The file is named <input_file_name >.h5. This output can be switched on or off with the option ENABLEHDF5.

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 is named <input_file_name >-trackOrbit.dat. The frequency of the data output can be set using the SPTDUMPFREQ option of OPTION statement.

The magnetic field map read by OPAL could be saved into an output file called gnu.out. This file stores the magnetic field \(B_z\) data in the median plane in cylindrical coordinates following the sequence shown in Figure 14. This output file is available in the cyclotron field types CARBONCYCL, BANDRF, FFA and AVFEQ.

In case of CARBONCYCL or BANDRF an additional output field file an additional output field file is saved (eb.out), storing three vectors with the position, electric field and the magnetic field.

Writing these output files is enabled when OPAL is run on single node, and it can be switched on or off with the INFO option of OPTION statement.

Moreover, DUMPFIELDS and DUMPEMFIELDS routines (see Chapter Field Output) can be employed to write out the external field used in OPAL-cycl in a different grid than the input field maps.

Simulation of the PSI Ring Cyclotron at 580 MeV and current 2.2 mA. 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.

The creation of the element map is based on the Hamiltonian Mechanic, more specifically on the Lie Operator.

Here \(H\), the time dependent Hamiltonian represents the total energy, consisting of the kinetic \(T\) and potential \(V\) energy. The lower equations are the Hamiltonian equations of motion, where the momenta \(p_i\) and the positions \(q_i\) form the canonical pairs. Using canonical transformations, the Hamiltonian can be adjusted to use the path length \(s\) as independent and the particle parameters as dependent variable(s).

Introducing the Lie operator, which acts similar to a "waiting" Poisson Bracket:

If a function \(f:= f \left( \mathbf{q_{(s)}},\mathbf{p_{(s)}}\right) \) describes one of the phase space variables \(v\) its total derivative to the independent variable, combined with the Hamiltonian equations of motions, is similar to the Lie operator times the indepedent variable, i. e. \(s\).

The integral over the independent variable \(\int \cdot ds\):

Where \(e^{-:\!H\!: s}\) is the Lie expansion.

For the derivative of the Hamiltonian, a Differential Algebra (DA) was used, in particular the Truncated Power Series Algebra (TPSA). This algebra uses the Taylor expansion as the equivalent function, which also is responsible for its name by creating truncated power series. Just form the definition of the Taylor expansion, it can be seen that a finite, to the order \(\Omega\), expansion is an approximation of the actual function, due to the error term \(\mathcal{O}\left( \mathbf{v}^{\,\Omega +1}\left( \Delta s\right)\right) \).

In OPAL-map the Hamiltonian gets Taylor expanded and its map derived (Creation of Map) in the TPSA using the OPAL DA package. The truncation length gets defined in TRACK setting the MAP_ORDER attribute.

OPAL-map is the new flavour of OPAL and currently just contains the fundamental beam line elements:

The file structure is analogous to OPAL-t Statistics output file.

OPAL-map computes the dispersion of the beam line according to:

where \(R\) is the Transfermatrix and \(\eta\) is the dispersion.

Here the transfer map of the whole beamline gets printed. Its format is analogous to a DA- FVPs (src/Classic/FixedAlgebra), where the polynomials (DA- FTPs) get represented in a 6 dimensional vector. The first number describes the coefficient and the following 6 the variables of the monomial.

OPAL recognizes some specific string names using a list of predefined strings. Consequentially, attributes of type PredefinedString should only accept strings that are contained in a predefined set of strings according to their internal description. If the provided string is not contained in the predefined set, an exception is thrown, since this string is marked as a protected keyword and cannot be used to name an object. Most of the string attributes in OPAL have been defined based on PredefinedString to prevent the user from entering variable names that are used internally for another functionality. All these predefined string attributes do not have to be enclosed in quotes, unlike regular strings. The list of accepted strings as well as the default value, if any, are added to the --help-command message (see Table 4).

Numerical values are entered like FORTRAN constants. Real values are accepted in INTEGER or REAL format. The use of a decimal exponent, marked by the letter D or E, is permitted.

Examples:

OPAL recognizes a few built-in mathematical and physical constants (see <<tab_constant>). Their names must not be used for user-defined labels. Additional symbolic constants may be defined to simplify their repeated use in statements and expressions.

Here the RANK represents the MPI-Rank of the process and \(N_{p}\) the total number of MPI processes.

Often a set of numerical values depends on a common variable parameter. Such a variable must be defined as a global variable by one of

When such a variable is used in an expression, OPAL uses the current value of the variable. When the value is a constant or an expression preceded by the delimiter = it is evaluated immediately and the result is assigned to the variable as a constant. When the value is an expression preceded by the delimiter := the expression is retained and re-evaluated whenever one of its operands changes.
Example:

When the value of X is changed, the lengths of the drift spaces are recalculated as X and 2-X respectively.

In arithmetic expressions the attributes of physical elements or commands can occur as operands. They are named respectively by

If they are arrays, they are denoted by

Values are assigned to attributes in element definitions or commands.

Example:

D1→L refers to the length L of the drift space D1.

Definition of random machine imperfections requires evaluation of expressions containing random functions. These are not evaluated like other expressions before a command begins execution, but sampled as required from the distributions indicated when errors are generated. Such an expression is known as a deferred expression. Its value cannot occur as an operand in another expression.

A place denotes a single element, or the position following that element. It can be specified by one of the choices

Either form may be qualified by one or more beam line names, as described by the formal syntax:

An omitted index defaults to one. Examples: assume the following definitions:

The line L is equivalent to the sequence of elements

Some possible place definitions are:

A range in a beam line (see Chapter Beam Lines) is selected by the following syntax:

This denotes the range of elements from the first`place` to the second place. Both positions are included. A few special cases are worth noting:

Examples: Assume the following definitions:

The line L is equivalent to the sequence of elements

Examples for range selections:

For the time being, logical arrays can only be given as a list. The formal syntax is:

Example:

Real arrays have the following syntax:

Example:

There are also three functions allowing the generation of real arrays:

String arrays can only be given as lists of single values. For permissible values String values see String Attributes.

Example:

Token list arrays are always lists of single token lists.

Example:

A user who is uncertain about the attributes of a command should try the command HELP, which has three formats:

label is an identifier (see Identifiers or Labels). If it is non-blank, OPAL prints the function of the object label and lists its attribute types. Entering HELP alone displays help on the HELP command.

Examples:

OPAL recognizes several types of variables.

OPAL recognizes a few build-in built-in mathematical and physical constants (see Table 19). Additional constants can be defined by the command

which defines a constant with the name label. The keyword REAL is optional, and label must be unique. An existing symbolic constant can never be redefined. The real-expression is evaluated at the time the CONST definition is read, and the result is stored as the value of the constant.

Example:

A vector of expressions is established by a statement

The keyword REAL is optional. It creates a new global vector vector-name and discards any old vector with the same name. Its value depends on all quantities occurring in vector-expression (see Arrays). Whenever an operand changes in vector-expression, a new value is calculated. The definition may be thought of as a mathematical equation. However, OPAL is not able to solve the equation for a quantity on the right-hand side.

Example:

Circular definitions are not allowed.

A value is assigned to a variable or vector by using the function EVAL(real- expression). When seen, this function is immediately evaluated and replaced by the result treated like a constant.

This statement acts like a FORTRAN or C assignment. The real-expression or vector-expression is evaluated, and the result is assigned as a constant to the variable or vector on the left-hand side. Finally the expression is discarded. The EVAL function can also be used within an expression, e. g.:

A sequence like the following is permitted:

The statement

evaluates a set of expressions using the most recent values of any operands and prints the results on the standard error file.

Example:

These commands give the results:

This commands serves mainly for printing one or more quantities which depend on matched attributes. It also allows use of OPAL as a programmable calculator. One may also tabulate functions.

The ECHO statement has two formats:

message is a string value (see String Attributes). It is immediately transmitted to the ECHO stream.

The command SYSTEM allows to execute operating system commands on a single core. After execution of the system command, successful or not, control returns to OPAL. It has two formats:

The string (see String Attributes) string must be a valid operating system command.

Most UNIX commands can be issued directly.

Example:

causes a listing of the current directory in long form on the terminal.

The SYSTEM command can also be used to execute an external script.

Below is an example that generates a random seed with python. The python script gen_seed.py is called with the SYSTEM command. This scripts writes the seed to a file gen_seed.opal that can be read with the CALL command.

Also the command PSYSTEM allows to execute operating system commands. But contrary to the command SYSTEM this command is executed on all cores.

The CALL command has two formats:

file-name is a string (see String Attributes). The statement causes the input to switch to the named file. Input continues on that file until a STOP or an end of file is encountered. Example:

Before entering the main input file, the file $HOME/init.opal is executed (if present). It can be thought of as a CALL, FILE=$HOME/init.opal statement at the start of the input file. This init file can be used to avoid writing the same lines in every input file. For example one can set some OPTION, define a MACRO, or define some physics constants.

Note that one should be well aware of the contents of the init.opal file as otherwise it might lead to confusion when switching to a different computer where the file is not present. Therefore, it might be advisable to use an explicit CALL statement instead.

Example:

An RBEND is a rectangular bending magnet. The key property of an RBEND is that it has parallel pole faces. Figure 16 shows an RBEND with a positive bend angle and a positive entrance edge angle.

An RBEND3D3D is a rectangular bending magnet. The key property of an RBEND3D is that it has parallel pole faces. Figure 16 shows an RBEND3D with a positive bend angle and a positive entrance edge angle.

Figure 17 illustrates a general sector bend(SBEND) with a positive bend angle \(\alpha\). In this example the entrance and exit edge angles \(E_{1}\) and \(E_{2}\) have positive values. For a positively charge particle, the magnetic field is directed out of the page.

An SBEND is a sector bending magnet. An SBEND can have independent entrance and exit edge angles. Figure 17 shows an SBEND with a positive bend angle, a positive entrance edge angle, and a positive exit edge angle.

Describing an RBEND or an SBEND in an OPAL-t simulation requires effectively identical commands. There are only slight differences between the two. The L attribute has a different definition for the two types of bends (see RBend (OPAL-t) and SBend (OPAL-t)), and an SBEND has an additional attribute E2 that has no effect on an RBEND (see SBend (OPAL-t)). Therefore, in this section, we will give several examples of how to implement a bend, using the RBEND and SBEND commands interchangeably. The understanding is that the command formats are essentially the same.

When implementing an RBEND or SBEND in an OPAL-t simulation, it is important to note the following:

For the rest of this section, we will give several examples of how to input bends in an OPAL-t simulation. We will start with a simple example using the ANGLE attribute to set the bend strength and using the default field map (see Default Field Map (OPAL-t)) for describing the magnet fringe fields (see 1DProfile1):

This is a definition of a simple RBEND that bends the beam in a positive direction 30 degrees (towards the negative x axis as if Figure 16). It has a design energy of 10 MeV, a length of 0.5 m, a vertical gap of 2 cm and a 0\(^{\circ}\) entrance edge angle. (Therefore the exit edge angle is 30\(^{\circ}\).) We are using the default, internal field map "1DPROFILE1-DEFAULT" see Default Field Map (OPAL-t) which describes the magnet fringe fields see 1DProfile1. When OPAL is run, you will get the following output (assuming an electron beam) for this RBEND definition:

The first section of this output gives the properties of the reference trajectory like that described in Figure 16. From the value of ANGLE and the length, L, of the magnet, OPAL calculates the 10 MeV reference particle trajectory radius, R. From the bend geometry and the entrance angle (0\(^{\circ}\) in this case), the exit angle is calculated.

The second section gives the field amplitude of the bend and its gradient (quadrupole focusing component), given the particle charge (\(-e\) in this case so the amplitude is negative to get a positive bend direction). Also listed is the rotation of the magnet about the various axes.

Of course, in the actual simulation the particles will not see a hard edge bend magnet, but rather a soft edge magnet with fringe fields described by the RBEND field map file FMAPFN (see 1DProfile1). So, once the hard edge bend/reference trajectory is determined, OPAL then includes the fringe fields in the calculation. When the user chooses to use the default field map, OPAL will automatically adjust the position of the fringe fields appropriately so that the soft edge magnet is equivalent to the hard edge magnet described by the reference trajectory. To check that this was done properly, OPAL integrates the reference particle through the final magnet description with the fringe fields included. The result is shown in the final part of the output. In this case we see that the soft edge bend does indeed bend our reference particle through the correct angle.

What is important to note from this first example, is that it is this final part of the bend output that tells you the actual bend angle of the reference particle.

In this next example, we merely rewrite the first example, but use K0 to set the field strength of the RBEND, rather than the ANGLE attribute:

The output from OPAL now reads as follows:

The output is effectively identical, to within a small numerical error.

Now, let us modify this first example so that we bend instead in the negative x direction. There are several ways to do this:

1.

2.

3.

4.

In each of these cases, we get the following output for the bend (to within small numerical errors).

In general, we suggest to always define a bend in the positive x direction (as in Figure 16) and then use the ROTATION attribute to bend in other directions in the x/y plane (as in examples 2 and 4 above).

As a final RBEND example, here is a suggested format for the four bend definitions if one where implementing a four dipole chicane:

Up to now, we have only given examples of RBEND definitions. If we replaced "RBend" in the above examples with "SBend", we would still be defining valid OPAL-t bends. In fact, by adjusting the L attribute according to RBend (OPAL-t) and SBend (OPAL-t), and by adding the appropriate definitions of the E2 attribute, we could even get identical results using `SBEND`s instead of `RBEND`s. (As we said, the two bends are very similar in command format.)

Up till now, we have only used the default field map. Custom field maps can also be used. There are two different options in this case see 1DProfile1:

The first case describes how field maps were used in previous versions of OPAL (and can still be used in the current version). The second option is new to OPAL OPALversion 1.2.00 and it has a couple of advantages:

We will now give examples of how to use a custom field map, starting with the first case where the field map describes the fringe fields and the magnet length. Assume we have the following 1DProfile1 field map:

We can use this field map to define the following bend (note we are now using the SBEND command):

Notice that we do not set the magnet length using the L attribute. (In fact, we don’t even include it. If we did and set it to a non-zero value, the exit fringe fields of the magnet would not be correct.) This input gives the following output:

Because we set the bend strength using the ANGLE attribute, the magnet field strength is automatically adjusted so that the reference particle is bent exactly ANGLE radians when the fringe fields are included. (Lower output.)

Now we will illustrate the case where the magnet length is set by the L attribute and only the fringe fields are described by the field map. We change the TEST-MAP.T7 file to:

and change the bend input to:

This results in the same output as the previous example, as we expect.

As a final example, let us now use the previous field map with the following input:

Instead of setting the bend strength using ANGLE, we use K0. This results in the following output:

In this case, the bend angle for the reference trajectory in the first section of the output no longer matches the reference trajectory bend angle from the lower section (although the difference is small). The reason is that the path of the reference particle through the real magnet (with fringe fields) no longer matches the ideal trajectory. (The effective length of the real magnet is not quite the same as the hard edged magnet for the reference trajectory.)

We can compensate for this by changing the field map file TEST-MAP.T7 file to:

We have moved the Enge function origins (see 1DProfile1) outward from the entrance and exit faces of the magnet by 0.3026 mm. This has the effect of making the effective length of the soft edge magnet longer. When we do this, the same input:

produces

Now we see that the bend angle for the ideal, hard edge magnet, matches the bend angle of the reference particle through the soft edge magnet. In other words, the effective length of the soft edge, real magnet is the same as the hard edge magnet described by the reference trajectory.

So far we have described how to setup an RBEND or SBEND element, but have not explained how OPAL-t uses this information to calculate the magnetic field. The field of both types of magnets is divided into three regions:

This can be seen clearly in Figure 45.

The purpose of the 1DProfile1 field map (see 1DProfile1) associated with the element is to define the Enge functions (Enge function) that model the entrance and exit fringe fields. To model a particular bend magnet, one must fit the field profile along the mid-plane of the magnet perpendicular to its face for the entrance and exit fringe fields to the Enge function:

where \(D\) is the full gap of the magnet, \(N_{order}\) is the Enge function order and \(z\) is the distance from the origin of the Enge function perpendicular to the edge of the dipole. The origin of the Enge function, the order of the Enge function, \(N_{order}\), and the constants \(c_0\) to \(c_{N_{order}}\) are free parameters that are chosen so that the function closely approximates the fringe region of the magnet being modeled. An example of the entrance fringe field is shown in Figure 18.

Let us assume we have a correctly defined positive RBEND or SBEND element as illustrated in Figure 16 and Figure 17. (As already stated, any bend can be described by a rotated positive bend.) OPAL-t then has the following information:

Here, we have defined an overall Enge function, \(F(z)\), with three parts: entrance, center and exit. The exit and entrance fringe field regions have the form of Enge function with parameters defined by the 1DProfile1 field map file given by the element parameter FMAPFN. Defining the coordinates:

using the conditions

and making the definitions:

we can expand the field off axis, with the result:

These expression are not well suited for numerical calculation, so, we expand them about \(y\) to \(O(y^2)\) to obtain:

These are the expressions OPAL-t uses to calculate the field inside an RBEND or SBEND. First, a particle’s position inside the bend is determined (entrance region, center region, or exit region). Depending on the region, OPAL-t then determines the values of \(\Delta_x\), \(y\) and \(\Delta_z\), and then calculates the field values using the above expressions.

Rather than force users to calculate the field of a dipole and then fit that field to find Enge coefficients for the dipoles in their simulation, we have a default set of values we use from [24] that are set when the default field map, 1DPROFILE1-DEFAULT is used:

The same values are used for both the entrance and exit regions of the magnet. In general they will give good results. (Of course, at some point as a beam line design becomes more advanced, one will want to find Enge coefficients that fit the actual magnets that will be used in a given design.)

The default field map is the equivalent of the following custom 1DProfile1 (see 1DProfile1 for an explanation of the field map format) map:

As one can see, the default magnet gap for 1DPROFILE1-DEFAULT is set to 2.0 cm. This value can be overridden by the GAP attribute of the magnet (see RBend (OPAL-t) and SBend (OPAL-t)).

The SBEND3D element enables definition of a bend from 3D field maps. This can be used in conjunction with the RINGDEFINITION element to make a ring for tracking through OPAL-cycl.

Field maps are defined using Cartesian coordinates but in a polar geometry. The following conventions have to be fulfilled:

SBEND3D assumes a dipole symmetry. In a dipole symmetry, fields below the symmetry plane Z=0 have the same field in the direction perpendicular to the symmetry plane, \(B_{z}\), but field components parallel to the symmetry plane have the opposite direction (sign).

Below are two examples of a SBEND3D which loads a field map file named “90degree_Dipole_Magnet.out” defining a hard edge model of 90\(^{\circ}\) dipole magnet with homogenous magnetic field. The first 8 lines are presumed to be header material and are ignored. Positions have units of m and fields units of Tesla. The corresponding 3D magnetic field map is shown in the following figure in the Cartesian coordinate system (x, y, z). A horizontal cross section of the 3D magnetic field map when z = 0 is also shown.

The first few lines of the field map file are as follows:

This is a restricted feature for OPAL-cycl.