A Super-Fast Free-Electron Laser Simulation Code for Online Optimization

Abstract

The X-ray free-electron lasers (FELs) have stimulated the growing interest of researchers in different fields. This gives rise to an increasing simulation work of design and optimization of FEL facilities and demonstrations of novel FEL ideas. Most of the multi-dimensional simulation codes in use require large computational resources, while one-dimensional simulation codes can merely give an acceptable description of the FEL amplification process. This paper presents the development of a super-fast time-dependent FEL simulation code, which is mainly designed for seeded FEL and its application on a seed FEL facility. The diffraction factors have been introduced in order to model the effects of radiation field diffraction and drastically simplified the working equations. Meanwhile, a specially designed module for seeded FEL has been added in order to study the laser-electron beam interactions with arbitrarily parameters. The code can also provide a fairly accurate, quasi-real-time assistant tool for online optimization.

1. Introduction

The X-ray free-electron laser (FEL) has been recognized as a unique X-ray source of choice for experts when exploring the structures and dynamics of atomic and molecular systems at femtosecond time scales. It has been widely adopted in the fields of physics [1], chemistry [2], and biology [3], due to its ultra-high intensity and ultra-fast time structure. With growing interests in FEL applications, increasing FEL facilities have been constructed and novel FEL ideas have been proposed in the last decades, leading to a growing demand for FEL simulations.

An increasing number of numerical models and simulation codes of the FEL have been developed since the late 1990s. The commonly used multi-dimensional FEL codes can be roughly categorized as slowly varying envelope approximation (SVEA) simulation codes (GENESIS 1.3 [4], GINGER [5], FAST [6], RON [7], SIMPLEX [8], MINERVA [9]), and particle-in-cell (PIC) simulation codes (Puffin [10]). The SVEA assumes that the radiation field is represented by a slowly-varying amplitude and phase with a fast oscillatory term. The averaging of the field envelope will obviate the need for calculating this rapid sinusoidal oscillation and simplifying the field equation. In contrast, the PIC codes do not average the Maxwell’s equation over the undulator period. Instead, it integrates the Newton–Lorentz equations for the particles and it is considered to be a more fundamental model of FEL physics. It is capable of analyzing effects with a broader bandwidth that cannot be modelled easily in SVEA, but it requires more computational resources.

While these codes have been used successfully to model the FEL process and they have been extensively benchmarked against each other and experiments [11,12,13], nearly all of them are rather time-consuming when optimizing FEL performances with the method of advanced optimization algorithms [14,15]. These methods require performing FEL simulations over ten thousand times, which would take couple of days or even longer when time-dependent results are expected. SIMPLEX is comparatively fast and user-friendly, but it is still somewhat of inconvenient since it is packaged.

Some one-dimensional (1D) simulation codes have been developed in order to reduce computation time (Perave [16], Perseo [17], AURORA [18], and so on [19,20]). Nevertheless, most of them omit significant effects of the diffraction of radiation, leading to an unacceptable increase of FEL power.

In this paper, the development of a super-fast time-dependent FEL simulation code, named FALCON, and its applications on the Shanghai soft X-ray FEL facility (SXFEL), have been reported. The code is especially designed for seeded FEL. The radiation field is calculated by 1D equation with the consideration of diffraction effects. This significantly reduces the computing requirements and it still gives an acceptable description of the FEL amplification process. Another important motivation in this code is to model different FEL schemes. Specifically, the laser-electron beam interaction module enables simulations for various kinds of seeded FEL schemes with arbitrary user-defined laser parameters, such as the incidence angle, transverse, and longitudinal distributions. The time-dependence effects are analyzed using the ’slice’ model and the electron slips one wavelength behind the EM wave while traveling one undulator period.

It should be noted that the primal aim for this simulation code is to provide a fairly accurate, flexible, and computational-time-saving tool that can be used to give guidance on the online optimization of FEL facilities or fast optimization of FEL parameters and to give an estimation of the FEL performance for theoretical researches.

2. Methods and Mathematical Model

2.1. The Initial Particle Generation from Shot-Noise

Macro particles are introduced in order to represent the electron beam in all dimensions of the 6D phase space ($x,y,p_x,p_y,\theta ,\gamma$). For convenience, the transverse momenta are normalized to $mc$ while assuming $v_z=c$.

