SPECTRA: a synchrotron radiation calculation code

Tanaka, T.; Kitamura, H.

doi:10.1107/S090904950101425X

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 8| Part 6| November 2001| Pages 1221-1228

doi:10.1107/S090904950101425X

SPECTRA: a synchrotron radiation calculation code

Takashi Tanaka ^a ^* and Hideo Kitamura ^a

^aCoherent Synchrotron Light Source Physics Laboratory, The Institute of Physical and Chemical Research, Koto 1-1-1, Mikazuki-cho, Sayo-gun, Hyogo 679-5148, Japan
^*Correspondence e-mail: ztanaka@spring8.or.jp

(Received 3 July 2001; accepted 28 August 2001)

An application program to calculate various characteristics of synchrotron radiation, called SPECTRA, is described. The program does not need any other commercial software and is equipped with a full graphical user interface which makes data input quite easy. Equations on synchrotron radiation from arbitrary-field sources in a near-field region are derived, as are simplified expressions for ideal devices using a far-field approximation. Effective numerical methods implemented in SPECTRA to reduce computation time are explained, and several examples are presented.

Keywords: computer code; synchrotron radiation sources; numerical methods.

1. Introduction

In the construction of a beamline at a synchrotron radiation (SR) facility, the understanding of the light source is important in order to design various components to be installed in the beamline. The power distribution of SR is important for the design of the front-end components. A knowledge of the photon-beam emittance (size and divergence) is necessary for ray-tracing of optical elements such as monochromators and mirrors. Spectral evaluation up to very high energy ( $[\leq]$ 300 keV) regions is necessary for designing components for radiation shielding.

In addition to information on the light source necessary for construction of the beamline, other photon-beam properties should be estimated for an accurate analysis of experimental data recorded at the beamline. For example, the degree of polarization should be known in order to analyze experimental data obtained using a polarized photon beam, and a profile of the energy spectrum of the undulator radiation (UR) provides information on electron-beam performances such as the emittance and coupling constant.

We have developed a computer code called SPECTRA to calculate the SR properties mentioned above. It can calculate properties of SR emitted from devices with arbitrary magnetic fields. For wigglers and bending magnets, the well known SR expressions are used. For UR, the so-called far-field approximation can be used for fast computation. For more accurate evaluation of SR, expressions on SR in a near-field region are used for numerical computation. In this case, properties of SR emitted from both the ideal- and arbitrary-field devices can be calculated.

To compute SR properties, many parameters concerning the electron beam, SR source and sampling range on the photon energy and observation position should be specified. In SPECTRA, a graphical user interface is adopted to help the user to specify the many parameters to be used in the computation.

The purpose of this article is to inform the SR users of the details and usefulness of SPECTRA. In the following sections, analytical expressions on SR, numerical methods implemented in SPECTRA and examples of computation obtained for various sources are described.

2. Equations on SR

In this section, analytical expressions on SR are derived. All equations are written using SI units.

2.1. SR in the near-field region

The scalar and vector potentials generated by an electron moving along an arbitrary trajectory (Liénard–Wiechert potential) are given by (Landau & Lifshits, 1971 ; Jackson, 1975 )

$[\eqalign{&\varphi\,({\bf r},t)={e\over{4\pi\varepsilon_0[R-{\bf v}\cdot{\bf R}/c]_{\rm ret}}},\cr&{\bf A}({\bf r},t)={{\bf v}\over{4\pi\varepsilon_0c^2[R-{\bf v}\cdot{\bf R}/c]_{\rm ret}}},}]$

with