Although a random distribution would be the most straightforward way to generate initial particles, a loading strategy is typically used to reduce the number of macro particles. It can reduce the computing time while keeping the numerical noise under the acceptable level.

The initial particle phase space is generated by the Hammersley sequences [21]. The electron phase is uniformly generated, while the electron energy and transverse variables are subject to Gaussian distribution by adopting the ’Box-Muller’ transformation [22]. With this method, the correlation between the longitudinal variables, electron phase ($\theta$), and electron energy ($\gamma$) has been completely removed.

In order to correctly simulate the initial fluctuation in the longitudinal position, special care has to be taken for the modification of electron phase $\theta$ [23,24,25]. The most common method is to group several macro particles into a beamlet. Subsequently, a small variation in the evenly distributed longitudinal position $\theta$ is applied. In this simulation code, the quite loading algorithm is implemented, as in Equation (1)

$${\theta }_n={\theta }_{n,ini}+A_{n}sin({\theta }_{n,ini}+{\varphi }_n).$$

(1)

The indices presentted in ${\theta }_n$ on the right side of the equation represent the initial electron phase before introducing the loading strategy. A random shot-noise of its amplitude and phase, for the nth beamlet, is obtained by $A_n=2\sqrt{-ln\left(x_n\right)/N}$ and ${\varphi }_n=2\pi y_n$, where $x_n$ and $y_n$ are independent uniform random numbers and N denotes the total number of electrons.

2.2. The FODO Lattice and Undulator Tapering

Although the full FODO dynamics can be incorporated into numerical simulations, some additional approximations are made in order to simplify the integration process and improve the algorithmic efficiency.

The true dynamics in the FODO lattice can be satisfied by approximating the particle trajectories by the simply harmonic motion for smooth focusing, as implied by the fact that the transverse motion is governed by the second-order differential equations. For this reason, the effect of FODO lattice is approximated by a smooth focusing field whose focusing strength (oscillator frequency) $k_{\beta }$ equals to $1/{\overline{\beta }}_x$ [26], where the $\overline{\beta }$ is the average beta-function that is defined by the FODO structure. Although the electron trajectories under such an approximation can only roughly mimic those of the FODO lattice, it can give quite accurate predictions in terms of FEL radiation characteristics [26]. This is because these radiation characteristics rarely depend on a single electron, but rather on the average features of the entire electron beam.

The undulator tapering is another issue when it comes to the FEL simulations. Linear taper, quadratic taper, and fully customized taper are provided in FALCON.

2.3. Modulator and Dispersion Section

FALCON provides two methods for analyzing the energy modulation: ‘Fast sinusoidal modulation’ and ‘single-particle tracking modulation’. The latter one makes it possible to model the laser-electron beam interaction with an arbitrary incidence angle and field distributions.

2.3.1. Fast Sinusoidal Modulation

The effects of modulation and dispersion section on the phase space are specified as Equation (2) [27]:

$$\begin{array}{cc}\hfill p_j^{\prime }& =p_j+A·sin\left({\zeta }_j\right),\hfill \\ \hfill {\zeta }_j^{\prime }& ={\zeta }_j+B·p_j,\hfill \end{array}$$

(2)

where $p_j=({\gamma }_j-{\gamma }_0)/{\sigma }_{{\gamma }_0}$ denotes the dimensionless energy deviation of jth particle, ${\zeta }_j=2\pi {\theta }_j$ is the dimensionless longitudinal coordinate for jth particle, and $A=\Delta E/{\sigma }_E$ and $B=\left(R_{56}{\omega }_s{\sigma }_{\gamma }\right)/\left(c{\gamma }_0\right)$ denote the energy modulation amplitude and dimensionless dispersion strength, respectively.

2.3.2. Single-Particle Tracking Modulation

A special modulator section has been introduced in order to model the laser-electron beam interaction with arbitrary parameters. The energy exchange due to the electron-radiation interaction can be calculated as

$$\frac{d\gamma }{dz}=-\frac{e}{m{c}^2}\overrightarrow{\beta }·\overrightarrow{E_s},$$

(3)

where $\overrightarrow{\beta }$ and $\overrightarrow{E_s}$ are the normalized velocity of the electron beam and the transverse electric field of seed laser, respectively. The instantaneous change of electron velocity that is caused by Lorentz force is

$$\frac{d\overrightarrow{\beta }}{dz}=-\frac{e}{\gamma m{c}^2}[c\overrightarrow{\beta }\times (\overrightarrow{B_u}+\overrightarrow{B_s})+\overrightarrow{E_s}],$$

(4)

where $\overrightarrow{B_u}$ and $\overrightarrow{B_s}$ are the undulator magnetic field and transverse magnetic field of seed laser, respectively.

The general electric field of seed laser is $\overrightarrow{E_s}=E_x\overrightarrow{i_x}+E_y\overrightarrow{i_y}$. The transverse components, $E_x$ and $E_y$, can be written as [28]:

$$\begin{array}{cc}\hfill E_x& =E_{i}sin(k_s(zcos\varphi -ct)+{\Phi }_x),\hfill \\ \hfill E_y& =E_{i}sin(k_s(zcos\varphi +xsin\varphi -ct)+{\Phi }_y),\hfill \end{array}$$

(5)

where ${\Phi }_x=k_s(z_{0}cos\varphi +y_{0}sin\varphi )$, ${\Phi }_y=k_s(z_{0}cos\varphi +x_{0}sin\varphi )$, and $\varphi$ denotes the oblique incidence angle of the seed laser.

Instead of averaging over the undulator period, each undulator period has been divided into $N_n$ steps. Integrating over the entire undulator length with the Equations (3)–(5), various kinds of seeded FEL schemes with arbitrary user-defined laser parameters, such as the incidence angle [28], transverse, and longitudinal distributions (that are especially important for the modulation with few-cycle laser pulses), can be performed.

2.4. Radiation Diffraction

The main three-dimensional (3D) effect that tends to inhibit the FEL gain process over those in one-dimensional (1D) simulation is the radiation field diffraction. Because this code is mainly designed for seeded FEL, a simple algorithm has been implemented to solve this problem in FALCON.

For seeded FELs, the transverse distribution of the radiation field is generally Gaussian. The electric field in phasor notation is given by:

$$E(r,z)=E_0\frac{{\omega }_0}{\omega \left(z\right)}exp(\frac{-r^2}{{\omega }^2\left(z\right)})exp(-ik\frac{r^2}{2R\left(z\right)}-ikz+i\varphi \left(z\right)),$$

(6)

where $E_0$ denotes the original field amplitude that is calculated by the input power of radiation field, ${\omega }_0$ is the beam radius, $R\left(z\right)$ denotes the radius of curvature of the beam wavefronts at z, and $\omega \left(z\right)$ is the beam width at z, which is shown as Equation (7):

$$\omega \left(z\right)={\omega }_0\sqrt{1+{(\frac{z}{z_R})}^2},$$

(7)

and $z_R=\pi {\omega }_0^2/{\lambda }_s$ is the Rayleigh length of the seeding radiation field.

When the radiation field propagating one integration step size $\Delta z$, no field is generated by electrons and, therefore, the field propagates in free space. The field radius should increase by a factor of ${\chi }_r$, while the amplitude is expected to decrease by a factor of ${\chi }_{amp}$. The two factors are equivalent to

$$\begin{array}{cc}\hfill {\chi }_r& =\frac{\omega (z+\Delta z)}{\omega \left(z\right)}=\sqrt{\frac{z_R^2+{(z+\Delta z)}^2}{z_R^2+z^2}},\hfill \\ \hfill {\chi }_{amp}& =\frac{\tilde{E}(r,z+\Delta z)}{\tilde{E}(r,z)}=\frac{\omega \left(z\right)}{\omega (z+\Delta z)}exp(\frac{r^2}{{\omega }^2(z+\Delta z)}-\frac{r^2}{{\omega }^2\left(z\right)})=\frac{1}{{\chi }_r}\Gamma ,\hfill \end{array}$$

(8)

where

$$\Gamma =exp(-\frac{r^2{z}_R^2}{{\omega }_0^2}\frac{2z\Delta z+\Delta z^2}{(z_R^2+{(z+\Delta z)}^2)(z_R^2+z^2)}).$$

(9)