$[{\bf R}={\bf r}-{\bf r'},]$

where $[{\bf r}]$ $[({\bf r}')]$ is the vector directing from the origin to the observer (electron), $[{\bf v}]$ is the velocity of the electron, e is the electron charge, c is the speed of light and $[\varepsilon_0]$ is the permittivity of a vacuum. The subscript `ret' denotes that the values in the square bracket should be calculated at the retarded time $[t\,']$ determined by the equation

$[t\,'+R(t\,')/c=t.\eqno(1)]$

Using the relation between the potentials and the electric field,

$[{\bf E}=-\nabla\varphi-{{\partial{\bf A}}\over{\partial t}},\eqno(2)]$

we have

$[{\bf E}(t)={{e}\over{4\pi\varepsilon_0}}\left [{{1-\beta^{\,2}}\over{R^{\,2}}}{{{\bf n}-{\boldbeta}}\over{(1-{\bf n}\cdot{\boldbeta})^3}}+{{{\bf n}\times({\bf n}-{\boldbeta})\times{\dot{\boldbeta}}}\over{cR(1-{\bf n}\cdot{\boldbeta})^3}}\right]_{\rm ret},\eqno(3)]$

where

$[{\bf n}={\bf R}/R.]$

Using this expression, a power density, or radiated power (dP) per unit surface area (dS) is calculated as

$[\eqalignno{{{{\rm d}P}\over{{\rm d}S}}&=\varepsilon_o c\int_{-\infty}^{\phantom{.}\infty}|{\bf E}(t)|^2\,\,{\rm d}t\cr& ={{e^{\,2}\gamma^{\,4}}\over{2\pi^{\,2}\varepsilon_{\,0}^{\,2}c^{\,2}}}\int_{-\infty}^{\phantom{.}\infty}{{{\rm d}t}\over{R^{\,2}}} \cr&\quad\times\left[{{(\gamma\dot{\beta_x})^2+(\gamma\dot{\beta_x})^2}\over{D^3}} -{{4\,(\xi\gamma\dot{\beta_x}+\psi\gamma\dot{\beta_y})^2}\over{D^5}}+{{4c^{\,2}(\xi^{\,2}+\psi^{\,2})}\over{R^{\,2}D^5}}\right],&(4)}]$

with

$[\eqalign{D&=1+\xi^{\,2}+\psi^{\,2},\cr \xi&=\gamma{n_x}-\gamma\beta_x,\cr \psi&=\gamma{n_y}-\gamma\beta_y,}\eqno(5)]$

where $[\gamma]$ and $[\beta]$ are the Lorentz factor and the relative velocity of the moving electron, respectively.

In order to calculate the spectral flux density, or the number of photons (dN_p) in a relative spectral interval ( $[{\rm d}\omega/\omega]$ ) per unit surface area, the Fourier components of $[{\bf E}]$ should be calculated. It is more convenient to use the Fourier transformation of (2) instead of (3) as (Chubar, 1995 ; Chubar & Smolyakov, 1993 )

$[{\bf E}_{\omega}=-\nabla\varphi_{\omega}+i\omega{\bf A}_{\omega}.\eqno(6)]$

Using (1), we have

$[\varphi_{\omega}={{e}\over{4\pi\varepsilon_0}}\int_{-\infty}^{\phantom{.}\infty}{{1}\over{R(t\,')}}\exp\left[i\omega t(t\,')\right]\,{\rm d}t\,',\eqno(7)]$

$[{\bf A}_{\omega}={{e}\over{4\pi\varepsilon_0c}}\int_{-\infty}^{\phantom{.}\infty}{{{\boldbeta}\left(t\,'\right)}\over{R(t\,')}}\exp\left[i\omega t(t\,')\right]\,{\rm d}t.\eqno(8)]$

Substituting (7) and (8) into (6), we have

$[{\bf E}_{\omega}={e\over{4\pi\varepsilon_0c}}\,\,{\bf F}_{\omega},\eqno(9)]$

with

$[{\bf F}_{\omega}= i\omega\int_{-\infty}^{\phantom{.}\infty}{{1}\over{R(t\,')}}\left[{\boldbeta}(t\,')-\left(1+{{ic}\over{\omega R(t\,')}}\right){\bf n}(t\,')\right] \exp\left[i\omega t(t\,')\right]\,{\rm d}t\,'.\eqno(10)]$

The spectral flux density is calculated as

$[{{{\rm d}^2N_p}\over{{\rm d}S\,{\rm d}\omega/\omega}}={{{\rm d}^2 P}\over{\hbar\,{\rm d}S\,{\rm d}\omega}}= {{\varepsilon_oc}\over{\hbar\pi}}|{\bf E}_{\omega}^2|={{\alpha}\over{4\pi^2}}|{\bf F}_{\omega}|^2,\eqno(11)]$

where $[\alpha]$ is the fine structure constant. The Stokes parameters specifying the polarization properties of the radiation are calculated as

$[\eqalign{&s_0=|F_{\omega{x}}|^2+|F_{\omega{y}}|^2,\cr& s_1=|F_{\omega{x}}|^2-|F_{\omega{y}}|^2,\cr& s_2=2\,{\rm Re}(F_{\omega{x}}F_{\omega{y}}^{\,*}),\cr& s_3=2\,{\rm Im}(F_{\omega{x}}F_{\omega{y}}^{\,*}).}]$

Using (4) and (11), most of the radiation characteristics, such as the power distribution, energy spectra and degree of polarization, can be calculated for an electron moving in arbitrary magnetic fields.

Let the origin of the coordinates be at the midpoint of the device generating the magnetic fields. In this case, the distance between the electron and observer, R, is rewritten as

$[\eqalignno{R&=\left[(X-x)^2+(Y-y)^2+(Z-z)^2\right]^{1/2}\cr& \simeq Z-z+{{(X-x)^2+(Y-y)^2}\over{2(Z-z)}},&(12)}]$

with

$[{\bf r}=(X,Y,Z),\qquad{\bf r'}=(x,y,z),]$

where we have used approximations

$[|Z-z|\,\,\,\gt\!\!\gt\,\,|X-x|,|Y-y|,]$

being usually satisfied in general radiation processes.

The retarded time $[t\,']$ is expressed in terms of the electron longitudinal position z as

$[t\,'=\int_0^{\,z}{{{\rm d}z}\over{c\beta_z}}= {{1}\over{c}}\int_0^{\,z}\left[1+{{1+(\gamma\beta_x)^2+(\gamma\beta_y)^2}\over{2\gamma^{\,2}}}\right]\,{\rm d}z.\eqno(13)]$

The time t in (10) is therefore expressed in terms of z as

$[t={{1}\over{2c}}\left[\,{{z}\over{\gamma^{\,2}}}+\int_0^{\,z}\left(\beta_{\,x}^{\,2}+\beta_{\,y}^{\,2}\right)\,{\rm d}z+{{(X-x)^2+(Y-y)^2}\over{Z-z}}\right]+t_0,\eqno(14)]$

where t₀ is a constant independent of z and can be omitted when performing the integration in (10).

2.2. UR in the far-field region

Let us consider radiation emitted from an electron moving in an undulator with completely periodic magnetic fields (ideal undulator). The periodic length and number of periods are assumed to be $[\lambda_u]$ and N, respectively. If the distance between the electron and observer, R, is much longer than the length of the undulator, we can apply a far-field approximation to simplify the expressions on SR. In this case, $[{\bf n}]$ is regarded as constant and R is simplified to

$[R=r-{\bf n}\cdot{\bf r}'.]$

The above equation is used to derive a simplified expression of t as

$[\eqalign{t\,(z)={}&{{1}\over{2c}}\left\{\,{{z}\over{\gamma^{\,2}}}+\int_0^{\,z}\left[\beta_{\,x}^{\,2}(z)+\beta_{\,y}^{\,2}(z)\right]\,{\rm d}z-2x(z)\theta_x-2y(z)\theta_y+\theta^{\,2}z\right\}\cr&+{\rm constant},}]$

with

$[\theta_x=\tan^{-1}\left(X/Z\right),\quad\theta_y=\tan^{-1}\left(Y/Z\right),\quad\theta^{\,2}=\theta_x^{\,2}+\theta_y^{\,2}]$

being angles of observation. In the far-field approximation, R is regarded as a constant except when it is used to calculate t (z).

Because the magnetic field is completely periodic, the velocity and position of the moving electron is also periodic. The time t therefore satisfies the condition

$[t\,(z+\lambda_u)=t\,(z)+T,\eqno(15)]$

with

$[T=2\pi/\omega_1,]$

$[\omega_1={{4\pi{c}\gamma^{\,2}}\over{\lambda_u}}{{1}\over{1+\left\langle(\gamma\beta)^2\right\rangle+(\gamma\theta)^2}},\eqno(16)]$

$[\left\langle\left(\gamma\beta\right)^2\right\rangle={{1}\over{\lambda_u}}\int_0^{\,\lambda_u}\left[\left(\gamma\beta_x\right)^2+\left(\gamma\beta_y\right)^2\right]\,{\rm d}z.\eqno(17)]$

Assuming that R is much longer than the wavelength of radiation $[2\pi c/\omega]$ , $[{\bf F}_\omega]$ is simplified to

$[{\bf F}_\omega= {{1}\over{R}}\int_0^{\,T\,'}\left[\,{\boldbeta}\left(t\,'\right)-{\bf n}\,\right]\,H\left[t\,(t\,')\right]\,{\rm d}t\,',\eqno(18)]$

with

$[H(t)=i\omega\sum_{n\,=\,0}^{N-1}\exp[i\omega(t+NT)]={{\rm d}\over{{\rm d}t}}\sum_{n\,=\,0}^{N-1}\exp[i\omega(t+NT)],]$

where $[T\,']$ is the period in the rest frame of the electron.

Expanding the function H(t) in a Fourier series between t (0) = 0 and $[t\,(T\,')]$ = T, we have

$[H(t)=\sum_{k\,=\,-\infty}^{\phantom{.}\infty}h_k\exp\left(ik\omega_1t\right),\eqno(19)]$

with

$[h_k={{1}\over{T}}\int_0^T\!\!H(t)\exp\left(-ik\omega_1t\right)\,{\rm d}t={{ik\omega_1}\over{T}}{{\exp\left[i\left(\omega-k\omega_1\right)NT\right]-1}\over{i(\omega-k\omega_1)}},\eqno(20)]$

where we have made a partial integration to derive (20).

Substituting (19) and (20) into (18), we have

$[{\bf F}_{\omega}=\sum_{k\,=\,-\infty}^{\phantom{.}\infty}{{\exp\left[2\pi{i}\left(\omega/\omega_1-k\right)N\right]-1}\over{i(\omega/\omega_1-k)T}}\,{\bf F}_k,]$

with

$[{\bf F}_k={{ik}\over{RT}}\int_0^{\,T\,'}\left[\,{\boldbeta(t\,')}-{\bf n}\,\right]\exp\left[ik\omega_1t\,(t\,')\right]\,{\rm d}t\,'.]$

For convenience, we change the variable of integration from $[t\,']$ to z to obtain

$[{\bf F}_k={{ik}\over{cRT}}\int_0^{\,\lambda_u}\left[\,{\boldbeta(z)}-{\bf n}\,\right]\exp\left[ik\omega_1t(z)\right]\,{\rm d}z.\eqno(21)]$

In order to calculate the flux density and Stokes parameters, products of components of vector $[{\bf F}_{\omega}]$ should be calculated. For example, the product $[F_{\omega x}\,F_{\omega x}]$ is calculated as

$[\eqalignno{F_{\omega{x}}\,F_{\omega{x}}={}&\sum_{k\,=\,-\infty}^{\phantom{.}\infty}F_{kx}\,F_{kx}\left[{{\sin\pi{N}(\omega/\omega_1-k)}\over{\pi(\omega/\omega_1-k)}}\right]^2\cr& +\sum_{k\,=\,-\infty}^{\phantom{.}\infty}\sum_{j\,\neq\,{k}}F_{kx}\,F_{jx}\,{{\sin\pi{N}(\omega/\omega_1-k)\sin\pi{N}(\omega/\omega_1-j\,)}\over{\pi^2(\omega/\omega_1-k)(\omega/\omega_1-j\,)}}.\cr&&(22)}]$

Other products are calculated in the same manner.

3. Electron-beam convolution

All equations described in the previous section concern the radiation emitted by a single electron, or a monochromatic electron beam with zero emittance. In order to evaluate practical properties of SR, the effects caused by an energy spread $[\sigma_E/E]$ and finite emittance $[\epsilon]$ of the electron beam should be taken into account.

3.1. Phase ellipse transfer

The effect due to the finite electron-beam emittance and energy spread is described as a beam envelope at the position of observation. The electron, after emitting radiation, is in general deflected by bending magnets. It is important, however, to know the beam envelope at the position of observation if the electron beam travels without any magnetic field, or in a drift space. Assuming the distance between the light source and observer to be l, the beam envelope $[\sigma]$ at the observation position is calculated as

$[\sigma=\left[\epsilon\left(\beta_0-2\alpha_0l+{{1+\alpha_0^2}\over{\beta_0}}l\right)+\left({{\sigma_E}\over{E}}\right)^2\left(\eta_0^2+\eta_{\,0}^{\,\prime\,2}l^{\,2}\right)\right]^{1/2}, \eqno(23)]$

with

$[\alpha_0=-\beta_0^{\,\prime}/2,]$

where $[\beta_0]$ , $[\eta_0]$ , $[\beta_0^{\,\prime}]$ and $[\eta_0^{\,\prime}]$ are the betatron and dispersion functions and their derivatives at the position of the light source, respectively. Both horizontal and vertical beam envelopes are calcuated using (23). The finite beam emittance is taken into account by a convolution of the result for a single electron with a two-dimensional electron distribution function.

3.2. Energy spread

The effect due to the energy spread of the electron beam has already been described in the previous section. There is another important effect for the UR calculation: because the energy of the first harmonic is proportional to $[\gamma^{\,2}]$ , the energy spread of the electron beam causes a spectral broadening. Its effect is therefore taken into account by an energy convolution with a Gaussian distribution function with an root mean square (r.m.s.) width of 2 $[\hbar\omega\sigma_E/E]$ , where $[\hbar\omega]$ is the photon energy. The effect is pronounced especially for UR energy spectra when the number of periods is very large, e.g. for those of the 25 m in-vacuum undulator installed at SPring-8, and/or when the spectral energy range includes very high harmonics.

3.3. General case

Let us consider a function F as functions of X and Y, or the transverse positions of observation. Assuming the electron distribution function to be Gaussian, convolution of F with the electron beam, F_c, is expressed as

$[\eqalignno{F_c(X,Y)={}&{{1}\over{2\pi\sigma_x\sigma_y}}\int_{-\infty}^{\phantom{.}\infty}\,{\rm d}x\int_{-\infty}^{\phantom{.}\infty}\,{\rm d}y\cr&\quad\times{F(x,y)}\exp\left[-{{(x-X)^2}\over{2\sigma_x^2}}-{{(y-Y)^2}\over{2\sigma_y^2}}\right],&(24)}]$

where $[\sigma_{x,y}]$ is the horizontal and vertical r.m.s. beam envelope calculated using (23).

In order to obtain a quantity passing though a slit $[\Delta F_c]$ , F_c(X,Y) should be integrated over the slit area S,

$[\eqalignno{\Delta{F_c}=&\int_SF_c(X,Y)\,{\rm d}S={{1}\over{2\pi\sigma_x\sigma_y}}\int_{-\infty}^{\phantom{.}\infty}\,{\rm d}x\int_{-\infty}^{\phantom{.}\infty}\,{\rm d}y\cr&\times{F(x,y)}\int_S\,{\rm d}S\exp\left[{-{{(x-X)^2}\over{2\sigma_x^2}}-{{(y-Y)^2}\over{2\sigma_y^2}}}\right].&(25)}]$

If the slit has a rectangular shape, the above equation is simplified to

$[\eqalign{\Delta{F_c}={}& {{1}\over{4}}\int_{-\infty}^{\phantom{.}\infty}\,{\rm d}x\int_{-\infty}^{\phantom{.}\infty}\,{\rm d}y\,F(x,y)\cr&\times\left[{\rm erf}\left({{X_2-x}\over{2^{1/2}\sigma_x}}\right)-{\rm erf}\left({{X_1-x}\over{2^{1/2}\sigma_x}}\right) \right]\cr&\times\left[{\rm erf}\left({{Y_2-y}\over{2^{1/2}\sigma_y}}\right)-{\rm erf}\left({{Y_1-y}\over{2^{1/2}\sigma_y}}\right)\right],}]$

with

$[X_2=X_c+\Delta{X/2},\quad\quad{X_1}=X_c-\Delta{X/2},]$

$[Y_2=Y_c+\Delta{Y/2},\quad\quad{Y_1}=Y_c-\Delta{Y/2},]$

where the slit has a horizontal and vertical aperture $[\Delta{x}]$ and $[\Delta{y}]$ located at (X_c,Y_c), and erf(x) is an error function defined as

$[{\rm erf}(x)={{2}\over{\pi^{1/2}}}\int_0^{\phantom{.}\infty}\exp\left(-x^2\right)\,{\rm d}x.]$

SPECTRA is applicable to slits with rectangular, circular or doughnut-shaped apertures. All the radiation properties are calculated with the schemes described above. For an energy spectrum of UR in the far-field region, however, the scheme described in the following section is applied to speed up the computation.

3.4. UR spectra in the far-field region

If N is much larger than unity, the second term in (22) can be neglected and the radiation can be expanded into harmonics with multiple energies of $[\hbar\omega_1]$ as

$[{{{\rm d}^2N_p(X,Y)}\over{{\rm d}S\,{\rm d}\omega/\omega}}=\sum_{k\,=\,-\infty}^{\phantom{.}\infty}f_k(X,Y)\left[{{\sin{N}\pi(\omega/\omega_1-k)}\over{N\pi(\omega/\omega_1-k)}}\right]^2,\eqno(26)]$

with

$[f_k(X,Y)={{\alpha}\over{4\pi^2}}|{\bf F}_k|^2.]$

Introducing a delta function, we can modify the above equation as

$[{{{\rm d}^2N_p(X,Y)}\over{{\rm d}S\,{\rm d}\omega/\omega}}=N^2\!\!\sum_{k\,=\,-\infty}^{\phantom{.}\infty}f_k(X,Y)\int_{-\infty}^{\phantom{.}\infty}\!\!\delta(\omega'-k\omega_1)P_N(\omega,\omega_1,\omega')\,{\rm d}\omega',\eqno(27)]$

with

$[P_N\left(\omega,\omega_1,\omega'\right)=\left[{{\sin N\pi(\omega-\omega')/\omega_1}\over{N\pi(\omega-\omega')/\omega_1}}\right]^2.\eqno(28)]$

The convolution of $[{\rm d}^2N_p/({\rm d}S\,{\rm d}\omega/\omega)]$ with the electron beam is expressed as

$[\eqalignno{{{{\rm d}^2N_p}\over{{\rm d}S\,{\rm d}\omega/\omega}}={}&{{1}\over{2\pi\sigma_x\sigma_y}}\int_{0}^{\phantom{.}\infty}\rho\,{\rm d}\rho\int_{0}^{\,2\pi}\,{\rm d}\varphi{{{\rm d}^2N(\rho\cos\varphi,\rho\sin\varphi)}\over{{\rm d}S\,{\rm d}\omega/\omega}}\cr&\times\exp\left[-{{(X-\rho\cos\varphi)^2}\over{2\sigma_x^2}}-{{(Y-\rho\sin\varphi)^2}\over{2\sigma_y^2}}\right].&(29)}]$

By using (27), (29) is simplified to

$[\eqalignno{{{{\rm d}^2N_p}\over{{\rm d}S\,{\rm d}\omega/\omega}}={}&{{N^2Z^2}\over{4\pi\sigma_x\sigma_y}}\int_{-\infty}^{\phantom{.}\infty}{{{\rm d}\omega'}\over{\omega'}}{{4\pi\gamma^{\,2}c/\lambda_u}\over{\omega'}}\sum_{k\,\gt\,\omega'/\omega_1}^{\phantom{.}\infty}kP_N(\omega,\omega'/k,\omega')\cr&\times\int_{0}^{\,2\pi}{\rm d}\varphi\,\,f_k\left[\rho_{k}(\omega')\cos\varphi,\rho_{k}(\omega')\sin\varphi\right]G_k(X,Y,\varphi,\omega'),\cr&&(30)}]$

with

$[G_k(X,Y,\varphi,\omega)= \exp\left\{{-{{\left[X-\rho_{k}(\omega)\cos\varphi\right]^2}\over{2\sigma_x^2}}-{{\left[Y-\rho_{k}(\omega)\sin\varphi\right]^2}\over{2\sigma_y^2}}}\right\},]$

$[\rho_k(\omega)= Z\left[{{{4\pi\gamma^{\,2}ck}\over{\omega\lambda_u}}-\left\langle(\gamma\beta)^2\right\rangle-1}\right]^{1/2}.]$

The integration along $[\rho]$ has been replaced by that along $[\omega']$ .

The function P_N in the integrand of the above equation has a sharp peak at $[\omega]$ = $[\omega']$ and a width proportional to $[\omega'/k]$ . In other words, the width is dependent on the integration variable $[\omega']$ . In SPECTRA, a replacement of the function P_N, i.e.

$[P_N(\omega,\omega'/k,\omega')\,\,\rightarrow\,\,\omega'/(k\omega_1)P_N(\omega,\omega_1,\omega'),]$

is made for fast computation. The coefficient $[\omega'/(k\omega_1)]$ is multiplied to avoid a discrepancy in area. Equation (30) therefore reduces to

$[{{{\rm d}^2N_{p}}\over{{\rm d}S\,{\rm d}\omega/\omega}}=\int_{-\infty}^{\phantom{.}\infty}\left[{{{\rm d}^2N_{p\infty}}\over{{\rm d}S\,{\rm d}\omega/\omega}}\right]_{\omega\,=\,\omega'}NS_N(\omega-\omega'){{{\rm d}\omega'}\over{\omega_1}},\eqno(31)]$

with

$[\eqalign{{{{\rm d}^2N_{p\infty}}\over{{\rm d}S\,{\rm d}\omega/\omega}}={}&{{N}\over{4\pi\sigma_x\sigma_y}}{{4\pi\gamma^{\,2}c/\lambda_u}\over{\omega}}\sum_{k\,\gt\,\omega/\omega_1}^{\phantom{.}\infty}\int_{0}^{\,2\pi}{\rm d}\varphi\cr&\times{f_k}\left[\rho_{k}(\omega)\cos\varphi,\rho_{k}(\omega)\sin\varphi\right]G_k(X,Y,\varphi,\omega),}]$

$[S_N(x)=\left({{\sin{N\pi}x}\over{{N\pi}x}}\right)^2.]$

As is easily understood, $[{\rm d}^2N_{p\infty}/({\rm d}S\,{\rm d}\omega/\omega)]$ is the flux density for an infinite number of N. Using this scheme, integration along $[\omega']$ in (30) reduces to a convolution with a function S_N, which enables us to apply a convolution scheme with a fast Fourier transform algorithm.

4. Numerical method

In SPECTRA, several numerical methods are implemented to save computation time, which are briefly explained in this section.

4.1. Fast computation of F_ω

Let us consider an integration of the form

$[G=\int{g(t)}\exp\left(i\omega{t}\right)\,{\rm d}t.\eqno(32)]$

Both $[{\bf F}_\omega]$ and $[{\bf F}_k]$ can be reduced to the above form by changing the integration variable from z to t. When $[\omega{t}]$ is much larger than unity, the integrand oscillates rapidly and should be computed many times for an accurate integration. By integrating (32) by parts n times, we have

$[G=\left(-{{1}\over{i\omega}}\right)^{\!n}\int{{{\rm d}^{\,n}g}\over{{\rm d}t^{\,n}}}\exp\left(i\omega{t}\right)\,{\rm d}t\,\,-\,\,\sum_{k\,=\,1}^{n}\left(-{{1}\over{i\omega}}\right)^{\!k}\,\left[{{{\rm d}^{\,k-1}g}\over{{\rm d}t^{\,k-1}}}\exp\left(i\omega{t}\right)\right].]$

Let us consider the integration over the range between t₁ and t₂. If g(t) can be approximated by an nth-order polynomial in this range, dⁿg/dtⁿ becomes a constant; therefore the first term in the above equation reduces to

$[\left(-{{1}\over{i\omega}}\right)^{\!\!n+1}\,\,{{{\rm d}^{\,n}g}\over{{\rm d}t^{\,n}}}\left[\exp\left(i\omega{t_1}\right)-\exp\left(i\omega{t_2}\right)\right],]$

which means that the integration reduces to a simple summation. It is expected that the time taken to compute $[F_{\omega}]$ based on the above expression is independent of $[\omega]$ or the photon energy, which is a great advantage for high-energy regions.

In SPECTRA, the whole region of integration is divided into several sections within which g(t) is well approximated by a third-order polynomial. The derivatives up to the third order are determined by a cubic-spline interpolation method (Stoer & Bulirsch, 1991 ). The integration G is therefore reduced to

$[\eqalign{G={}&\sum_{k\,=\,1}^{4}\left(-{{1}\over{i\omega}}\right)^{\!k}\,\left[g^{(k-1)}(t_1)\exp\left(i\omega{t_1}\right)-g^{(k-1)}(t_{m+1})\exp\left(i\omega{t_{m+1}}\right)\right]\cr&+{{1}\over{\omega^4}}\sum_{j\,=\,1}^{m-1}\left[g_{j+1}^{(3)}-g_{j}^{(3)}\right]\exp\left(i\omega{t_{j+1}}\right),}]$

where t₁ and t_m+1 show boundaries of the whole integration range, g^(k) is the kth-order derivative of g, and g⁽³⁾_j shows the third-order derivative of g between t_j and t_j+1.

4.2. Facile evaluation of the total flux

By integrating (11) over the whole solid angle, the total flux or the total number of photons emitted per unit time in a relative spectral interval is obtained. This quantity is frequently used to design radiation shielding of an experimental hutch. In this case, the spectral data up to an extremely high energy is necessary, which imposes calculation of quite high harmonic radiation such as $[k\simeq1000]$ or higher, resulting in much computation time. It should be noted, however, that a sharp peak specific to UR is no longer observed in such a high-energy region and the spectrum is similar to that of bending-magnet radiation (BR). Therefore, it is expected that evaluation of the total flux based on an expression on BR will be a good approximation of that of UR. The linear flux density of BR, a quantity obtained by integrating the flux density over Y, or the vertical position of observation, is analytically expressed as

$[{{{\rm d}^2N_{pBM}}\over{{\rm d}X\,{\rm d}\omega/\omega}}={3^{1/2}\over{2\pi}}{{\alpha\gamma}\over{Z}}{{\omega}\over{\omega_c}}\int_{\omega/\omega_c}^{\phantom{.}\infty}K_{5/3}(\eta)\,{\rm d}\eta,\eqno(33)]$

where $[\omega_c]$ is a critical frequency defined as

$[\omega_c=3\gamma^{\,3}c/(2\rho),\eqno(34)]$

with

$[\rho=\gamma{mc}/(eB),\eqno(35)]$

where $[\rho]$ is a bending radius and B is the uniform magnetic field of the bending magnet.

If the magnetic field is not uniform, then (33) cannot be applied. We have to divide the whole length into many regions and regard the magnetic field as constant in each region. The flux emitted by an electron passing through the ith region is expressed as

$[{{{\rm d}N_{pi}}\over{{\rm d}\omega/\omega}}=\left[{{{\rm d}^2N_{BM}}\over{{\rm d}X{\rm d}\omega/\omega}}\right]_{B\,=\,B_i}\Delta\beta_iZ,]$

where B_i and $[\Delta\beta_i]$ are the magnetic field and change of relative velocity in the ith region. Summation of $[{\rm d}N_{pi}/({\rm d}\omega/\omega)]$ reduces to the total flux. The factor $[\Delta\beta_iZ]$ corresponds to a length of the region over which the radiation from the ith bending magnet is emitted.

It is easy to show that the above equation can be expressed in an integral form as

$[{{{\rm d}N_p}\over{{\rm d}\omega/\omega}}=Z\int\left[{{{\rm d}^2N_{BM}}\over{{\rm d}X\,{\rm d}\omega/\omega}}\right]_{B\,=\,B(z)}\left|{{{\rm d}{\boldbeta}}\over{{\rm d}z}}\right|\,{\rm d}z,\eqno(36)]$

with

$[B(z)=\left[B_y^2(z)+B_x^2(z)\right]^{1/2}.]$

Substituting (33) and considering the reaction between $[{\boldbeta}]$ and $[{\bf B}]$ , we have

$[{{{\rm d}N_p}\over{{\rm d}\omega/\omega}}={3^{1/2}\over{2\pi}}{{\alpha{e}}\over{mc}}\int{\rm d}z\,{{\omega}\over{\omega_c(z)}}\int^{\phantom{.}\infty}_{\omega/\omega_c(z)}K_{5/3}(\eta)\,{\rm d}\eta,\eqno(37)]$

where $[\omega_c(z)]$ is the critical frequency calculated with B(z).

4.3. Omission of radial integration

In order to obtain a precise spectrum of UR with finite beam emittance, two-dimensional integration along an azimuthal ( $[\varphi]$ ) and radial ( $[\rho]$ ) direction is necessary. As shown in §3.3, the radial integration can be reduced to an energy convolution. For a precise calculation, $[{\rm d}^2N_{p\infty}/({\rm d}S\,{\rm d}\omega/\omega)]$ should be calculated with fine steps of $[\omega]$ . In an energy region of the spectrum composed of high harmonics, however, the convolution with the S_N function changes little the shape of the spectrum. In such a case, the convolution, or the radial integration, can be omitted. In SPECTRA, the energy region to omit the radial integration can be specified to save computation time without sacrificing the accuracy of the calculation.

5. Program details

5.1. Types of SR sources

In SPECTRA, SR sources can be classified into four types, i.e. bending magnets, wigglers, undulators and arbitrary-field devices. The numerical method to be applied depends on the type of the SR source.

If the SR source is regarded as a bending magnet, the well known expressions on SR (Schwinger, 1949 ; Kim, 1986 ) are applied. In the case of wigglers, the magnetic fields are varied according to the horizontal observation angle. The undulator includes familiar devices, such as planar, helical, elliptical and figure-8 undulators. In addition, a device with arbitrary but periodic magnetic fields is also regarded as an undulator.

In fact, the magnetic field of a practical device is not completely periodic. In such a case, SPECTRA calculates SR properties in the near-field region without any approximation and assumption. Using this scheme, degradation of peak intensity at each harmonic of UR can be calculated. Edge radiation, i.e. SR emitted from fringing fields of bending magnets, can also be characterized. Evaluation of edge radiation is important for estimating performances of an X-ray beam-position monitor installed in the undulator beamline because it causes a significant background.

Radiation properties can be calculated as functions of photon energy and observation position. In the case of undulators, the K-value (deflection parameter) dependence can also be obtained.

5.2. Graphical user interface

In addition to the numerical methods implemented in a computer code, the user interface which helps the user to input their desired parameters is very important. SPECTRA adopts a graphical user interface (GUI) with which the user can specify computation parameters, as shown in Fig. 1. Storage-ring, SR-source and sampling parameters can be edited by the GUI. The storage-ring parameters include the electron energy, average current, emittance, twiss parameters and so forth. The SR-source parameters include a strength and periodic length of the magnetic field, or the field profile data for arbitrary-field devices. The sampling parameters include a computation range of the photon energy or observation positions, harmonic numbers of UR and so forth.

Figure 1
Graphical user interface as a pre-processor of SPECTRA.

For operating-system portability, a Tcl/Tk programming language has been adopted for building the GUI. Because the main code for numerical computation is written in a simple C language, SPECTRA can be executed on almost all of the platforms which Tcl/Tk supports, such as MS Windows, Macintosh, Linux and other kinds of UNIX.

6. Examples

In this section, examples of computation of various properties of SR are presented to show the usefulness of SPECTRA. The storage-ring parameters are shown in Table 1.

Table 1
Storage-ring parameters used in the calculation

Electron energy	8 GeV
Natural emittance	6 × 10⁻⁹ m rad
Coupling constant	0.3%
Horizontal betatron value	24 m
Vertical betatron value	6 m
Energy spread	0.1%

6.1. Energy spectrum

First of all, energy-spectrum computations are shown. As an example, a planar undulator with a periodic length of 32 mm and number of periods 140 is considered.

Fig. 2 shows spectra of the total flux, on-axis flux density and partial flux passing through a slit with dimensions of 1 mm × 1 mm located 30 m from the source. The far-field approximations are used in this example. The K value is assumed to be 2.56. The dotted line is the total flux spectrum evaluated by the scheme described in §4.2. The difference between the two total-flux spectra is small above 10 keV, indicating the validity of the scheme.

Figure 2
Energy spectra of the on-axis flux density, partial flux passing through a 1 mm × 1 mm slit located 30 m from the source, and total flux. The dotted line is the total flux calculated with the scheme described in §4.2

Fig. 3 shows the near-field effect on spectra of UR. The K value is assumed to be 1.34 for the first-harmonic energy to be at 10 keV. Fig. 3(a) shows on-axis spectra of the angular flux density, i.e. the flux density normalized by the square of the distance between the light source and observer, while Fig. 3(b) shows off-axis spectra observed at $[\theta_x]$ = 0 and $[\theta_y =]$ 60 µrad. Comparing the two figures, the near-field effect is found to be more remarkable for the off-axis spectra, as is described by Hirai et al. (1984 ) and Walker (1988 ).

Figure 3
Energy spectra observed (a) on axis and (b) at $[\theta_y]$ = 60 µrad calculated with the near-field scheme for various distances between the light source and observer.

Fig. 4 shows effects due to magnetic field errors, which are inevitable for a practical device. Fig. 4(a) shows electron orbits in both transverse directions assuming a permanent-magnet device composed of magnet blocks whose r.m.s. magnetization error is 1% in magnitude and 1° in angle. An on-axis spectrum calculated with the orbit data is shown in Fig. 4(b) together with an ideal case (without error). The peak, corresponding to harmonics up to the 11th, are found in each case. The difference is larger for higher harmonics. For example, the 11th-harmonic (= 49 keV) peak intensity for the error-field case is about half of that of the ideal.

Figure 4
Effects due to magnetization errors of undulator magnets. The electron trajectories are shown in (a) and the energy spectra with and without magnetization errors are shown in (b).

6.2. Spatial distribution

Next, spatial-distribution computations are shown. As an example, the so-called edge radiation from fringe fields of bending magnets is considered. Fig. 5 shows a typical magnetic field distribution along the longitudinal direction generated by two bending magnets. The flux density calculated with the field distribution is shown as functions of the horizontal and vertical observation positions in Figs. 6(a) and 6(b), respectively. The photon energy and the distance between the observer and midpoint of the straight section are assumed to be 5 eV and 20 m, respectively. The dotted line shows a zero-emittance case. A clear interference pattern is seen in each direction, while the effect due to the finite emittance is more pronounced for the horizontal case. The power density, also shown in Figs. 6(a) and 6(b), is found to have an asymmetry in the horizontal direction due to the difference in distance from the two bending magnets to the observer.

Figure 5
A fringe magnetic field distribution generated by two bending magnets.

Figure 6
Edge radiation spatial profiles along (a) horizontal and (b) vertical axes.

6.3. Filtering effect

In SR beamlines, several kinds of filters are used for various purposes. SPECTRA can calculate an energy spectrum and power after passing through a filter made of various materials (filtering). As an example, let us consider a filter made of copper with a thickness of 10 µm. Fig. 7 shows on-axis power densities before and after filtering as functions of the first-harmonic energy when the magnetic field is varied in order to tune the photon energy. The light source is the same as that described in §6.1. The dip found at 8.9 keV is due to the K-absorption edge of Cu. Figs. 8(a) and 8(b) show power densities as functions of the horizontal and vertical observation positions when the first-harmonic energy is set at 10 keV. The difference of profile within the range $[|\theta_{x,y}|]$ < 0.03 mrad is also caused by the Cu K-edge. The first-harmonic energy observed at $[|\theta_{x,y}|]$ = 0.03 mrad is equal to 8.9 keV in this case.

Figure 7
On-axis power densities before and after the 10 µm Cu filter as a function of the first-harmonic energy.

Figure 8
Power density spatial profiles along (a) horizontal and (b) vertical axes before and after the Cu filter.

In SPECTRA, filtering material frequently used in SR beamlines, such as aluminium, carbon and beryllium, are implemented. Other materials are also available by creating a new parameter set specifying its atomic number, density and mass ratio.

7. Summary

We have described the details of SPECTRA and shown several examples of computation. It is freely available from the SPECTRA homepage (https://radiant.harima.riken.go.jp/spectra/index_e.html ) on the SPring-8 web site. It does not require any other commercial software or libraries. In other words, it is a stand-alone application program.

Besides several numerical methods described in §4, the expansion of $[{\bf F}_k]$ by Bessel functions (Motz, 1951 ; Kincaid, 1977 ; Yamamoto & Kitamura, 1987 ) are implemented for a planar, helical and elliptical undulator in the far-field region. Spatial symmetries of radiation from these devices are also considered, while no symmetries are assumed for arbitrary-field devices because they are not ensured generally.

The computation time depends on the type of the computer and calculation. For example, the partial-flux spectrum calculation shown in Fig. 2 is performed with 10000 points of photon energy in 17 s using a computer with a Pentium 750 MHz processor. The calculation of the on-axis spectrum observed at the point 30 m from the source shown in Fig. 3 takes about 2 min with 5000 points of photon energy. The spatial dependence of the power density after filtering, as shown in Fig. 8, is calculated in 11 s with 100 points of observation position.

Finally, let us make a short comparison with other programs. According to our experiences, the computation time of SPECTRA is almost equal to those of other programs such as SRW (Chubar & Elleaume, 1998 ). It should be noted, however, that SPECTRA is more user-friendly, especially in the near-field computation including the electron-beam effects. In the case of SRW, the user must try many sets of parameters such as grid points and spatial range and use post-processing to obtain a reasonable result, while SPECTRA does not require such trial and error.

References

Chubar, O. V. (1995). Rev. Sci. Instrum. 66, 1872–1874. CrossRef CAS Web of Science
Chubar, O. V. & Elleaume, P. (1998). Proceedings of the 1998 European Particle Accelerator Conference, Stockholm, Sweden, pp. 1177–1179.
Chubar, O. V. & Smolyakov, N. V. (1993). Proceedings of the 1993 Particle Accelerator Conference, Washington DC, USA, pp. 1626–1628.
Hirai, Y., Luccio, A. & Yu, L. (1984). J. Appl. Phys. 55, 25–32. CrossRef Web of Science
Jackson, J. D. (1975). Classical Electrodynamics. New York: Wiley.
Kim, K. J. (1986). Proceedings of the 1986 US Particle Accelerator Summer School. Berlin: Springer.
Kincaid, B. M. (1977). J. Appl. Phys. 48, 2684–2691. CrossRef Web of Science
Landau, L. D. & Lifshits, E. M. (1971). The Classical Theory of Fields. Oxford: Pergamon.
Motz, H. (1951). J. Appl. Phys. 22, 527–535. CrossRef Web of Science
Schwinger, J. (1949). Phys. Rev. 75, 1912–1925. CrossRef Web of Science
Stoer, J. & Bulirsch, R. (1991). Introduction to Numerical Analysis. New York: Springer-Verlag.
Walker, R. P. (1988). Nucl. Instrum. Methods, A267, 537–546. CrossRef CAS
Yamamoto, S. & Kitamura, H. (1987). Jpn. J. Appl. Phys. 26, L1613–L1615. CrossRef

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 8| Part 6| November 2001| Pages 1221-1228

doi:10.1107/S090904950101425X