The bending of the wavefronts is another factor that worth mentioning during the field propagation. A numerical calculation has been carried out based on Equation (6). For typical FEL simulations, the electron beam radius is less than 100 $\mathsf{\mu }$m and the FEL Pierce parameter, denoted as $\rho$, is less than 0.001. Taking the SXFEL as an example, the radiation wavelength and the undulator period are $8.8$ $\mathrm{n}$m and $2.35$ $\mathrm{c}$m. respectively. The Rayleigh length is 15 m and the FEL gain length, denoted as $L_{gain}$, can be approximated by ${\lambda }_{\mu }/4\pi \sqrt{3}\rho$. The field phase difference between the central and edge of the electron beam at $z=L_{gain}$ is $1.9^{\circ }$. Normally, the electron beam radius is smaller than 100 $\mathsf{\mu }$m and the field phase difference will drastically decrease as the electron beam radius decreases. Various examples have been considered and it is safe to conclude that the field phase difference between the central and the edge of the electron beam after one FEL gain length is smaller than $3^{\circ }$. Therefore, practically, the field phase differences that are caused by the bending of the wavefronts can be ignored in this research.

In FALCON, the factors ${\chi }_r$ and ${\chi }_{amp}$ are involved in modifying the radiation field in every integration step. Such modifications can partially solve the diffraction effects in the FEL process. Figure 1 shows the FEL power gain curves with and without these diffraction modifications under the conditions of different electron beam energies and corresponding FEL wavelengths.

In this figure, the blue lines and red lines represent the FEL power with and without diffraction modifications, respectively. The black stars indicate the saturation points that are calculated by Xie’s formulas [29]. The errors of saturation power and saturation length are approximately $10\%$, showing significant improvements as compared to those without diffraction modifications.

2.5. The Working Equations

With the above assumptions, the particle equations are shown as [26]:

$$\begin{array}{cc}\hfill \frac{dx}{dz}& =p,\hfill \\ \hfill \frac{dp}{dz}& =-k_{\beta }^{2}x,\hfill \\ \hfill \frac{d\theta }{dz}& =k_u(1-\frac{{\gamma }_r^2}{{\gamma }^2})-\frac{k_s}{2}(p^2+k_{\beta }^2{x}^2),\hfill \\ \hfill \frac{d\gamma }{dz}& =\frac{e{a}_w\left[JJ\right]}{{\gamma }_{r}m{c}^2}\mathfrak{Re}\left(E(x,z)e^{i\theta }\right),\hfill \end{array}$$

(10)

where $x=(x,y)$ is the electron’s transverse position, $p=(p_x,p_y)$ is the electron’s transverse momenta, $k_u$ is the undulator wavenumber ($=2\pi /{\lambda }_{\mu }$), $k_s$ denotes the fundamental radiation wavenumber ($=2\pi /{\lambda }_s$), and the Bessel factor $\left[JJ\right]$ is equal to $J_0\left(\xi \right)-J_1\left(\xi \right)$ with $\xi =a_w^2/(2+2a_w^2)$.

Additionally, the radiation field equation that is based on SVEA can be derived from Ref. [30]:

$$(\frac{\partial }{\partial z}+\frac{{\nabla }_{\perp }^2}{2i{k}_s})E(x,z)=-\frac{{\mu }_0{c}^{2}e{a}_w\left[JJ\right]}{2}\sum _{j=1}^{N_e}\frac{e^{-i{\theta }_j}}{{\gamma }_j}\delta (x\mathbf{-}{x}_j,z-z_j).$$

(11)

In FALCON, the diffraction effect is modelled, as described in Section 2.4, and the diffraction term, ${\nabla }_{\perp }^2$, has been accordingly dropped to simplify the radiation field equation and increase the computation efficiency. Therefore, the radiation field equation is simplified as

$$\frac{d}{dz}E(x_0,z)=-\frac{{\mu }_{0}c}{2}\frac{I_0}{{\mathcal{A}}_{tr}}{a}_w\left[JJ\right]\langle \frac{e^{-i{\theta }_j}}{{\gamma }_j}\rangle ,$$

(12)

where $I_0$ denotes the beam current, ${\mathcal{A}}_{tr}$ is the beam cross section that equals to $2\pi {\sigma }_e^2$, and ${\mu }_0$ is the vacuum permeability. After each integration step, the amplitude of radiation field has been modified by the diffraction factor, ${\chi }_{amp}$. It is worth noting that there is no explicit dependence on t in the differential equation. The time-dependence effects are analyzed while using the ‘slice’ model.

2.6. Validation

Figure 2 shows a comparison of the simulation results based on the self-amplified spontaneous emission (SASE) scheme between FALCON and GENESIS.

Table 1. Simulation Parameters for the self-amplified spontaneous emission (SASE).

Parameters	Value	Unit
Beam energy	15	GeV
Normalized emittance	0.3	mm $·$ mrad
Relative energy spread	${10}^{-4}$	–
Peak current	4000	A
FEL wavelength	0.1	nm
Normalised undulator parameter	2.178	–
Undulator period	3	cm

shows the simulation parameters. The electron beam energy is 15 GeV with a relative slice energy spread of ${10}^{-4}$. The normalized emittance of the electron beam is about $0.3$ mm · mrad and the peak current is about 4000 A. The electron bunch length is approximately 60 $\mathrm{f}$s (FWHM). The undulator period is 3 cm with a normalised undulator parameter of 2.178. The FEL radiation wavelength is $0.1$ $\mathrm{n}$m.

The emitted pulse energies are about $1.83$ $\mathrm{m}$J in GENESIS and $1.47$ $\mathrm{m}$J in FALCON at the first saturation points. The difference in the power level of these codes seems apparent. The saturation length of FALCON is shorter than that of GENESIS, while the saturation powers are in line with each other. According to the assumption proposed in Section 2.4, a lack of mode competition process is probably responsible for the shortening of saturation length. Although FALCON is mainly designed for seeded FELs, the longitudinal properties for the SASE mode show good consistency between these two codes. Significantly, the simulation time is approximately 200 times longer in GENESIS 1.3 (more than 4 h) than in FALCON (less than 2 min) while using the same computer (Inter(R) Core(TM) i7-6700 [email protected]).

3. Results and Discussion

3.1. Example Simulations

The echo-enabled harmonic generation (EEHG) [27] is chosen as an example to show the typical simulation capabilities of FALCON. The example is not exhaustive at all. However, FALCON can be easily modified and extended to meet researchers’ various requirements.

Table 2. Simulation Parameters for the echo-enabled harmonic generation (EEHG).

Parameters	Value	Unit
Beam energy	840	MeV
Slice energy spread	50	keV
Normalized emittance	1.5	mm $·$ mrad
Peak current	800	A
Wavelength of seed laser	265	nm
Energy modulation amplitudes(A1/A2)	4/2	–
Dimensionless dispersion strength(B1/B2)	16.35/0.56	–
Normalised undulator parameter	1.012	–
Undulator period	23.5	cm
FEL wavelength	8.8	nm

lists the simulation parameters.

The electron beam energy is 840 MeV, with a slice energy spread of about 50 keV. The normalized emittance is $1.5$ mm · mrad and the peak current is 800 A. The FEL radiation wavelength is $8.8$ $\mathrm{n}$m, which is the 30th harmonic of the seed laser. The seed laser pulse that is used in the simulation is about 1 $\mathrm{p}$s (FWHM), which is much longer than the electron bunch length of about 320 $\mathrm{f}$s (FWHM). The energy modulation amplitudes in the modulator 1 and 2 are chosen to be $A_1=4$ and $A_2=2$, respectively. The optimized dispersion strengths for the 30th harmonic are found to be $B_1=16.35$ and $B_2=0.56$, corresponding to $R_{56}^{\left(1\right)}=11.58\mathrm{mm}$ and $R_{56}^{\left(2\right)}=0.40\mathrm{mm}$, respectively. The undulator period is $23.5$ cm with a normalised undulator parameter of 1.012.

The bunching factor at 30th harmonic of the seed is around 12% at the entrance of the radiator. Figure 3 shows the FEL power along radiator, FEL power profile, and spectrum at saturation. The pulse energies at saturation are $108.62$ $\mathsf{\mu }$J and $118.88$ $\mathsf{\mu }$J in GENESIS and FALCON, respectively.

For the seeded FEL operation mode, the seed laser should propagate co-axially with the electron beam in the modulator in order to ensure the formation of the desired microbunch. However, the drift of laser incident angle due to the mechanical vibration or the transverse position instability of the seed laser will bring about a small angle between the propagation direction of seed laser and electron beam. Such obliquely incident laser brings not only energy modulations, but also angular modulations [28]. Figure 4a,b demonstrate the transverse divergences ($\Delta y\prime$) and energy deviations ($\Delta \delta =\Delta \gamma /\gamma$) after the first modulator with incident angle of 1 $\mathrm{m}$rad in the $z-y$ plane. Figure 4c illustrates the influences of seed laser incident angle on the bunching factor and saturation power. For simplicity, the vertical laser angular drifts in both modulators are assumed to be the same. It can be seen from Figure 4c that the bunching factor at the entrance of the radiator and the saturation power decrease drastically as the incident angle increases.

3.2. Implementation at SXFEL

The primal target of this simulation code is to provide an assistant optimization tool at the Shanghai Soft X-ray FEL facility (SXFEL) [31], which aims at generating the soft X-ray radiation from a conventional seed laser through the frequency up-conversion scheme. An X-band RF deflector and a dipole magnet have been installed downstream of the undulator section, which provide a diagnostic tool for obtaining the image of time-energy longitudinal phase space of the electron beam [32]. In this section, the experiment and simulation results for the HGHG scheme at the SXFEL have been shown to demonstrate the feasibility of FALCON as an assistant online optimization tool.

The electron beam parameters can be acquired by FALCON through various diagnostic tools in real time. This enables FALCON to provide real-time predictions on the FEL performance and make comparisons between the measurements and simulations in order to analyze the experimental phenomena. The “lasing-off” longitudinal phase space images from the experiment are collected by turning off the seed laser. An example is shown in Figure 5, with the x-axis corresponding to energy and y-axis corresponding to time. The beam current (red line), central energy (black line), and energy spread (orange line) along the time axis can be derived from the image, as shown in the figure. These parameters combined with calibrated time can be used in order to generate realistic electron beam input file for FALCON.

Figure 6 shows one of the possible applications of FALCON. The longitudinal phase spaces on the screen after lasing through the experiment, GENESIS simulations, and FALCON simulations are presented in Figure 6a–c. The induced additional energy spread leaves a conspicuous imprint of the FEL radiation on the longitudinal phase space. Shouldering the main peak, there are two small sidebands that have been observed during the experiment and simulations. By analyzing the real-time simulation results, this part of electron beam developed very strong microbunching after lasing in the radiator, which was overbunched by the following chicane and emitted 266 nm coherent radiation in the downstream magnets. With these phase spaces, the longitudinal proprieties of FEL pulses can be determined by the method that was proposed in Refs. [33,34]. Figure 6d shows the reconstructed and simulated power profiles of FEL pulses. The pulse energies are about $16.22$ $\mathsf{\mu }$J, $16.58$ $\mathsf{\mu }$J and $17.71$ $\mathsf{\mu }$J in experiment, GENESIS, and FALCON, respectively.

It is notable that the way we use FALCON for on-line optimization is not to calculating every single shot. Instead, it plays a useful role in commissioning and operating the FEL machine under specific circumstances. FALCON can help researchers to simultaneously estimate some FEL properties (like pulse energy and spectrum) with the realistic parameters of the electron beam (import from the emittance and longitudinal phase space measurements) and give the correlation analysis. Moreover, FALCON can give some valuable information, even without the specific diagnostics tools. For example, the code can be applied in order to estimate the possible longitudinal phase space after lasing when a combination of the RF deflector and dipole magnet installed downstream the undulator is unavailable. Additionally, FALCON can help the researchers to understand the physics behind the experiment phenomena, such as the sidebands that are discussed above. These can give the operators and researchers more information about the FEL evolution process and the machine commissioning. The fairly accurate simulation results and greatly-improved computing efficiency have made it possible for FALCON to give guidance on the commissioning of the FEL.

4. Conclusions

This paper gives an introduction to a super-fast time-dependent simulation code, named FALCON, and its applications on a seeded FEL facility. A special care has been taken for the evolution of the radiation field envelope. The introduction of the radius factor (${\chi }_r$) and amplitude factor (${\chi }_{amp}$) has simplified the field equation and drastically reduced the required computational resources. The special designed modulation section removes the constraints on the laser-electron beam interaction with arbitrarily parameters. The present code is capable of simulating different FEL schemes within several minutes, which is about two orders of magnitude less than the time that is required by GENESIS 1.3. The validation of FALCON against GENESIS 1.3 has shown quite acceptable simulation results. FALCON shows a great potential as a fast optimization tool for FEL studies and machine operations due to the relatively high accuracy and computation efficiency.