Raycing backend¶
Package raycing
provides the internal backend of xrt. It
defines beam sources in the module sources
,
rectangular and round apertures in apertures
,
optical elements in oes
, material properties
(essentially reflectivity, transmittivity and absorption coefficient) for
interfaces and crystals in materials
and screens
in screens
.
Coordinate systems¶
The following coordinate systems are considered (always righthanded):
The global coordinate system. It is arbitrary (userdefined) with one requirement driven by code simplification: Zaxis is vertical. For example, the system origin of Alba synchrotron is in the center of the ring at the ground level with Yaxis northward, Z upright and the units in mm.
Note
The positions of all optical elements, sources, screens etc. are given in the global coordinate system. This feature simplifies the beamline alignment when 3D CAD models are available.
The local systems.
of the beamline. The local Y direction (the direction of the source) is determined by azimuth parameter of
BeamLine
– the angle measured cw from the global Y axis. The local beamline Z is also vertical and upward. The local beamline X is to the right. At azimuth = 0 the global system and the local beamline system are parallel to each other. In most of the supplied examples the global system and the local beamline system coincide.of an optical element. The origin is on the optical surface. Z is outofsurface. At pitch, roll and yaw all zeros the local oe system and the local beamline system are parallel to each other.
Note
Pitch, roll and yaw rotations (correspondingly: Rx, Ry and Rz) are defined relative to the local axes of the optical element. The local axes rotate together with the optical element!
Note
The rotations are done in the following default sequence: yaw, roll, pitch. It can be changed by the user for any particular optical element. Sometimes it is necessary to define misalignment angles in addition to the positional angles. Because rotations do not commute, an extra set of angles may become unavoidable, which are applied after the positional rotations. See
OE
.The usersupplied functions for the surface height (z) and the normal as functions of (x, y) are defined in the local oe system.
of other beamline elements: sources, apertures, screens. Z is upward and Y is along the beam line. The origin is given by the user. Usually it is on the original beam line.
xrt sequentially transforms beams (instances of
Beam
) – containers of arrays which hold
beam properties for each ray. Geometrical beam properties such as x, y, z
(ray origins) and a, b, c (directional cosines) as well as polarization
characteristics depend on the above coordinate systems. Therefore, beams are
usually represented by two different objects: one in the global and one in a
local system.
Units¶
For the internal calculations, lengths are assumed to be in mm, although for reflection geometries and simple Bragg cases (thick crystals) this convention is not used. Angles are unitless (radians). Energy is in eV.
For plotting, the user may select units and conversion factors. The latter are usually automatically deduced from the units.
Beam categories¶
xrt discriminates rays by several categories:
good
: reflected within the working optical surface;out
: reflected outside of the working optical surface, i.e. outside of a metal stripe on a mirror;over
: propagated over the surface without intersection;dead
: arrived below the optical surface and thus absorbed by the OE.
This distinction simplifies the adjustment of entrance and exit slits. The
user supplies physical and optical limits, where the latter is used to
define the out
category (for rays between physical and optical limits).
An alarm is triggered if the fraction of dead rays exceeds a specified level.
Scripting in python¶
The user of raycing
must do the following:
 Instantiate class
BeamLine
and fill it with sources, optical elements, screens etc.  Create a modulelevel function that returns a dictionary of beams – the
instances of
Beam
. Assign this function to the module variable xrt.backends.raycing.run.run_process. The beams should be obtained by the methods shine() of a source, expose() of a screen, reflect() or multiple_reflect() of an optical element, propagate() of an aperture.  Use the keys in this dictionary for creating the plots (instances of
XYCPlot
). Note that at the time of instantiation the plots are just empty placeholders for the future 2D and 1D histograms.  Run
run_ray_tracing()
function for the created plots.
Additionally, the user may define a generator that will run a loop of ray
tracing for changing geometry (mimics a real scan) or for different material
properties etc. The generator should modify the beamline elements and output
file names of the plots before yield. After the yield the plots are ready
and the generator may use their fields, e.g. intensity or dE or dy or
others to prepare a scan plot. Typically, this sequence is within a loop; after
the loop the user may prepare the final scan plot using matplotlib
functionality. The generator is given to run_ray_tracing()
as a parameter.
See the supplied examples.

class
xrt.backends.raycing.
BeamLine
(azimuth=0.0, height=0.0, alignE='auto')¶ Container class for beamline components. It also defines the beam line direction and height.

__init__
(azimuth=0.0, height=0.0, alignE='auto')¶  azimuth: float
 Is counted in cw direction from the global Y axis. At azimuth = 0 the local Y coincides with the global Y.
 height: float
 Beamline height in the global system.
 alignE: float or ‘auto’
 Energy for automatic alignment in [eV]. If ‘auto’, alignment energy is defined as the middle of the Source energy range. Plays a role if the pitch or bragg parameters of the energy dispersive optical elements were set to ‘auto’.


xrt.backends.raycing.run.
run_process
(beamLine)¶
Sources¶
Geometric sources¶
Module sources
defines the container class
Beam
that holds beam properties as numpy arrays and defines the beam
sources. The sources are geometric sources in continuous and mesh forms and
synchrotron sources.
The intensity and polarization of each ray are determined via coherency matrix. This way is appropriate for incoherent addition of rays, in which the components of the coherency matrix of different rays are simply added.

class
xrt.backends.raycing.sources.
Beam
(nrays=100000, copyFrom=None, forceState=False, withNumberOfReflections=False, withAmplitudes=False, xyzOnly=False, bl=None)¶ Container for the beam arrays. x, y, z give the starting points. a, b, c give normalized vectors of ray directions (the source must take care about the normalization). E is energy. Jss, Jpp and Jsp are the components of the coherency matrix. The latter one is complex. Es and Ep are s and p field amplitudes (not always used). path is the total path length from the source to the last impact point. theta is the incidence angle. order is the order of grating diffraction. If multiple reflections are considered: nRefl is the number of reflections, elevationD is the maximum elevation distance between the rays and the surface as the ray travels from one impact point to the next one, elevationX, elevationY, elevationZ are the coordinates of the highest elevation points. If an OE uses a parametric representation, s, phi, r arrays store the impact points in the parametric coordinates.

class
xrt.backends.raycing.sources.
GeometricSource
¶ Implements a geometric source  a source with the ray origin, divergence and energy sampled with the given distribution laws.

__init__
(bl=None, name='', center=(0, 0, 0), nrays=100000, distx='normal', dx=0.32, disty=None, dy=0, distz='normal', dz=0.018, distxprime='normal', dxprime=0.001, distzprime='normal', dzprime=0.0001, distE='lines', energies=(9000.0, ), energyWeights=None, polarization='horizontal', filamentBeam=False, uniformRayDensity=False, pitch=0, yaw=0)¶ bl: instance of
BeamLine
name: str
 center: tuple of 3 floats
 3D point in global system
nrays: int
 distx, disty, distz, distxprime, distzprime:
 ‘normal’, ‘flat’, ‘annulus’ or None.
If is None, the corresponding arrays remain with the values got at
the instantiation of
Beam
. ‘annulus’ sets a uniform distribution for (x and z) or for (xprime and zprime) pairs. You can assign ‘annulus’ to only one member in the pair.  dx, dy, dz, dxprime, dzprime: float
 for normal distribution is sigma or (sigma, cut_limit), for flat is full width or tuple (min, max), for annulus is tuple (rMin, rMax), otherwise is ignored.
distE: ‘normal’, ‘flat’, ‘lines’, None
 energies: all in eV. (centerE, sigmaE) for distE = ‘normal’,
 (minE, maxE) for distE = ‘flat’, a sequence of E values for distE = ‘lines’
 energyWeights: 1D arraylike
 Can be used together with distE = ‘lines’ to specify the weight of each line. Must be of the shape of energies.
 polarization:
 ‘h[orizontal]’, ‘v[ertical]’, ‘+45’, ‘45’, ‘r[ight]’, ‘l[eft]’, None, custom. In the latter case the polarization is given by a tuple of 4 components of the coherency matrix: (Jss, Jpp, Re(Jsp), Im(Jsp)).
 filamentBeam: bool
 If True the source generates coherent monochromatic wavefronts. Required for the wave propagation calculations.
 uniformRayDensity: bool
 If True, the radiation is sampled uniformly but with varying amplitudes, otherwise with the density proportional to intensity and with constant amplitudes. Required as True for wave propagation calculations. False is usual for raytracing. This parameter only affects normal distributions, as for flat and annulus distributions the density is already uniform. If you set it True, the size parameter (dx or dz) must be given as (sigma, cut_limit).
 pitch, yaw: float
 rotation angles around x and z axis. Useful for canted sources.


class
xrt.backends.raycing.sources.
MeshSource
¶ Implements a point source representing a rectangular angular mesh of rays. Primarily, it is meant for internal usage for matching the maximum divergence to the optical sizes of optical elements.

__init__
(bl=None, name='', center=(0, 0, 0), minxprime=0.0001, maxxprime=0.0001, minzprime=0.0001, maxzprime=0.0001, nx=11, nz=11, distE='lines', energies=(9000.0, ), energyWeights=None, polarization='horizontal', withCentralRay=True, autoAppendToBL=False)¶ bl: instance of
BeamLine
name: str
 center: tuple of 3 floats
 3D point in global system
 minxprime, maxxprime, minzprime, maxzprime: float
 limits for the ungular distributions
 nx, nz: int
 numbers of points in x and z dircetions
distE: ‘normal’, ‘flat’, ‘lines’, None
 energies, all in eV: (centerE, sigmaE) for distE = ‘normal’,
 (minE, maxE) for distE = ‘flat’, a sequence of E values for distE = ‘lines’
 energyWeights: 1D arraylike
 Can be used together with distE = ‘lines’ to specify the weight of each line. Must be of the shape of energies.
 polarization: ‘h[orizontal]’, ‘v[ertical]’, ‘+45’, ‘45’, ‘r[ight]’,
 ‘l[eft]’, None, custom. In the latter case the polarization is given by a tuple of 4 components of the coherency matrix: (Jss, Jpp, Re(Jsp), Im(Jsp)).
 withCentralRay: bool
 if True, the 1st ray in the beam is along the nominal beamline direction
 autoAppendToBL: bool
 if True, the source is added to the list of beamline sources.
Otherwise the user must manually start it with
shine()
.


class
xrt.backends.raycing.sources.
CollimatedMeshSource
(bl=None, name='', center=(0, 0, 0), dx=1.0, dz=1.0, nx=11, nz=11, distE='lines', energies=(9000.0, ), energyWeights=None, polarization='horizontal', withCentralRay=True, autoAppendToBL=False)¶ Implements a source representing a mesh of collimated rays. Is similar to
MeshSource
.

class
xrt.backends.raycing.sources.
GaussianBeam
(bl=None, name='', center=(0, 0, 0), w0=0.1, distE='lines', energies=(9000.0, ), energyWeights=None, polarization='horizontal', pitch=0, yaw=0)¶ Implements a Gaussian beam https://en.wikipedia.org/wiki/Gaussian_beam. It must be used for an already available set of 3D points which are obtained by
prepare_wave()
of a slit, oe or screen. See a usage example in\tests\raycing\laguerre_hermite_gaussian_beam.py
.
__init__
(bl=None, name='', center=(0, 0, 0), w0=0.1, distE='lines', energies=(9000.0, ), energyWeights=None, polarization='horizontal', pitch=0, yaw=0)¶ bl: instance of
BeamLine
name: str
 center: tuple of 3 floats
 3D point in global system
 w0: float or 2sequence
 Gaussian beam waist size. If a 2sequence, the sizes refer to the horizontal and the vertical axes.
distE: ‘normal’, ‘flat’, ‘lines’, None
 energies: all in eV. (centerE, sigmaE) for distE = ‘normal’,
 (minE, maxE) for distE = ‘flat’, a sequence of E values for distE = ‘lines’
 energyWeights: 1D arraylike
 Can be used together with distE = ‘lines’ to specify the weight of each line. Must be of the shape of energies.
 polarization:
 ‘h[orizontal]’, ‘v[ertical]’, ‘+45’, ‘45’, ‘r[ight]’, ‘l[eft]’, None, custom. In the latter case the polarization is given by a tuple of 4 components of the coherency matrix: (Jss, Jpp, Re(Jsp), Im(Jsp)).
 pitch, yaw: float
 rotation angles around x and z axis. Useful for canted sources.


class
xrt.backends.raycing.sources.
LaguerreGaussianBeam
(*args, **kwargs)¶ Implements LaguerreGaussian beam https://en.wikipedia.org/wiki/Gaussian_beam. It must be used for an already available set of 3D points which are obtained by
prepare_wave()
of a slit, oe or screen. See a usage example in\tests\raycing\laguerre_hermite_gaussian_beam.py
.
__init__
(*args, **kwargs)¶  vortex: None or tuple(l, p)
 specifies LaguerreGaussian beam with l the azimuthal index and p >= 0 the radial index.


class
xrt.backends.raycing.sources.
HermiteGaussianBeam
(*args, **kwargs)¶ Implements HermiteGaussian beam https://en.wikipedia.org/wiki/Gaussian_beam. It must be used for an already available set of 3D points which are obtained by
prepare_wave()
of a slit, oe or screen. See a usage example in\tests\raycing\laguerre_hermite_gaussian_beam.py
.
__init__
(*args, **kwargs)¶  TEM: None or tuple(m, n)
 specifies HermiteGaussian beam of order (m, n) referring to the x and y directions.


xrt.backends.raycing.sources.
make_energy
(distE, energies, nrays, filamentBeam=False, energyWeights=None)¶ Creates energy distributions with the distribution law given by distE. energies either determine the limits or is a sequence of discrete energies. If distE is ‘lines’, energyWeights can define the relative weights of the lines.

xrt.backends.raycing.sources.
make_polarization
(polarization, bo, nrays=100000)¶ Initializes the coherency matrix. The following polarizations are supported:
horizontal (polarization is a string started with ‘h’):
\[\begin{split}J = \left( \begin{array}{ccc}1 & 0 \\ 0 & 0\end{array} \right)\end{split}\]vertical (polarization is a string started with ‘v’):
\[\begin{split}J = \left( \begin{array}{ccc}0 & 0 \\ 0 & 1\end{array} \right)\end{split}\]at +45º (polarization = ‘+45’):
\[\begin{split}J = \frac{1}{2} \left( \begin{array}{ccc}1 & 1 \\ 1 & 1\end{array} \right)\end{split}\]at 45º (polarization = ‘45’):
\[\begin{split}J = \frac{1}{2} \left( \begin{array}{ccc}1 & 1 \\ 1 & 1\end{array} \right)\end{split}\]right (polarization is a string started with ‘r’):
\[\begin{split}J = \frac{1}{2} \left( \begin{array}{ccc}1 & i \\ i & 1\end{array} \right)\end{split}\]left (polarization is a string started with ‘l’):
\[\begin{split}J = \frac{1}{2} \left( \begin{array}{ccc}1 & i \\ i & 1\end{array} \right)\end{split}\]
unpolarized (polarization is None):
\[\begin{split}J = \frac{1}{2} \left( \begin{array}{ccc}1 & 0 \\ 0 & 1\end{array} \right)\end{split}\]userdefined (polarization is 4sequence with):
\[\begin{split}J = \left( \begin{array}{ccc} {\rm polarization[0]} & {\rm polarization[2]} + i * {\rm polarization[3]} \\ {\rm polarization[2]}  i * {\rm polarization[3]} & {\rm polarization[1]}\end{array} \right)\end{split}\]

xrt.backends.raycing.sources.
rotate_coherency_matrix
(beam, indarr, roll)¶ Rotates the coherency matrix \(J\):
\[\begin{split}J = \left( \begin{array}{ccc} J_{ss} & J_{sp} \\ J^*_{sp} & J_{pp}\end{array} \right)\end{split}\]by angle \(\phi\) around the beam direction as \(J' = R_{\phi} J R^{1}_{\phi}\) with the rotation matrix \(R_{\phi}\) defined as:
\[\begin{split}R_{\phi} = \left( \begin{array}{ccc} \cos{\phi} & \sin{\phi} \\ \sin{\phi} & \cos{\phi}\end{array} \right)\end{split}\]

xrt.backends.raycing.sources.
shrink_source
(beamLine, beams, minxprime, maxxprime, minzprime, maxzprime, nx, nz)¶ Utility function that does ray tracing with a mesh source and shrinks its divergence until the footprint beams match to the optical surface. Parameters:
beamline: instance of
BeamLine
beams: tuple of str
Dictionary keys in the result ofrun_process()
corresponding to the wanted footprints.minxprime, maxxprime, minzprime, maxzprime: float
Determines the size of the mesh source. This size can only be shrunk, not expanded. Therefore, you should provide it sufficiently big for your needs. Typically, min values are negative and max values are positive.nx, nz: int
Sizes of the 2D mesh grid in x and z direction.Returns an instance of
MeshSource
which can be used then for getting the divergence values.
Synchrotron sources¶
Note
In this section we consider zaxis to be directed along the beamline in order to be compatible with the cited works. Elsewhere in xrt zaxis looks upwards.
The synchrotron sources have two implementations: based on own fully pythonic or OpenCL aided calculations and based on external codes [Urgent] and [WS]. The latter codes have some drawbacks, as demonstrated in the section Comparison of synchrotron source codes, but nonetheless can be used for comparison purposes. If you are going to use them, the codes are freely available as parts of [XOP] package.
[Urgent]  (1, 2, 3, 4) R. P. Walker, B. Diviacco, URGENT, A computer program for calculating undulator radiation spectral, angular, polarization and power density properties, STM9112B, July 1991, Presented at 4th Int. Conf. on Synchrotron Radiation Instrumentation, Chester, England, Jul 1519, 1991 
[WS]  R. J. Dejus, Computer simulations of the wiggler spectrum, Nucl. Instrum. Methods Phys. Res., Sect. A, 347 (1994) 5660 
[XOP]  (1, 2, 3, 4, 5) M. Sanchez del Rio, R. J. Dejus, XOP: A Multiplatform Graphical User Interface for Synchrotron Radiation Spectral and Optics Calculations, SPIE Proc. 3152 (1997) 148; web page: http://www.esrf.eu/Instrumentation/software/dataanalysis/xop2.3. 
The internal synchrotron sources are based on the following works: [Kim] and [Walker]. We use the general formulation for the flux distribution in 3dimensional phase space (solid angle and energy) [Kim]:
\[\begin{split}\mathcal{F}(\theta,\psi,E) &= \alpha\frac{\Delta \omega}{\omega} \frac{I_e}{e^{}}(A_{\sigma}^2 + A_{\pi}^2)\\\end{split}\]
For the bending magnets the amplitudes can be calculated analytically using the modified Bessel functions \(K_v(y)\):
\[\begin{split}\begin{bmatrix} A_{\sigma}\\ A_{\pi} \end{bmatrix} &= \frac{\sqrt{3}}{2\pi}\gamma\frac{\omega}{\omega_c} (1+\gamma^2\psi^2)\begin{bmatrix}i K_{2/3}(\eta)\\ \frac{\gamma\psi}{\sqrt{1+\gamma^2\psi^2}} K_{1/3}(\eta)\end{bmatrix}\end{split}\]
 where
 \[\begin{split}\gamma &= \frac{E_e}{m_{e}c^2} = 1957E_e[{\rm GeV}]\\ \eta &= \frac{1}{2}\frac{\omega}{\omega_c}(1+\gamma^2\psi^2)^{3/2}\\ \omega_c &= \frac{3\gamma^{3}c}{2\rho}\\ \rho &= \frac{m_{e}c\gamma}{e^{}B}\end{split}\]
with \(I_e\)  the current in the synchrotron ring, \(B\)  magnetic field strength, \(e^{}, m_e\) the electron charge and mass, \(c\)  the speed of light.
Wiggler radiation relies on the same equation considering each pole as a bending magnet with variable magnetic field/curvature radius: \(\rho(\theta) = \sin(\arccos(\theta\gamma/K))\), where \(K\) is deflection parameter. Total flux is multiplied then by \(2N\), where \(N\) is the number of wiggler periods.
For the undulator sources the amplitude integrals must be calculated numerically, starting from the magnetic field.
\[\begin{split}\begin{bmatrix} A_{\sigma}\\ A_{\pi} \end{bmatrix} &= \frac{1}{2\pi}\int\limits_{\infty}^{+\infty}dt' \begin{bmatrix}\frac{[\textbf{n}\times[(\textbf{n} \boldsymbol{\beta})\times\dot{\boldsymbol{\beta}}]]} {(1  \textbf{n}\cdot\boldsymbol{\beta})^2}\end{bmatrix}_{x, y} e^{i\omega (t' + R(t')/c)}\end{split}\]\[\begin{split}B_x &= B_{x0}\sin(2\pi z /\lambda_u + \phi),\\ B_y &= B_{y0}\sin(2\pi z /\lambda_u)\end{split}\]
the corresponding velosity components are
\[\begin{split}\beta_x &= \frac{K_y}{\gamma}\cos(\omega_u t),\\ \beta_y &= \frac{K_x}{\gamma}\cos(\omega_u t + \phi)\\ \beta_z &= \sqrt{1\frac{1}{\gamma^2}\beta_{x}^{2}\beta_{y}^{2}},\end{split}\]
where \(\omega_u = 2\pi c /\lambda_u\)  undulator frequency, \(\phi\)  phase shift between the magnetic field components. In this simple case one can consider only one period in the integral phase term replacing the exponential series by a factor \(\frac{\sin(N\pi\omega/\omega_1)}{\sin(\pi\omega/\omega_1)}\), where \(\omega_1 = \frac{2\gamma^2}{1+K_x^2/2+K_y^2/2+\gamma^2(\theta^2+\psi^2)} \omega_u\).
In the case of tapered undulator, the vertical magnetic field is multiplied by an additional factor \(1  \alpha z\), that in turn results in modification of horizontal velocity and coordinate.
In the farfield approximation we consider the undulator as a point source and replace the distance \(R\) by a projection \(\mathbf{n}\cdot\mathbf{r}\), where \(\mathbf{n} = [\theta,\psi,1\frac{1}{2}(\theta^2+\psi^2)]\)  direction to the observation point, \(\mathbf{r} = [\frac{K_y}{\gamma}\sin(\omega_u t'), \frac{K_x}{\gamma}\sin(\omega_u t' + \phi)\), \(\beta_m\omega t'\frac{1}{8\gamma^2} (K_y^2\sin(2\omega_u t')+K_x^2\sin(2\omega_u t'+2\phi))]\)  electron trajectory, \(\beta_m = 1\frac{1}{2\gamma^2}\frac{K_x^2}{4\gamma^2} \frac{K_y^2}{4\gamma^2}\).
Configurations with nonequivalent undulator periods i.e. tapered undulator require integration over full undulator length, similar approach is used for the near field calculations, where the undulator extension is taken into account and the phase component in the integral is taken in its initial form \(i\omega (t' + R(t')/c)\).
For the custom field configuratiuons, where the magnetic field components are tabulated as functions of longitudinal coordinate \(\textbf{B}=[B_{x}(z), B_{y}(z), B_{z}(z))]\), preliminary numerical calculation of the velocity and coordinate is nesessary. For that we solve the system of differential equations in the trajectory coordinate \(s\):
\[\begin{split}\frac{d^2}{ds^2} \begin{bmatrix}x\\ y\\ z\end{bmatrix} &= \frac{e^{}}{\gamma m_{e} c} \begin{bmatrix}\beta_{y}B_{z}B_{y}\\ \beta_{x}B_{z}+B_{x}\\ \beta_{y}B_{x}+\beta_{x}B_{y}\end{bmatrix}\end{split}\]
using the classical RungeKutta fourthorder method.
For the Undulator and custom field models we directly calculate the integral
using the ClenshawCurtis quadrature, it proves
to converge as quickly as the previously used GaussLegendre method, but the
nodes and weights are calculated significantly faster. The size
of the integration grid is evaluated at the points of slowest convergence
(highest energy, maximum angular deviation i.e. a corner of the plot) before
the start of intensity map calculation and then applied to all points.
This approach creates certain computational overhead for the onaxis/low energy
parts of the distribution but enables efficient parallelization and gives
significant overall gain in performance. Initial evaluation typically takes
just a few seconds, but might get much longer for custom magnetic fields and
near edge calculations. If such heavy task is repeated many times for the given
angular and energy limits it might make sense to note the evaluated size of
the grid during the first run or call the test_convergence()
method once,
and then use the fixed grid by defining the gNodes at the init.
Note also that the grid size will be automatically reevaluated if any of the
limits/electron energy/undulator deflection parameter or period length are
redefined dynamically in the script.
Typical convergence threshold is defined by machine precision multiplied by the
size of the integration grid. Default integration parameters proved to work
very well in most cases, but may fail if the angular and/or energy limits are
unreasonably wide. If in doubt, check the convergence
with test_convergence()
. See also
a convergence study that justifies our automatic grid
evaluation.
For the purpose of ray tracing (and this is not necessary for wave propagation) the undulator source size is calculated following [TanakaKitamura]. Their formulation includes dependence on electron beam energy spread. The effective linear and angular source sizes are given by
\[\begin{split}\Sigma &= \left(\sigma_e^2 + \sigma_r^2\right)^{1/2} =\left(\varepsilon\beta + \frac{\lambda L}{2\pi^2} [Q(\sigma_\epsilon/4)]^{4/3}\right)^{1/2}\\ \Sigma' &= \left({\sigma'_e}^2 + {\sigma'_r}^2\right)^{1/2} =\left(\varepsilon/\beta + \frac{\lambda}{2L} [Q(\sigma_\epsilon)]^2\right)^{1/2},\end{split}\]
where \(\varepsilon\) and \(\beta\) are the electron beam emittance and betatron function, the scaling function \(Q\) is defined as
\[Q(x) = \left(\frac{2x^2}{\exp(2x^2)+(2\pi)^{1/2}x\ {\rm erf} (2^{1/2}x)1}\right)^{1/2}\]
(notice \(Q(0)=1\)) and \(\sigma_\epsilon\) is the normalized energy spread
\[\sigma_\epsilon = 2\pi nN\sigma_E\]
i.e. the energy spread \(\sigma_E\) divided by the undulator bandwidth \(1/nN\) of the nth harmonic, with an extra factor \(2\pi\). See an application example here.
Note
If you want to compare the source size with that by [SPECTRA], note that their radiation source size \(\sigma_r\) is by a factor of 2 smaller in order to be compatible with the traditional formula by [Kim]. In this aspect SPECTRA contradicts to their own paper [TanakaKitamura], see the paragraph after Eq(23).
[Kim]  (1, 2, 3) K.J. Kim, Characteristics of Synchrotron Radiation, AIP Conference Proceedings, 184 (AIP, 1989) 565. 
[Walker]  R. Walker, Insertion devices: undulators and wigglers, CAS  CERN Accelerator School: Synchrotron Radiation and Free Electron Lasers, Grenoble, France, 2227 Apr 1996: proceedings (CERN. Geneva, 1998) 129190. 
[TanakaKitamura]  (1, 2, 3) T. Tanaka and H. Kitamura, Universal function for the brilliance of undulator radiation considering the energy spread effect, J. Synchrotron Rad. 16 (2009) 380–6. 

class
xrt.backends.raycing.sources.
UndulatorUrgent
¶ Undulator source that uses the external code Urgent. It has some drawbacks, as demonstrated in the section Comparison of synchrotron source codes, but nonetheless can be used for comparison purposes. If you are going to use it, the code is freely available as part of XOP package.

__init__
(bl=None, name='UrgentU', center=(0, 0, 0), nrays=100000, period=32.0, K=2.668, Kx=0.0, Ky=0.0, n=12, eE=6.0, eI=0.1, eSigmaX=134.2, eSigmaZ=6.325, eEpsilonX=1.0, eEpsilonZ=0.01, uniformRayDensity=False, eMin=1500, eMax=81500, eN=1000, eMinRays=None, eMaxRays=None, xPrimeMax=0.25, zPrimeMax=0.1, nx=25, nz=25, path=None, mode=4, icalc=1, useZip=True, order=3, processes='auto')¶ The 1st instantiation of this class runs the Urgent code and saves its output into a “.pickle” file. The temporary directory “tmp_urgent” can then be deleted. If any of the Urgent parameters has changed since the previous run, the Urgent code is forced to redo the calculations.
 bl: instance of
BeamLine
 Container for beamline elements. Sourcess are added to its sources list.
 name: str
 Userspecified name, can be used for diagnostics output.
 center: tuple of 3 floats
 3D point in global system
 nrays: int
 the number of rays sampled in one iteration
 period: float
 Magnet period (mm).
 K or Ky: float
 Magnet deflection parameter (Ky) in the vertical field.
 Kx: float
 Magnet deflection parameter in the horizontal field.
 n: int
 Number of magnet periods.
 eE: float
 Electron beam energy (GeV).
 eI: float
 Electron beam current (A).
 eSigmaX, eSigmaZ: float
 rms horizontal and vertical electron beam sizes (µm).
 eEpsilonX, eEpsilonZ: float
 Horizontal and vertical electron beam emittance (nm rad).
 eMin, eMax: float
 Minimum and maximum photon energy (eV) used by Urgent.
 eN: int
 Number of photon energy intervals used by Urgent.
 eMinRays, eMaxRays: float
 The range of energies for rays. If None, are set equal to eMin and eMax. These two parameters are useful for playing with the energy axis without having to force Urgent to redo the calculations each time.
 xPrimeMax, zPrimeMax: float
 Half of horizontal and vertical acceptance (mrad).
 nx, nz: int
 Number of intervals in the horizontal and vertical directions from zero to maximum.
 path: str
 Full path to Urgent executable. If None, it is set automatically from the module variable xopBinDir.
 mode: 1, 2 or 4
 the MODE parameter of Urgent. If =1, UndulatorUrgent scans energy and reads the xz distribution from Urgent. If =2 or 4, UndulatorUrgent scans x and z and reads energy spectrum (angular density for 2 or flux through a window for 4) from Urgent. The meshes for x, z, and E are restricted in Urgent: nx,nz<50 and nE<5000. You may overcome these restrictions if you scan the corresponding quantities outside of Urgent, i.e. inside of this class UndulatorUrgent. mode = 4 is by far most preferable.
 icalc: int
 The ICALC parameter of Urgent.
 useZip: bool
 Use gzip module to compress the output files of Urgent. If True, the temporary storage takes much less space but a slightly bit more time.
 order: 1 or 3
 the order of the spline interpolation. 3 is recommended.
 processes: int or any other type as ‘auto’
 the number of worker processes to use. If the type is not int then the number returned by multiprocessing.cpu_count()/2 is used.
 bl: instance of


class
xrt.backends.raycing.sources.
WigglerWS
¶ Wiggler source that uses the external code ws. It has some drawbacks, as demonstrated in the section Comparison of synchrotron source codes, but nonetheless can be used for comparison purposes. If you are going to use it, the code is freely available as part of XOP package.

__init__
(*args, **kwargs)¶ Uses WS code. All the parameters are the same as in UndulatorUrgent.


class
xrt.backends.raycing.sources.
BendingMagnetWS
¶ Bending magnet source that uses the external code ws. It has some drawbacks, as demonstrated in the section Comparison of synchrotron source codes, but nonetheless can be used for comparison purposes. If you are going to use it, the code is freely available as parts of XOP package.

__init__
(*args, **kwargs)¶ Uses WS code.
 B0: float
 Field in Tesla.
 K, n, period and nx:
 Are set internally.
The other parameters are the same as in UndulatorUrgent.


class
xrt.backends.raycing.sources.
SourceBase
¶ Base class for the Synchrotron Sources. Not to be called explicitly.

__init__
(bl=None, name='GenericSource', center=(0, 0, 0), nrays=100000, eE=6.0, eI=0.1, eEspread=0.0, eSigmaX=None, eSigmaZ=None, eEpsilonX=1.0, eEpsilonZ=0.01, betaX=9.0, betaZ=2.0, eMin=5000.0, eMax=15000.0, distE='eV', xPrimeMax=0.5, zPrimeMax=0.5, R0=None, uniformRayDensity=False, filamentBeam=False, pitch=0, yaw=0, eN=51, nx=25, nz=25)¶  bl: instance of
BeamLine
 Container for beamline elements. Sourcess are added to its sources list.
 name: str
 Userspecified name, can be used for diagnostics output.
 center: tuple of 3 floats
 3D point in global system.
 nrays: int
 The number of rays sampled in one iteration.
 eE: float
 Electron beam energy (GeV).
 eI: float
 Electron beam current (A).
 eEspread: float
 Energy spread relative to the beam energy, rms.
 eSigmaX, eSigmaZ: float
 rms horizontal and vertical electron beam sizes (µm). Alternatively, betatron functions can be specified instead of the electron beam sizes.
 eEpsilonX, eEpsilonZ: float
 Horizontal and vertical electron beam emittance (nm rad).
 betaX, betaZ:
 Betatron function (m). Alternatively, beam size can be specified.
 R0: float
 Distance centertoscreen for the near field calculations (mm). If None, the far field approximation (i.e. “usual” calculations) is used.
 eMin, eMax: float
 Minimum and maximum photon energy (eV). Used as band width for flux calculation.
 distE: ‘eV’ or ‘BW’
 The resulted flux density is per 1 eV or 0.1% bandwidth. For ray tracing ‘eV’ is used.
 xPrimeMax, zPrimeMax:
Horizontal and vertical acceptance (mrad).
Note
The Monte Carlo sampling of the rays having their density proportional to the beam intensity can be extremely inefficient for sharply peaked distributions, like the undulator angular density distribution. It is therefore very important to restrict the sampled angular acceptance down to very small angles. Use this source only with reasonably small xPrimeMax and zPrimeMax!
 uniformRayDensity: bool
 If True, the radiation is sampled uniformly but with varying amplitudes, otherwise with the density proportional to intensity and with constant amplitudes. Required as True for wave propagation calculations. False is usual for raytracing.
 filamentBeam: bool
 If True the source generates coherent monochromatic wavefronts. Required as True for the wave propagation calculations in partially coherent regime.
 pitch, yaw: float
 rotation angles around x and z axis. Useful for canted sources.
 bl: instance of

intensities_on_mesh
(energy='auto', theta='auto', psi='auto', harmonic=None, eSpreadSigmas=3.5, eSpreadNSamples=36)¶ Returns the Stokes parameters in the shape (energy, theta, psi, [harmonic]), with theta being the horizontal mesh angles and psi the vertical mesh angles. Each one of the input parameters is a 1D array of an individually selectable length.
Note
We do not provide any internal mesh optimization, as mesh functions are not our core objectives. In particular, the angular meshes must be wider than the electron beam divergences in order to convolve the field distribution with the electron distribution. A warning will be printed (new in version 1.3.4) if the requested meshes are too narrow.

multi_electron_stack
(energy='auto', theta='auto', psi='auto', harmonic=None, withElectronDivergence=True)¶ Returns Es and Ep in the shape (energy, theta, psi, [harmonic]). Along the 0th axis (energy) are stored “macroelectrons” that emit at the photon energy given by energy (constant or variable) onto the angular mesh given by theta and psi. The transverse field from each macroelectron gets individual random angular offsets dtheta and dpsi within the emittance distribution if withElectronDivergence is True and an individual random shift to gamma within the energy spread. The parameter self.filamentBeam is irrelevant for this method.

real_photon_source_sizes
(energy='auto', theta='auto', psi='auto', method='rms')¶ Returns energy dependent arrays: flux, (dx’)², (dz’)², dx², dz². Depending on distE being ‘eV’ or ‘BW’, the flux is either in ph/s or in ph/s/0.1%BW, being integrated over the specified theta and psi ranges. The squared angular and linear photon source sizes are variances, i.e. squared sigmas. The latter two (linear sizes) are in mm**2.


class
xrt.backends.raycing.sources.
IntegratedSource
¶ Base class for the Sources with numerically integrated amplitudes:
SourceFromField
andUndulator
. Not to be called explicitly.
__init__
(*args, **kwargs)¶  gp: float
 Defines the relative precision of the integration (last significant digit). Undulator model converges down to 1e6 and below. Custom field calculation may require setting the precision of 1e3.
 gNodes: int
 Number of integration nodes in each of the integration intervals. If not provided at init, will be defined automatically.
 targetOpenCL: None, str, 2tuple or tuple of 2tuples
assigns the device(s) for OpenCL accelerated calculations. None, if pyopencl is not wanted. Ignored if pyopencl is not installed. Accepts the following values:
 a tuple (iPlatform, iDevice) of indices in the
lists
cl.get_platforms()
andplatform.get_devices()
, see the section Calculations on GPU.  a tuple of tuples ((iP1, iD1), …, (iPn, iDn)) to assign specific devices from one or multiple platforms.
 int iPlatform  assigns all devices found at the given platform.
 ‘GPU’  lets the program scan the system and select all found GPUs.
 ‘CPU’  similar to ‘GPU’. If one CPU exists in multiple platforms the program tries to select the vendorspecific driver.
 ‘other’  similar to ‘GPU’, used for Intel PHI and other OpenCL capable accelerator boards.
 ‘all’  lets the program scan the system and assign all found devices. Not recommended, since the performance will be limited by the slowest device.
 ‘auto’  lets the program scan the system and make an assignment according to the priority list: ‘GPU’, ‘other’, ‘CPU’ or None if no devices were found. Used by default.
 ‘SERVER_ADRESS:PORT’  calculations will be run on remote
server. See
tests/raycing/RemoteOpenCLCalculation
.
 a tuple (iPlatform, iDevice) of indices in the
lists
Warning
A good graphics or dedicated accelerator card is highly recommended! Special cases as wigglers by the undulator code, near field, wide angles and tapering are hardly doable on CPU.
Note
Consider the warnings and tips on using xrt with GPUs.
 precisionOpenCL: ‘float32’ or ‘float64’, only for GPU.
 Single precision (float32) should be enough in most cases. The calculations with doube precision are much slower. Double precision may be unavailable on your system. Tapering and Near Field calculations require double precision.

shine
(toGlobal=True, withAmplitudes=True, fixedEnergy=False, wave=None, accuBeam=None)¶ Returns the source beam. If toGlobal is True, the output is in the global system. If withAmplitudes is True, the resulted beam contains arrays Es and Ep with the s and p components of the electric field.
fixedEnergy is either None or a value in eV. If fixedEnergy is specified, the energy band is not 0.1%BW relative to fixedEnergy, as probably expected but is given by (eMax  eMin) of the constructor.
wave and accuBeam are used in wave diffraction. wave is a Beam object and determines the positions of the wave samples. It must be obtained by a previous
prepare_wave
run. accuBeam is only needed with several repeats of diffraction integrals when the parameters of the filament beam must be preserved for all the repeats.

test_convergence
(nMax=500000, withPlots=True, overStep=100)¶ This function evaluates the length of the integration grid required for convergence.
 nMax: int
 Maximum number of nodes.
 withPlots: bool
 Enables visualization.
 overStep: int
 Defines the number of extra points to calculate when the convergence is found. If None, calculation will proceed till nMax.


class
xrt.backends.raycing.sources.
Undulator
¶ Undulator source. The computation is volumnous and thus a decent GPU is highly recommended.

__init__
(*args, **kwargs)¶  period, n:
 Magnetic period (mm) length and number of periods.
 K, Kx, Ky: float
 Deflection parameter for the vertical field or for an elliptical undulator.
 B0x, B0y: float
 Maximum magnetic field. If both K and B provided at the init, K value will be used.
 phaseDeg: float
 Phase difference between horizontal and vertical magnetic arrays. Used in the elliptical case where it should be equal to 90 or 90.
 taper: tuple(dgap(mm), gap(mm))
 Linear variation in undulator gap. None if tapering is not used. Pyopencl is recommended for tapering.
 targetE: a tuple (Energy, harmonic{, isElliptical})
 Can be given for automatic calculation of the deflection parameter. If isElliptical is not given, it is assumed as False (as planar).
 xPrimeMaxAutoReduce, zPrimeMaxAutoReduce: bool
 Whether to reduce too large angular ranges down to the feasible values in order to improve efficiency. It is highly recommended to keep them True.

get_SIGMA
(E, onlyOddHarmonics=True, with0eSpread=False)¶ Calculates total linear source size, also including the effect of electron beam energy spread. Uses Tanaka and Kitamura, J. Synchrotron Rad. 16 (2009) 380–6.
E can be a value or an array. Returns a 2tuple with x and y sizes.

get_SIGMAP
(E, onlyOddHarmonics=True, with0eSpread=False)¶ Calculates total angular source size, also including the effect of electron beam energy spread. Uses Tanaka and Kitamura, J. Synchrotron Rad. 16 (2009) 380–6.
E can be a value or an array. Returns a 2tuple with x and y sizes.

power_vs_K
(energy, theta, psi, harmonics, Ks)¶ Calculates power curve – total power in W for all harmomonics at given K values (Ks). The power is calculated through the aperture defined by theta and psi opening angles within the energy range.
The result of this numerical integration depends on the used angular and energy meshes; you should check convergence. Internally, electron beam energy spread is also sampled by adding another dimension to the intensity array and making it 5dimensional. You therefore may want to set energy spread to zero, it doesn’t affect the resulting power anyway.
Returns a 1D array corresponding to Ks.

tuning_curves
(energy, theta, psi, harmonics, Ks)¶ Calculates tuning curves – maximum flux of given harmomonics at given K values (Ks). The flux is calculated through the aperture defined by theta and psi opening angles (1D arrays).
Returns two 2D arrays: energy positions and flux values. The rows correspond to Ks, the colums correspond to harmomonics.


class
xrt.backends.raycing.sources.
SourceFromField
¶ Dedicated class for the sources based on a custom field table.

__init__
(*args, **kwargs)¶  customField: float or str or tuple(fileName, kwargs) or numpy array.
 If float, adds a constant longitudinal field.
If str or tuple, expects table of field samples given as an Excel
file or as text file. If given as a tuple or list, the 2nd member
is a key word dictionary for reading Excel by
pandas.read_excel()
or reading text file bynumpy.loadtxt()
, e.g.dict(skiprows=4)
for skipping the file header. The file must contain the columns with longitudinal coordinate in mm, {B_hor,} B_ver {, B_long}, all in T. The field can be provided as a numpy array with the same structure as the table from file.


class
xrt.backends.raycing.sources.
Wiggler
¶ Wiggler source. The computation is reasonably fast and thus a GPU is not required and is not implemented.

__init__
(*args, **kwargs)¶ Parameters are the same as in BendingMagnet except B0 and rho which are not required and additionally:
 K: float
 Deflection parameter
 period: float
 period length in mm.
 n: int
 Number of periods.


class
xrt.backends.raycing.sources.
BendingMagnet
¶ Bending magnet source. The computation is reasonably fast and thus a GPU is not required and is not implemented.

__init__
(*args, **kwargs)¶  B0: float
 Magnetic field (T). Alternatively, specify rho.
 rho: float
 Curvature radius (m). Alternatively, specify B0.

Comparison of synchrotron source codes¶
Using xrt synchrotron sources on a mesh¶
The main modus operandi of xrt synchrotron sources is to provide Monte Carlo rays or wave samples. For comparing our sources with other codes – all of them are fully deterministic, being defined on certain meshes – we also supply mesh methods such as intensities_on_mesh, power_vs_K and tuning_curves. Note that we do not provide any internal mesh optimization, as these mesh functions are not our core objectives. Instead, the user themself should care about the proper mesh limits and step sizes. In particular, the angular meshes must be wider than the electron beam divergences in order to convolve the field distribution with the electron distribution of nonzero emittance. The above mentioned mesh methods will print a warning (new in version 1.3.4) if the requested meshes are too narrow.
If you want to calculate flux through a narrow aperture, you first calculate intensities_on_mesh on wide enough angular meshes and then slice the intensity down to the needed aperture size. An example of such calculations is given in tests/raycing/test_undulator_on_mesh.py which produces the following plot (for a BESSY undulator, zero energy spread, as Urgent cannot take account of it):
Undulator spectrum across a harmonic¶
The above six classes result in ray distributions with the density of rays proportional to intensity. This requires an algorithm of Monte Carlo ray sampling that needs a 3D (two directions + energy) intensity distribution. The classes using the external meshbased codes interpolate the source intensity precalculated over the user specified mesh. The latter three classes (internal implementation of synchrotron sources) do not require interpolation, which eliminates two problems: artefacts of interpolation and the need for the mesh optimization. However, this requires the calculations of intensity for each ray.
For bending magnet and wiggler sources these calculations are not heavy and are actually faster than 3D interpolation. See the absence of interpolation artefacts in Synchrotron sources in the gallery.
For an undulator the calculations are much more demanding and for a wide angular acceptance the Monte Carlo ray sampling can be extremely inefficient. To improve the efficiency, a reasonably small acceptance should be considered.
There are several codes that can calculate undulator spectra: [Urgent], [US], [SPECTRA]. There is a common problem about them that the energy spectrum may get strong distortions if calculated with a sparse spatial and energy mesh. SPECTRA code seems to provide the best reference for undulator spectra, which was used to optimize the meshes of the other codes. The final comparison of the resulted undulator spectra around a single harmonic is shown below.
[US]  (1, 2) R. J. Dejus, US: Program to calculate undulator spectra. Part of XOP. 
[SPECTRA]  (1, 2, 3, 4, 5) T. Tanaka and H. Kitamura, SPECTRA  a synchrotron radiation calculation code, J. Synchrotron Radiation 8 (2001) 12218. 
Note
If you are going to use UndulatorUrgent, you should optimize the spatial and energy meshes! The resulted ray distribution is strongly dependent on them, especially on the energy mesh. Try different numbers of points and various energy ranges.
Undulator spectrum at very high energies¶
The codes [Urgent] and [SPECTRA] result in saturation at high energies (see the image below) thus leading to a divergent total power integral. The false radiation has a circular offaxis shape. To the contrary, xrt and [SRW] flux at high energies vanishes and follows the wiggler approximation. More discussion will follow in a future journal article about xrt.
Singleelectron and multielectron undulator radiation¶
Here we compare singleelectron and multielectron (i.e. with a finite electron beam size and energy spread) undulator radiation, as calculated by xrt and [SRW]. The calculations are done on a 3D mesh of energy (the long axis on the images below) and two transverse angles. Notice also the duration of execution given below each image. The 3D mesh was the following: theta = 321 point, 0.3 to +0.3 mrad, psi = 161 point, 0.15 to +0.15 mrad, energy: 301 point 1.5 to 4.5 keV.
[SRW]  (1, 2, 3) O. Chubar, P. Elleaume, Accurate And Efficient Computation Of Synchrotron Radiation In The Near Field Region, proc. of the EPAC98 Conference, 2226 June 1998, p.11771179. 
Undulator spectrum with tapering¶
The spectrum can be broadened by tapering the magnetic gap. The figure below shows a comparison of xrt with [SPECTRA], [YAUP] and experimental measurements [exp_taper]. The gap values and the taper were slightly varied in all three codes to reach the best match with the experimental curve. We had to do so because in the other codes taper is not clearly defined (where is the gap invariant – at the center or at one of the ends?) and also because the nominal experimental gap and taper are not fully trustful.
[YAUP]  B. I. Boyanov, G. Bunker, J. M. Lee, and T. I. Morrison, Numerical Modeling of Tapered Undulators, Nucl. Instr. Meth. A339 (1994) 596603. 
[exp_taper]  Measured on 27 Nov 2013 on P06 beamline at Petra 3, R. Chernikov and O. Müller, unpublished. 
The source code is in examples\withRaycing\01_SynchrotronSources
Notice that not only the band width is affected by tapering. Also the transverse distribution attains inhomogeneity which varies with energy, see the animation below. Notice also that such a picture is sensitive to emittance; the one below was calculated for the emittance of Petra III ring.
Undulator spectrum in transverse plane¶
The codes calculating undulator spectra – [Urgent], [US], [SPECTRA] – calculate either the spectrum of flux through a given aperture or the transversal distribution at a fixed energy. It is not possible to simultaneously have two dependencies: on energy and on transversal coordinates.
Whereas xrt gives equal results to other codes in such univariate distributions as flux through an aperture:
… and transversal distribution at a fixed energy:
SPECTRA  xrt  

E = 4850 eV (3rd harmonic)  
E = 11350 eV (7th harmonic) 
…, xrt can combine the two distributions in one image and thus be more informative:
zoom in  zoom out  

E ≈ 4850 eV (3rd harmonic)  
E ≈ 11350 eV (7th harmonic) 
In particular, it is seen that divergence strongly depends on energy, even for such a narrow energy band within one harmonic. It is also seen that the maximum flux corresponds to slightly offaxis radiation (greenish color) but not to the onaxis radiation (bluish color).
Undulator radiation in near field¶
Notice that on the following pictures the ppolarized flux is only ~6% of the total flux.
at 5 m  SPECTRA  xrt 

far field at 5 m, full flux  
far field at 5 m, ppolarized  
near field at 5 m, full flux  
near field at 5 m, ppolarized 
at 25 m  SPECTRA  xrt 

far field at 25 m, full flux  
far field at 25 m, ppolarized  
near field at 25 m full flux  
near field at 25 m, ppolarized 
Field phase in near field¶
The phase of the radiation field on a flat screen as calculated by the three codes is compared below. Notice that the visualization (brightness=intensity, color=phase) is not by SRW and SPECTRA but was done by us.
SRW [SRW]  SPECTRA  xrt 

Undulator source size dependent on energy spread and detuning¶
The linear and angular source sizes, as calculated with equations from [TanakaKitamura] (summarized here) for a U19 undulator in MAX IV 3 GeV ring (\(E_1\) = 1429 eV) with \(\varepsilon_x\) = 263 pmrad, \(\varepsilon_y\) = 8 pmrad, \(\beta_x\) = 9 m and \(\beta_y\) = 2 m, are shown below. Energy spread mainly affects the angular sizes and not the linear ones.
The calculated sizes were further compared with those of the sampled field (circles) at variable energies around the nominal harmonic positions, i.e. at so called undulator detuning. To get the photon source size distribution, the angular distributions of Es and Ep field amplitudes were Fourier transformed, as described in [Coïsson]. The sampled field sizes strongly vary due to undulator detuning, as is better seen on the magnified insets. The size variation by detuning is the underlying reason for the size dependence on energy spread: with a nonzero energy spread the undulator becomes effectively detuned for some electrons depending on their velocity.
The size of the circles is proportional to the total flux normalized to the maximum at each respective harmonic. It sharply decreases at the higher energy end of a harmonic and has a long tail at the lower energy end, in accordance with the above examples.
The effect of energy detuning from the nominal undulator harmonic energy on the photon source size is compared to the results by Coïsson [Coïsson] (crosses in the figures below). He calculated the sizes for a single electron field, and thus without emittance and energy spread. For comparison, we also sampled the undulator field at zero energy spread and emittance.
[Coïsson]  (1, 2) R. Coïsson, Effective phase space widths of undulator radiation, Opt. Eng. 27 (1988) 250–2. 
Optical elements¶
Module oes
defines a generic optical element in
class OE
. Its methods serve mainly for propagating the beam downstream
the beamline. This is done in the following sequence: for each ray transform
the beam from global to local coordinate system, find the intersection point
with the OE surface, define the corresponding state of the ray, calculate the
new direction, rotate the coherency matrix to the local sp basis, calculate
reflectivity or transmittivity and apply it to the coherency matrix, rotate the
coherency matrix back and transform the beam from local to the global
coordinate system.
Module oes
defines also several other optical
elements with various geometries.

class
xrt.backends.raycing.oes.
OE
¶ The main base class for an optical element. It implements a generic flat mirror, crystal, multilayer or grating.

__init__
(bl=None, name='', center=[0, 0, 0], pitch=0, roll=0, yaw=0, positionRoll=0, rotationSequence='RzRyRx', extraPitch=0, extraRoll=0, extraYaw=0, extraRotationSequence='RzRyRx', alarmLevel=None, surface=None, material=None, alpha=None, limPhysX=[1000.0, 1000.0], limOptX=None, limPhysY=[1000.0, 1000.0], limOptY=None, isParametric=False, shape='rect', gratingDensity=None, order=None, shouldCheckCenter=False, targetOpenCL=None, precisionOpenCL='float64')¶  bl: instance of
BeamLine
 Container for beamline elements. Optical elements are added to its oes list.
 name: str
 Userspecified name, occasionally used for diagnostics output.
 center: 3sequence of floats
 3D point in global system. Any two coordinates can be ‘auto’ for automatic alignment.
 pitch, roll, yaw: floats
 Rotations Rx, Ry, Rz, correspondingly, defined in the local system. If the material belongs to Crystal, pitch can be calculated automatically if alignment energy is given as a single element list [energy]. If ‘auto’, the alignment energy will be taken from beamLine.alignE.
 positionRoll: float
 A global roll used for putting the OE upside down (=np.pi) or at horizontal deflection (=[]np.pi/2). This parameter does the same rotation as roll. It is introduced for holding large angles, as π or π/2 whereas roll is meant for smaller [mis]alignment angles.
 rotationSequence: str, any combination of ‘Rx’, ‘Ry’ and ‘Rz’
 Gives the sequence of rotations of the OE around the local axes. The sequence is read from left to right (do not consider it as an operator). When rotations are more than one, the final position of the optical element depends on this parameter.
 extraPitch, extraRoll, extraYaw, extraRotationSequence:
 Similar to pitch, roll, yaw, rotationSequence but applied after them. This is sometimes necessary because rotations do not commute. The extra angles were introduced for easier misalignment after the initial positioning of the OE.
 alarmLevel: float or None
 Allowed fraction of incident rays to be absorbed by OE. If exceeded, an alarm output is printed in the console.
 surface: None or sequence of str
 If there are several optical surfaces, such as metalized stripes on a mirror, these are listed here as names; then also the optical limits must all be given by sequences of the same length if not None.
 material: None or sequence of material objects
 The material(s) must have
get_amplitude()
orget_refractive_index()
method. If not None, must correspond to surface. If None, the reflectivities are equal to 1.  alpha: float
 Asymmetry angle for a crystal OE (rad). Positive sign is for the atomic planes’ normal looking towards positive y.
 limPhysX and limPhysY: [min, max] where min, max are
 floats or sequences of floats Physical dimension = local coordinate of the corresponding edge. Can be given by sequences of the length of surface. You do not have to provide the limits, although they may help in finding intersection points, especially for (strongly) curved surfaces.
 limOptX and limOptY: [min, max] where min, max are
 floats or sequences of floats Optical dimension = local coordinate of the corresponding edge. Useful when the optical surface is smaller than the whole surface, e.g. for metalized stripes on a mirror.
 isParametric: bool
 If True, the OE is defined by parametric equations rather than by
z(x, y) function. For example, parametric representation is
useful for describing closed surfaces, such as capillaries. The
user must supply the transformation functions
param_to_xyz()
andxyz_to_param()
between local (x, y, z) and (s, phi, r) and the parametric surface local_r dependent on (s, phi). The exact meaning of these three new parameters is up to the user because this meaning is selfcontained in the above mentioned usersupplied functions. For example, these can be viewed as cylindricallike coordinates, where s is a running coordinate on a 3D axial curve, phi and r are polar coordinates in planes normal to the axial curve and crossing that curve at point s. ClassSurfaceOfRevolution
gives an example of the transformation functions and represents a useful kind of parametric surface. The methodslocal_n()
(surface normal) andlocal_g()
(grating vector, if used for this OE) return 3D vectors in local xyz space but now the two input coordinate parameters are s and phi. The limits [limPhysX, limOptX] and [limPhysY, limOptY] still define, correspondingly, the limits in local x and y. The local beams (footprints) will additionally contain s, phi and r arrays.  shape: str or list of [x, y] pairs
 The shape of OE. Supported: ‘rect’, ‘round’ or a list of [x, y] pairs for an arbitrary shape.
 gratingDensity: None or list
If material kind = ‘grating’, its density can be defined as list [axis, ρ_{0}, P_{0}, P_{1}, …], where ρ_{0} is the constant line density in inverse mm, P_{0} – P_{n} are polynom coefficients defining the line density variation, so that for a given axis
\[\rho_x = \rho_0\cdot(P_0 + 2 P_1 x + 3 P_2 x^2 + ...).\]Example: [‘y’, 800, 1] for the grating with constant spacing along the ‘y’ direction; [‘y’, 1200, 1, 1e6, 3.1e7] for a VLS grating. The length of the list determines the polynomial order.
 order: int or sequence of ints
 The order(s) of grating, FZP or BraggFresnel diffraction.
 shouldCheckCenter: bool
 This is a leagcy parameter designed to work together with alignment
stages – classes in module
stages
– which modify the orientation of an optical element. if True, invokes checkCenter method for checking whether the oe center lies on the original beam line. checkCenter implies vertical deflection and ignores any difference in height. You should override this method for OEs of horizontal deflection.  targetOpenCL: None, str, 2tuple or tuple of 2tuples
 pyopencl can accelerate the search for the intersections of rays with the OE. If pyopencl is used, targetOpenCL is a tuple (iPlatform, iDevice) of indices in the lists cl.get_platforms() and platform.get_devices(), see the section Calculations on GPU. None, if pyopencl is not wanted. Ignored if pyopencl is not installed.
 precisionOpenCL: ‘float32’ or ‘float64’, only for GPU.
 Single precision (float32) should be enough. So far, we do not see any example where double precision is required. The calculations with double precision are much slower. Double precision may be unavailable on your system.
 bl: instance of

local_g
(x, y, rho=100.0)¶ For a grating, gives the local reciprocal groove vector (without 2pi!) in 1/mm at (x, y) position. The vector must lie on the surface, i.e. be orthogonal to the normal. Typically is overridden in the derived classes or defined in Material class. Returns a 3tuple of floats or of arrays of the length of x and y.

local_n
(x, y)¶ Determines the normal vector of OE at (x, y) position. Typically is overridden in the derived classes. If OE is an asymmetric crystal, local_n must return 2 normals as a 6sequence: the 1st one of the atomic planes and the 2nd one of the surface. Note the order!
If isParametric in the constructor is True,
local_n()
still returns 3D vector(s) in local xyz space but now the two input coordinate parameters are s and phi.The result is a 3tuple or a 6tuple. Each element is either a scalar or an array of the length of x and y.

local_n_distorted
(x, y)¶ Distortion to the local normal. If isParametric in the constructor is True, the input arrays are understood as (s, phi).
Distortion can be given in two ways and is signaled by the length of the returned tuple:
 As d_pitch and d_roll rotation angles of the normal (i.e. rotations Rx and Ry). A tuple of the two arrays must be returned. This option is also suitable for parametric coordinates because the two rotations will be around Cartesian axes and the local normal (local_n) is also a 3D vector in local xyz space.
 As a 3D vector that will be added to the local normal calculated at the same coordinates. The returned vector can have any length, not necessarily unity. As for local_n, the 3D vector is in local xyz space even for a parametric surface. The resulted vector local_n + local_n_distorted will be normalized internally before calculating the reflected beam direction. A tuple of 3 arrays must be returned.

local_z
(x, y)¶ Determines the surface of OE at (x, y) position. Typically is overridden in the derived classes. Must return either a scalar or an array of the length of x and y.

multiple_reflect
(beam=None, maxReflections=1000, needElevationMap=False, returnLocalAbsorbed=None)¶ Does the same as
reflect()
but with up to maxReflections reflection on the same surface.The returned beam has additional fields: nRefl for the number of reflections, elevationD for the maximum elevation distance between the rays and the surface as the ray travels between the impact points, elevationX, elevationY, elevationZ for the coordinates of the maximum elevation points.
 returnLocalAbsorbed: None or int
 If not None, returns the absorbed intensity in local beam.

prepare_wave
(prevOE, nrays, shape='auto', area='auto', rw=None)¶ Creates the beam arrays used in wave diffraction calculations. prevOE is the diffracting element: a descendant from
OE
,RectangularAperture
orRoundAperture
. nrays: if int, specifies the number of randomly distributed samples the surface withinself.limPhysX
limits; if 2tuple of ints, specifies (nx, ny) sizes of a uniform mesh of samples.

reflect
(beam=None, needLocal=True, noIntersectionSearch=False, returnLocalAbsorbed=None)¶ Returns the reflected or transmitted beam as \(\vec{out}\) in global and local (if needLocal is true) systems.
Mirror [wikiSnell]:
\[\begin{split}\vec{out}_{\rm reflect} &= \vec{in} + 2\cos{\theta_1}\vec{n}\\ \vec{out}_{\rm refract} &= \frac{n_1}{n_2}\vec{in} + \left(\frac{n_1}{n_2}\cos{\theta_1}  \cos{\theta_2}\right)\vec{n},\end{split}\]where
\[\begin{split}\cos{\theta_1} &= \vec{n}\cdot\vec{in}\\ \cos{\theta_2} &= sign(\cos{\theta_1})\sqrt{1  \left(\frac{n_1}{n_2}\right)^2\left(1\cos^2{\theta_1}\right)}.\end{split}\]Grating or FZP [SpencerMurty]:
For the reciprocal grating vector \(\vec{g}\) and the \(m\)th diffraction order:
\[\vec{out} = \vec{in}  dn\vec{n} + \vec{g}m\lambda\]where
\[dn = \cos{\theta_1} \pm \sqrt{\cos^2{\theta_1}  2(\vec{g}\cdot\vec{in})m\lambda  \vec{g}^2 m^2\lambda^2}\]Crystal [SanchezDelRioCerrina]:
Crystal (generally asymmetrically cut) is considered a grating with the reciprocal grating vector equal to
\[\vec{g} = \left(\vec{n_\vec{H}}  (\vec{n_\vec{H}}\cdot\vec{n})\vec{n})\right) / d_\vec{H}.\]Note that \(\vec{g}\) is along the line of the intersection of the crystal surface with the plane formed by the two normals \(\vec{n_\vec{H}}\) and \(\vec{n}\) and its length is \(\vec{g}=\sin{\alpha}/d_\vec{H}\), with \(\alpha\) being the asymmetry angle.
[wikiSnell] http://en.wikipedia.org/wiki/Snell%27s_law . [SpencerMurty] G. H. Spencer and M. V. R. K. Murty, J. Opt. Soc. Am. 52 (1962) 672. [SanchezDelRioCerrina] M. Sánchez del Río and F. Cerrina, Rev. Sci. Instrum. 63 (1992) 936.  returnLocalAbsorbed: None or int
 If not None, returns the absorbed intensity in local beam.
 noIntersectionSearch: bool
 Used in wave propagation, normally should be False.


class
xrt.backends.raycing.oes.
DicedOE
(OE)¶ Base class for a diced optical element. It implements a flat diced mirror.

__init__
(*args, **kwargs)¶  dxFacet, dyFacet: float
 Size of the facets.
 dxGap, dyGat: float
 Width of the gap between facets.

facet_center_n
(x, y)¶ Surface normal or (Bragg normal and surface normal).

facet_center_z
(x, y)¶ Z of the facet centers at (x, y).

facet_delta_n
(u, v)¶ Local surface normal (always without Bragg normal!) in the facet coordinates. In the asymmetry case the lattice normal is taken as constant over the facet and is given by
facet_center_n()
.

facet_delta_z
(u, v)¶ Local Z in the facet coordinates.


class
xrt.backends.raycing.oes.
JohannCylinder
(OE)¶ Simply bent reflective crystal.

__init__
(*args, **kwargs)¶  Rm: float
 Meridional radius.
 crossSection: str
 Determines the bending shape: either ‘circular’ or ‘parabolic’.


class
xrt.backends.raycing.oes.
JohanssonCylinder
(JohannCylinder)¶ Groundbent (Johansson) reflective crystal.

class
xrt.backends.raycing.oes.
JohannToroid
(OE)¶ 2D bent reflective crystal.

__init__
(*args, **kwargs)¶  Rm and Rs: float
 Meridional and sagittal radii.


class
xrt.backends.raycing.oes.
JohanssonToroid
(JohannToroid)¶ Ground2Dbent (Johansson) reflective optical element.

class
xrt.backends.raycing.oes.
GeneralBraggToroid
(JohannToroid)¶ Ground2Dbent reflective optical element with 4 independent radii: meridional and sagittal for the surface (Rm and Rs) and the atomic planes (RmBragg and RsBragg).

class
xrt.backends.raycing.oes.
DicedJohannToroid
(DicedOE, JohannToroid)¶ A diced version of
JohannToroid
.

class
xrt.backends.raycing.oes.
DicedJohanssonToroid
(DicedJohannToroid, JohanssonToroid)¶ A diced version of
JohanssonToroid
.

class
xrt.backends.raycing.oes.
LauePlate
(OE)¶ A flat Laue plate. The thickness is defined in its material part.

class
xrt.backends.raycing.oes.
BentLaueCylinder
(OE)¶ Simply bent reflective optical element in Laue geometry (duMond).

__init__
(*args, **kwargs)¶  R: float or 2tuple.
 Meridional radius. Can be given as (p, q) for automatic calculation based the “Coddington” equations.
 crossSection: str
 Determines the bending shape: either ‘circular’ or ‘parabolic’.


class
xrt.backends.raycing.oes.
GroundBentLaueCylinder
(BentLaueCylinder)¶ Groundbent reflective optical element in Laue geometry.

class
xrt.backends.raycing.oes.
BentLaueSphere
(BentLaueCylinder)¶ Spherically bent reflective optical element in Laue geometry.

class
xrt.backends.raycing.oes.
BentFlatMirror
(OE)¶ Implements cylindrical parabolic mirror. Exemplifies inclusion of a new parameter (here, R) without the need of explicit repetition of all the parameters of the parent class.

class
xrt.backends.raycing.oes.
ToroidMirror
(OE)¶ Implements toroidal mirror. Exemplifies inclusion of new parameters (here, R and r) without the need of explicit repetition of all the parameters of the parent class.

class
xrt.backends.raycing.oes.
EllipticalMirrorParam
(OE)¶ The elliptical mirror is implemented as a parametric surface. The parameterization is the following: s  is local coordinate along the major axis with origin at the ellipse center. phi and r are local polar coordinates in planes normal to the major axis at every point s. The polar axis is upwards.
The user supplies two foci either by focal distances p and q (both are positive) or as f1 and f2 points in the global coordinate system (3sequences). Any combination of (p or f1) and (q or f2) is allowed. If p is supplied, not f1, the incoming optical axis is assumed to be along the global Y axis. For a general orientation of the ellipse axes f1 or pAxis – the p arm direction in global coordinates – should be supplied.
If isCylindrical is True, the figure is an elliptical cylinder, otherwise it is an ellipsoid of revolution around the major axis.
Warning
If you want to change any of p, q, f1, f2 or pAxis after the creation of the OE, you must invoke the method
reset_pq()
to recalculate the ellipsoid parameters.The usage is exemplified in test_param_mirror.py.

class
xrt.backends.raycing.oes.
ParabolicalMirrorParam
(EllipticalMirrorParam)¶ The parabolical mirror is implemented as a parametric surface. The parameterization is the following: s  is local coordinate along the paraboloid axis with origin at the focus. phi and r are local polar coordinates in planes normal to the axis at every point s. The polar axis is upwards.
The user supplies one (and only one) focal distance p or q as a positive value. Alternatively, instead of p one can specify f1 (3sequence) as a 3D point in the global coordinate system and instead of q – f2. If p or q is supplied, the paraboloid axis isassumed to be along the global Y axis, otherwise supply parabolaAxis as a vector in global coordinates.
If isCylindrical is True, the figure is an parabolical cylinder, otherwise it is a paraboloid of revolution around the major axis.
Warning
If you want to change any of p, q, f1, f2 or parabolaAxis after the creation of the OE, you must invoke the method
reset_pq()
to recalculate the paraboloid parameters.The usage is exemplified in test_param_mirror.py.

class
xrt.backends.raycing.oes.
DCM
(OE)¶ Implements a Double Crystal Monochromator with flat crystals.

__init__
(*args, **kwargs)¶  bragg: float, str, list
 Bragg angle in rad. Can be calculated automatically if alignment energy is given as a single element list [energy]. If ‘auto’, the alignment energy will be taken from beamLine.alignE.
 cryst1roll, cryst2roll, cryst2pitch, cryst2finePitch: float
 Misalignment angles in rad.
 cryst2perpTransl, cryst2longTransl: float
 perpendicular and longitudinal translations of the 2nd crystal in respect to the 1st one.
 limPhysX2, limPhysY2, limOptX2, limOptY2, material2:
 refer to the 2nd crystal and are similar to the same parameters
of the parent class
OE
without the trailing “2”.  fixedOffset: float
 Offset between the incoming and outcoming beams in mm. If not None or zero the value of cryst2perpTransl is replaced by fixedOffset/2/cos(bragg)

double_reflect
(beam=None, needLocal=True, fromVacuum1=True, fromVacuum2=True, returnLocalAbsorbed=None)¶ Returns the reflected beam in global and two local (if needLocal is true) systems.
 returnLocalAbsorbed: None or int
 If not None, returns the absorbed intensity in local beam. If equals zero, total absorbed intensity is return in the last local beam, otherwise the Nth local beam returns the absorbed intensity on Nth surface of the optical element.


class
xrt.backends.raycing.oes.
DCMwithSagittalFocusing
(DCM)¶ Creates a DCM with a horizontally focusing 2nd crystal.

__init__
(*args, **kwargs)¶ Assume Bragg planes and physical surface planes are parallel (no miscut angle).
 Rs: float
 Sagittal radius of second crystal.


class
xrt.backends.raycing.oes.
Plate
(DCM)¶ Implements a body with two surfaces. Is derived from
DCM
because it also has two interfaces but the parameters referring to the 2nd crystal should be ignored.
__init__
(*args, **kwargs)¶  t: float
 Tthickness in mm.

double_refract
(beam=None, needLocal=True, returnLocalAbsorbed=None)¶ Returns the refracted beam in global and two local (if needLocal is true) systems.
 returnLocalAbsorbed: None, int
 If not None, returned local beam represents the absorbed intensity instead of transmitted. The parameter defines the ordinal number of the surface in multisurface element to return the absorbed intensity, i.e. 1 for the entrance surface of the plate, 2 for the exit, 0 for total intensity absorbed in the element.


class
xrt.backends.raycing.oes.
ParaboloidFlatLens
(Plate)¶ Implements a refractive lens or a stack of lenses (CRL) with one side as paraboloid and the other one flat.

__init__
(*args, **kwargs)¶  focus: float
The focal distance of the of paraboloid in mm. The paraboloid is then defined by the equation:
\[z = (x^2 + y^2) / (4 * \mathit{focus})\]Note
This is not the focal distance of the lens but of the parabola! The former also depends on the refractive index. focus is only a shape parameter!
 pitch: float
 the default value is set to π/2, i.e. to normal incidence.
 zmax: float
 If given, limits the z coordinate; the object becomes then a plate of the thickness zmax + t with a paraboloid hole at the origin.
 nCRL: int or tuple (focalDistance, E)
 If used as CRL (a stack of several lenslets), the number of the
lenslets nCRL is either given by the user directly or calculated
for focalDistance at energy E and then rounded. The lenses are
stacked along the local [0, 0, 1] direction with the step equal to
zmax + t for curvedflat lenses or 2*zmax + t for
double curved lenses. For propagation with nCRL > 1 please use
multiple_refract()
.


class
xrt.backends.raycing.oes.
ParabolicCylinderFlatLens
(ParaboloidFlatLens)¶ Implements a refractive lens or a stack of lenses (CRL) with one side as parabolic cylinder and the other one flat. If used as a CRL, the lenslets are arranged such that they alternatively focus in the 45° and +45° planes. Therefore the total number of lenslets is doubled as compared to ParaboloidFlatLens case.

class
xrt.backends.raycing.oes.
DoubleParaboloidLens
(ParaboloidFlatLens)¶ Implements a refractive lens or a stack of lenses (CRL) with two equal paraboloids from both sides.

class
xrt.backends.raycing.oes.
DoubleParabolicCylinderLens
(ParabolicCylinderFlatLens)¶ Implements a refractive lens or a stack of lenses (CRL) with two equal parabolic cylinders from both sides.

class
xrt.backends.raycing.oes.
SurfaceOfRevolution
(OE)¶ Base class for parametric surfaces of revolution. The parameterization implements cylindrical coordinates, where s is y (along the beamline), and phi and r are polar coordinates in planes normal to s.

class
xrt.backends.raycing.oes.
NormalFZP
(OE)¶ Implements a circular Fresnel Zone Plate, as it is described in XRay Data Booklet, Section 4.4.
Warning
Do not forget to specify
kind='FZP'
in the material!
__init__
(*args, **kwargs)¶  f: float
 The focal distance (mm) calculated for waves of energy E.
 E: float
 Energy (eV) for which f is calculated.
 N: int
 The number of zones. Is either directly given or calculated from thinnestZone (mm).
 thinnestZone: float
 In mm; can be given to calculate N.
 isCentralZoneBlack: bool
 if False, the zones are inverted.
 order: int or sequence of ints
 Needed diffraction order(s).

rays_good
(x, y, z, is2ndXtal=False)¶ Returns state value for a ray with the given intersection point (x, y) with the surface of OE: 1: good (intersected) 2: reflected outside of working area (“out”), 3: transmitted without intersection (“over”), NN: lost (absorbed) at OE#NN  OE numbering starts from 1 !!!
Note, x, y, z are local Cartesian coordinates, even for a parametric OE.


class
xrt.backends.raycing.oes.
GeneralFZPin0YZ
(OE)¶ Implements a general Fresnel Zone Plate, where the zones are determined by two foci and the surface shape of the OE.
Warning
Do not forget to specify
kind='FZP'
in the material!
__init__
(*args, **kwargs)¶  f1 and f2: 3 or 4sequence or str
 The two foci given by 3sequences representing 3D points in _local_ coordinates or ‘inf’ for infinite position. The 4th member in the sequence can be 1 to give the negative sign to the path if both foci are on the same side of the FZP.
 E: float
 Energy (eV) for which f is calculated.
 N: int
 The number of zones.
 grazingAngle: float
 The angle of the main optical axis to the surface. Defaults to self.pitch.
 phaseShift: float
 The zones can be phase shifted, which affects the zone structure but does not affect the focusing. if phaseShift is 0, the central zone is at the constructive interference.
 order: int or sequence of ints
 Needed diffraction order(s).


class
xrt.backends.raycing.oes.
BlazedGrating
(OE)¶ Implements a grating of triangular shape given by two angles. The front side of the triangle (the one looking towards the source) is at blaze angle to the base plane. The back side is at antiblaze angle.
Note
In contrast to the geometric implementation of the grating diffraction when the deflection is calculated by the grating equation, the diffraction by
BlazedGrating
is meant to be used by the wave propagation methods, see Gallery of plots and scripts 3. Wave propagation. In those methods, the diffraction is not given by the grating equation but by the surface itself through the calculation of the Kirchhoff integral. Therefore the surface material should not have the propertykind='grating'
but ratherkind='mirror'
.A usual optical element (of class
OE
) with such a developed surface would have troubles in finding correct intersection points because for each ray there are several solutions and we implicitly assume only one. The classBlazedGrating
implements an ad hoc methodfind_intersection()
for selecting the 1st intersection point among the several possible ones. The left picture below illustrates the behavior ofOE
(the footprint shown by circles is partially in the shadowed area). The right picture demonstrates the correct behavior ofBlazedGrating
in respect to illumination and shadowing. Notice that wave propagation gives the same result for the two cases, apart from a small vertical shift. The difference is purely esthetic.

class
xrt.backends.raycing.oes.
LaminarGrating
(OE)¶ Implements a grating of rectangular profile.

class
xrt.backends.raycing.oes.
VLSLaminarGrating
(OE)¶ Implements a grating of rectangular profile with variable period.
Distorted surfaces¶
For introducing an error to an ideal surface you must define two methods in
your descendant of the OE
: local_z_distorted
(or
local_r_distorted
for a parametric surface) and local_n_distorted
. The
latter method returns two angles d_pitch and d_roll or a 3D vector that will be
added to the local normal. See the docstrings of OE.local_n_distorted()
and the example ‘Defocusing by a distorted mirror’.
Materials¶
Module materials
defines atomic and material
properties related to xray scattering, diffraction and propagation:
reflectivity, transmittivity, refractive index, absorption coefficient etc.

xrt.backends.raycing.materials.
read_atomic_data
(elem)¶ Reads atomic data from
AtomicData.dat
file adopted from XOP [XOP]. It has the following data: 0 AtomicRadius[Å] CovalentRadius[Å] AtomicMass BoilingPoint[K] MeltingPoint[K] Density[g/ccm] AtomicVolume CoherentScatteringLength[1E12cm] IncoherentXsection[barn] Absorption@1.8Å[barn] DebyeTemperature[K] ThermalConductivity[W/cmK]In
read_atomic_data()
only the mass is inquired. The user may extend the method to get the other values by simply adding the corresponding array elements to the returned value.

class
xrt.backends.raycing.materials.
Element
¶ This class serves for accessing the scattering factors f0, f1 and f2 of a chemical element. It can also report other atomic data listed in
AtomicData.dat
file adopted from XOP [XOP].
__init__
(elem=None, table='Chantler')¶  elem: str or int
 The element can be specified by its name (case sensitive) or its ordinal number.
 table: str
 This parameter is explained in the description of
Material
.

get_f0
(qOver4pi=0)¶ Calculates f0 for the given qOver4pi.

get_f1f2
(E)¶ Calculates (interpolates) f1 and f2 for the given array E.

read_f0_Kissel
()¶ Reads f0 scattering factors from the tabulation of XOP [XOP]. These were calculated by [Kissel] and then parameterized as [Waasmaier]:
\[f_0\left(\frac{q}{4\pi}\right) = c + \sum_{i=1}^5{a_i\exp\left(b_i \left(q/(4\pi)\right)^2\right)}\]where \(q/(4\pi) = \sin{\theta} / \lambda\) and \(a_i\), \(b_i\) and \(c\) are the coefficients tabulated in the file
f0_xop.dat
.[Kissel] L. Kissel, Radiation physics and chemistry 59 (2000) 185200, http://wwwphys.llnl.gov/Research/scattering/RTAB.html [Waasmaier] D. Waasmaier & A. Kirfel, Acta Cryst. A51 (1995) 416413

read_f1f2_vs_E
(table)¶ Reads f1 and f2 scattering factors from the given table at the instantiation time.


class
xrt.backends.raycing.materials.
Material
¶ Material
serves for getting reflectivity, transmittivity, refractive index and absorption coefficient of a material specified by its chemical formula and density.
__init__
(elements=None, quantities=None, kind='auto', rho=0, t=None, table='Chantler total', efficiency=None, efficiencyFile=None, name='')¶  elements: str or sequence of str
 Contains all the constituent elements (symbols)
 quantities: None or sequence of floats of length of elements
 Coefficients in the chemical formula. If None, the coefficients are all equal to 1.
 kind: str
 One of ‘mirror’, ‘thin mirror’, ‘plate’, ‘lens’, ‘grating’, ‘FZP’.
If ‘auto’, the optical element will decide which material kind to
use via its method
assign_auto_material_kind()
.  rho: float
 Density in g/cm^{3}.
 t: float
 Thickness in mm, required only for ‘thin mirror’.
 table: str
At the time of instantiation the tabulated scattering factors of each element are read and then interpolated at the requested q value and energy. table can be ‘Henke’ (10 eV < E < 30 keV) [Henke], ‘Chantler’ (11 eV < E < 405 keV) [Chantler] or ‘BrCo’ (30 eV < E < 509 keV) [BrCo].
The tables of f2 factors consider only photoelectric crosssections. The tabulation by Chantler can optionally have total absorption crosssections. This option is enabled by table = ‘Chantler total’.
[Henke] http://henke.lbl.gov/optical_constants/asf.html B.L. Henke, E.M. Gullikson, and J.C. Davis, Xray interactions: photoabsorption, scattering, transmission, and reflection at E=5030000 eV, Z=192, Atomic Data and Nuclear Data Tables 54 (no.2) (1993) 181342. [Chantler] http://physics.nist.gov/PhysRefData/FFast/Text/cover.html http://physics.nist.gov/PhysRefData/FFast/html/form.html C. T. Chantler, Theoretical Form Factor, Attenuation, and Scattering Tabulation for Z = 1  92 from E = 1  10 eV to E = 0.4  1.0 MeV, J. Phys. Chem. Ref. Data 24 (1995) 71643. [BrCo] http://www.bmsc.washington.edu/scatter/periodictable.html ftp://ftpa.aps.anl.gov/pub/crosssection_codes/ S. Brennan and P.L. Cowan, A suite of programs for calculating xray absorption, reflection and diffraction performance for a variety of materials at arbitrary wavelengths, Rev. Sci. Instrum. 63 (1992) 850853.  efficiency: sequence of pairs [order, value]
 Can be given for kind = ‘grating’ and kind = ‘FZP’. It must
correspond to the field order of the OE. It can be given as a
constant per diffraction order or as an energy dependence, also per
diffraction order. It is a sequence of pairs [order, value],
where value is either the efficiency itself or an index in the
data file given by efficiencyFile. The data file can either be
(1) a pickle file with energy and efficiency arrays as two
first dump elements and efficiency shape as (len(energy),
orders) or (2) a column file with energy in the leftmost column
and the order efficiencies in the next columns. The value is a
corresponding array index (zerobased) or a column number (also
zerobased, the 0th column is energy). An example of the efficiency
calculation can be found in
\examples\withRaycing\11_Wave\waveGrating.py
.  efficiencyFile: str
 See the definition of efficiency.
 name: str
 Material name. Not used by xrt. Can be used by the user for annotations of graphs or other output purposes. If empty, the name is constructed from the elements and the quantities.

get_absorption_coefficient
(E)¶ Calculates the linear absorption coefficient from the imaginary part of refractive index. E can be an array. The result is in cm^{1}.
\[\mu = 2 \Im(n) k.\]

get_amplitude
(E, beamInDotNormal, fromVacuum=True)¶ Calculates amplitude of reflectivity (for ‘mirror’ and ‘thin mirror’) or transmittivity (for ‘plate’ and ‘lens’) [wikiFresnelEq], [AlsNielsen]. E is energy, beamInDotNormal is cosine of the angle between the incoming beam and the normal (\(\theta_1\) below), both can be scalars or arrays. The interface of the material is assumed to be with vacuum; the direction is given by boolean fromVacuum. Returns a tuple of the amplitudes of s and p polarizations and the absorption coefficient in cm^{1}.
\[\begin{split}r_s^{\rm mirror} &= \frac{n_1\cos{\theta_1}  n_2\cos{\theta_2}} {n_1\cos{\theta_1} + n_2\cos{\theta_2}}\\ r_p^{\rm mirror} &= \frac{n_2\cos{\theta_1}  n_1\cos{\theta_2}} {n_2\cos{\theta_1} + n_1\cos{\theta_2}}\\ r_{s,p}^{\rm thin\ mirror} &= r_{s,p}^{\rm mirror}\frac{1  p^2} {1  (r_{s,p}^{\rm mirror})^2p^2},\end{split}\]where the phase factor \(p^2 = \exp(2iEtn_2\cos{\theta_2}/c\hbar)\).
\[\begin{split}t_s^{\rm plate,\ lens} &= 2\frac{n_1\cos{\theta_1}} {n_1\cos{\theta_1} + n_2\cos{\theta_2}}t_f\\ t_p^{\rm plate,\ lens} &= 2\frac{n_1\cos{\theta_1}} {n_2\cos{\theta_1} + n_1\cos{\theta_2}}t_f\\\end{split}\]where \(t_f = \sqrt{\frac{\Re(n_2n_1)\cos{\theta_2}} {cos{\theta_1}}}/n_1\).
[wikiFresnelEq] http://en.wikipedia.org/wiki/Fresnel_equations . [AlsNielsen] (1, 2) Jens AlsNielsen, Des McMorrow, Elements of Modern Xray Physics, John Wiley and Sons, 2001.

get_refractive_index
(E)¶ Calculates refractive index at given E. E can be an array.
\[n = 1  \frac{r_0\lambda^2 N_A \rho}{2\pi M}\sum_i{x_i f_i(0)}\]where \(r_0\) is the classical electron radius, \(\lambda\) is the wavelength, \(N_A\) is Avogadro’s number, \(\rho\) is the material density, M is molar mass, \(x_i\) are atomic concentrations (coefficients in the chemical formula) and \(f_i(0)\) are the complex atomic scattering factor for the forward scattering.


class
xrt.backends.raycing.materials.
Multilayer
¶ Multilayer
serves for getting reflectivity of a multilayer. The multilayer may have variable thicknesses of the two alternating layers as functions of local x and y and/or as a function of the layer number.
__init__
(tLayer=None, tThickness=0.0, bLayer=None, bThickness=0.0, nPairs=0.0, substrate=None, tThicknessLow=0.0, bThicknessLow=0.0, idThickness=0.0, power=2.0, substRoughness=0, substThickness=inf, name='', geom='reflected')¶  tLayer, bLayer, substrate: instance of
Material
 The top layer material, the bottom layer material and the substrate material.
 tThickness and bThickness: float in Å
 The thicknesses of the layers. If the multilayer is depth graded, tThickness and bThickness are at the top and tThicknessLow and bThicknessLow are at the substrate. If you need laterally graded thicknesses, modify get_t_thickness and/or get_b_thickness in a subclass.
 power: float
Defines the exponent of the layer thickness power law, if the multilayer is depth graded:
\[d_n = A / (B + n)^{power}.\] tThicknessLow and bThicknessLow: float
 Are ignored (left as zeros) if not depth graded.
 nPairs: int
 The number of layer pairs.
 idThickness: float in Å
 RMS thickness \(\sigma_{j,j1}\) of the interdiffusion/roughness interface.
 substThickness: float in Å
 Is only relevant in transmission if substrate is present.
 geom: str
 Either ‘transmitted’ or ‘reflected’.
 tLayer, bLayer, substrate: instance of

get_amplitude
(E, beamInDotNormal, x=None, y=None, ucl=None)¶ Calculates amplitude of reflectivity [AlsNielsen]. E is energy, beamInDotNormal is cosine of the angle between the incoming beam and the normal (\(\theta_0\) below), both can be scalars or arrays. The top interface of the multilayer is assumed to be with vacuum. Returns a tuple of the amplitudes of s and p polarizations.
The calculation starts from the bottommost layer (with index \(N\)). The reflectivity from its top into the adjacent layer (\(N1\)) is:
\[R_N = \frac{r_{N1, N} + r_{N, N+1} p_N^2} {1 + r_{N1, N} r_{N, N+1} p_N^2},\]where the capital \(R\) denotes the net reflectivity of the layer and the small letters \(r\) denote the interface reflectivity (Fresnel equations):
\[r_{j, j+1} = \frac{Q_j  Q_{j+1}}{Q_j + Q_{j+1}},\]here \(N+1\) refers to the substrate material and
\[Q_j = \sqrt{Q^2  8k^2\delta_j + i8k^2\beta_j}, \quad Q = 2k\sin{\theta_0}\]and \(\delta_j\) and \(\beta_j\) are parts of the refractive index \(n_j = 1  \delta_j + i\beta_j\). The phase factor \(p_j^2\) is \(\exp(i\Delta_j Q_j)\), \(\Delta_j\) being the layer thickness. The calculation proceeds recursively upwards by layers as
\[R_j = \frac{r_{j1, j} + R_{j+1} p_j^2} {1 + r_{j1, j} R_{j+1} p_j^2},\]until \(R_1\) is reached, where the 0th layer is vacuum and \(Q_0 = Q\).
If the interdiffusion thickness is not zero, the reflectivity at each interface is attenuated by a factor of \(exp(2k_{j,z}k_{j1,z}\sigma^{2}_{j,j1})\), where \(k_{j,z}\) is longitudinal component of the wave vector in jth layer [NevotCroce].
The above formulas refer to s polarization. The p part differs at the interface:
\[r^p_{j, j+1} = \frac{Q_j\frac{n_{j+1}}{n_j}  Q_{j+1}\frac{n_{j}}{n_{j+1}}}{Q_j\frac{n_{j+1}}{n_j} + Q_{j+1}\frac{n_{j}}{n_{j+1}}}\]and thus the p polarization part requires a separate recursive chain.
In transmission, the recursion is the following:
\[T_{N+1} = \frac{t_{N, N+1}t_{N+1, N+2}p_{N+1}} {1 + r_{N, N+1} r_{N+1, N+2} p_N^2}, \quad T_j = \frac{T_{j+1}t_{j1, j}p_j}{1 + r_{j1, j} R_{j+1} p_j^2},\]where the layer \(N+2\) is vacuum and the interface transmittivities for the two polarizations are equal to:
\[t^s_{j, j+1} = \frac{2Q_j}{Q_j + Q_{j+1}}, \quad t^p_{j, j+1} = \frac{2Q_j\frac{n_{j+1}}{n_j}} {Q_j\frac{n_{j+1}}{n_j} + Q_{j+1}\frac{n_{j}}{n_{j+1}}}\][NevotCroce] L. Nevot and P. Croce, Rev. Phys. Appl. 15, (1980) 761

get_b_thickness
(x, y, iPair)¶ The bottom (the lower in the period pair) layer thickness in Å as a function of local coordinates x and y and the index (zero at vacuum) of the period pair.
For parametric surfaces, the x and y local coordinates are assumed to be s and phi of the parametric representation.

get_dtheta_symmetric_Bragg
(E, order=1)¶ The angle correction for the symmetric Bragg case:
\[\delta\theta = \theta_B  \arcsin(\sqrt{m^2\lambda^2 + 8 d^2 \overline\delta} / 2d),\]where \(\overline\delta\) is the periodaveraged real part of the refractive index.

get_t_thickness
(x, y, iPair)¶ The top (the upper in the period pair) layer thickness in Å as a function of local coordinates x and y and the index (zero at vacuum) of the period pair.
For parametric surfaces, the x and y local coordinates are assumed to be s and phi of the parametric representation.


class
xrt.backends.raycing.materials.
Coated
¶ Derivative class from
Mutilayer
with a single reflective layer on a substrate.
__init__
(*args, **kwargs)¶  coating, substrate: instance of
Material
 Material of the mirror coating layer, and the substrate material.
 cThickness: float
 The thicknesses of mirror coating in Å.
 surfaceRoughness: float
 RMS rougness of the mirror surface in Å.
 substRoughness: float
 RMS rougness of the mirror substrate in Å.
 coating, substrate: instance of


class
xrt.backends.raycing.materials.
Crystal
(Material)¶ The parent class for crystals. The descendants must define
get_structure_factor()
.Crystal
gives reflectivity and transmittivity of a crystal in Bragg and Laue cases.
__init__
(hkl=[1, 1, 1], d=0, V=None, elements='Si', quantities=None, rho=0, t=None, factDW=1.0, geom='Bragg reflected', table='Chantler', name='', nuPoisson=0.0, calcBorrmann=None, useTT=False, mosaicity=0)¶  hkl: sequence
 hkl indices.
 d: float
 Interatomic spacing in Å.
 V: float
 Unit cell volume in Å^{3}. If not given, is calculated from d assuming a cubic symmetry.
 factDW: float
 DebyeWaller factor applied to the structure factor.
 geom: str
 The 1st word is either ‘Bragg’ or ‘Laue’, the 2nd word is either ‘transmitted’ or ‘reflected’ or ‘Fresnel’ (the optical element must then provide local_g method that gives the grating vector).
 table: str
 This parameter is explained in the description of the parent class
Material
.  nuPoisson: float
 Poisson’s ratio. Used to calculate the properties of bent crystals.
 calcBorrmann: str
 Controls the origin of the ray leaving the crystal. Can be ‘None’,
‘uniform’, ‘Bessel’ or ‘TT’. If ‘None’, the point of reflection
is located on the surface of incidence. In all other cases the
coordinate of the exit point is sampled according to the
corresponding distribution: ‘uniform’ is a fast approximation for
thick crystals, ‘Bessel’ is exact solution for the flat crystals,
‘TT’ is exact solution of TakagiTaupin equations for bent and flat
crystals (‘TT’ requires targetOpenCL in the Optical Element to be
not ‘None’ and useTT in the
Crystal
to be ‘True’. Not recommended for crystals thicker than 100 µm due to heavy computational load).  useTT: bool
 Specifies whether the reflectivity will by calculated by analytical formula or by solution of the TakagiTaupin equations (so far only for the Laue geometry). Must be set to ‘True’ in order to calculate the reflectivity of bent crystals.
 mosaicity: float, radians
The sigma of the normal distribution of the crystallite normals.
xrt follows the concept of mosaic crystals from [SanchezDelRioMosaic]. This concept has three main parts: (i) a random distribution of the crystallite normals results in a distribution in the reflected directions, (ii) the secondary extinction results in a mean free path distribution in the new ray origins and (iii) the reflectivity is calculated following the work [BaconLowde].
Note
The mosaicity is assumed large compared with the Darwin width. Therefore, there is no continuous transition mosaictoperfect crystal at a continuously reduced mosaicity parameter.
See the tests here.
[SanchezDelRioMosaic] M. Sánchez del Río et al., Rev. Sci. Instrum. 63 (1992) 932. [BaconLowde] G. E. Bacon and R. D. Lowde, Acta Crystallogr. 1, (1948) 303.

get_Darwin_width
(E, b=1.0, polarization='s')¶ Calculates the Darwin width as
\[2\delta = C\sqrt{\chi_h\chi_{\overline{h}} / b}/\sin{2\theta}\]

get_amplitude
(E, beamInDotNormal, beamOutDotNormal=None, beamInDotHNormal=None, xd=None, yd=None)¶ Calculates complex amplitude reflectivity and transmittivity for s and ppolarizations (\(\gamma = s, p\)) in Bragg and Laue cases for the crystal of thickness L, based upon Belyakov & Dmitrienko [BD]:
\[\begin{split}R_{\gamma}^{\rm Bragg} &= \chi_{\vec{H}}C_{\gamma}\left(\alpha + i\Delta_{\gamma}\cot{l_{\gamma}}\right)^{1}b^{\frac{1}{2}}\\ T_{\gamma}^{\rm Bragg} &= \left(\cos{l{_\gamma}}  i\alpha\Delta {_\gamma}^{1}\sin{l_{\gamma}}\right)^{1} \exp{\left(i\vec{\kappa}_0^2 L (\chi_0  \alpha b) (2\vec{\kappa}_0\vec{s})^{1}\right)}\\ R_{\gamma}^{\rm Laue} &= \chi_{\vec{H}}C_{\gamma} \Delta_{\gamma}^{1}\sin{l_{\gamma}}\exp{\left(i\vec{\kappa}_0^2 L (\chi_0  \alpha b) (2\vec{\kappa}_0\vec{s})^{1}\right)} b^{\frac{1}{2}}\\ T_{\gamma}^{\rm Laue} &= \left(\cos{l_{\gamma}} + i\alpha \Delta_{\gamma}^{1}\sin{l_{\gamma}}\right) \exp{\left(i\vec{\kappa}_0^2 L (\chi_0  \alpha b) (2\vec{\kappa}_0\vec{s})^{1}\right)}\end{split}\]where
\[\begin{split}\alpha &= \frac{\vec{H}^2 + 2\vec{\kappa}_0\vec{H}} {2\vec{\kappa}_0^2}+\frac{\chi_0(1b)}{2b}\\ \Delta_{\gamma} &= \left(\alpha^2 +\frac{C_{\gamma}^2\chi_{\vec{H}} \chi_{\overline{\vec{H}}}}{b}\right)^{\frac{1}{2}}\\ l_{\gamma} &= \frac{\Delta_{\gamma}\vec{\kappa}_0^2L} {2\vec{\kappa}_{\vec{H}}\vec{s}}\\ b &= \frac{\vec{\kappa}_0\vec{s}}{\vec{\kappa}_{\vec{H}}\vec{s}}\\ C_s &= 1, \quad C_p = \cos{2\theta_B}\end{split}\]In the case of thick crystal in Bragg geometry:
\[R_{\gamma}^{\rm Bragg} = \frac{\chi_{\vec{H}} C_{\gamma}} {\alpha\pm\Delta_{\gamma}}b^{\frac{1}{2}}\]with the sign in the denominator that gives the smaller modulus of \(R_\gamma\).
\(\chi_{\vec{H}}\) is the Fourier harmonic of the xray susceptibility, and \(\vec{H}\) is the reciprocal lattice vector of the crystal. \(\vec{\kappa}_0\) and \(\vec{\kappa}_{\vec{H}}\) are the wave vectors of the direct and diffracted waves. \(\chi_{\vec{H}}\) is calculated as:
\[\chi_{\vec{H}} =  \frac{r_0\lambda^2}{\pi V}F_{\vec{H}},\]where \(r_e = e^2 / mc^2\) is the classical radius of the electron, \(\lambda\) is the wavelength, V is the volume of the unit cell.
Notice \(b^{\frac{1}{2}}\) added to the formulas of Belyakov & Dmitrienko in the cases of Bragg and Laue reflections. This is needed because ray tracing deals not with wave fields but with rays and therefore not with intensities (i.e. per crosssection) but with flux.
[BD] V. A. Belyakov and V. E. Dmitrienko, Polarization phenomena in xray optics, Uspekhi Fiz. Nauk. 158 (1989) 679–721, Sov. Phys. Usp. 32 (1989) 697–719. xd and yd are local coordinates of the corresponding optical element. If they are not None and crystal’s get_d method exists, the d spacing is given by the get_d method, otherwise it equals to self.d. In a parametric representation, xd and yd are the same parametric coordinates used in local_r and local_n` methods of the corresponding optical element.

get_dtheta
(E, alpha=None)¶ The angle correction for the general asymmetric case:
\[\begin{split}\delta\theta = \frac{\mp \gamma_0 \pm \sqrt{\gamma_0^2 \mp (\gamma_0  \gamma_h) \sqrt{1  \gamma_0^2} \chi_0 / \sin{2\theta_B}}}{\sqrt{1  \gamma_0^2}}\\\end{split}\]where \(\gamma_0 = \sin(\theta_B + \alpha)\), \(\gamma_h = \mp \sin(\theta_B  \alpha)\) and the upper sign is for Bragg and the lower sign is for Laue geometry.

get_dtheta_regular
(E, alpha=None)¶ The angle correction for the general asymmetric case in its simpler version:
\[\begin{split}\delta\theta = (1  b)/2 \cdot \chi_0 / \sin{2\theta_B}\\ b = \sin(\theta_B + \alpha) / \sin(\theta_B  \alpha)\end{split}\]For the symmetric Bragg b = 1 and for the symmetric Laue b = +1.

get_dtheta_symmetric_Bragg
(E)¶ The angle correction for the symmetric Bragg case:
\[\delta\theta = \chi_0 / \sin{2\theta_B}\]

get_refractive_correction
(E, beamInDotNormal=None, alpha=None)¶ The difference in the glancing angle of incidence for incident and exit waves, Eqs. (2.152) and (2.112) in [Shvydko_XRO]:
\[\theta_c  \theta'_c = \frac{w_H^{(s)}}{2} \left(b  \frac{1}{b} \right) \tan{\theta_c}\]Note
Not valid close to backscattering.
[Shvydko_XRO] Yu. Shvyd’ko, XRay Optics HighEnergyResolution Applications, SpringerVerlag Berlin Heidelberg, 2004.


class
xrt.backends.raycing.materials.
CrystalFcc
(Crystal)¶ A derivative class from
Crystal
that defines the structure factor for an fcc crystal as:\[\begin{split}F_{hkl}^{fcc} = f \times \left\{ \begin{array}{rl} 4 &\mbox{if $h,k,l$ are all even or all odd} \\ 0 &\mbox{ otherwise} \end{array} \right.\end{split}\]

class
xrt.backends.raycing.materials.
CrystalDiamond
(CrystalFcc)¶ A derivative class from
Crystal
that defines the structure factor for a diamondlike crystal as:\[F_{hkl}^{\rm diamond} = F_{hkl}^{fcc}\left(1 + e^{i\frac{\pi}{2} (h + k + l)}\right).\]

class
xrt.backends.raycing.materials.
CrystalSi
(CrystalDiamond)¶ A derivative class from
CrystalDiamond
that defines the crystal dspacing as a function of temperature.
__init__
(*args, **kwargs)¶  tK: float
 Temperature in Kelvin.
 hkl: sequence
 hkl indices.

dl_l
(t=None)¶ Calculates the crystal elongation at temperature t. Uses the parameterization from [Swenson]. Less than 1% error; the reference temperature is 19.9C; data is in units of unitless; t must be in degrees Kelvin.
[Swenson] C.A. Swenson, J. Phys. Chem. Ref. Data 12 (1983) 179

get_Bragg_offset
(E, Eref)¶ Calculates the Bragg angle offset due to a mechanical (mounting) misalignment.
E is the calculated energy of a spectrum feature, typically the edge position.
Eref is the tabulated position of the same feature.

get_a
()¶ Gives the lattice parameter.


class
xrt.backends.raycing.materials.
CrystalFromCell
(Crystal)¶ CrystalFromCell
builds a crystal from cell parameters and atomic positions which can be found e.g. in Crystals.dat of XOP [XOP] or xraylib. Examples:
>>> xtalQu = rm.CrystalFromCell( >>> 'alphaQuartz', (1, 0, 2), a=4.91304, c=5.40463, gamma=120, >>> atoms=[14]*3 + [8]*6, >>> atomsXYZ=[[0.4697, 0., 0.], >>> [0.4697, 0.4697, 1./3], >>> [0., 0.4697, 2./3], >>> [0.4125, 0.2662, 0.1188], >>> [0.1463, 0.4125, 0.4521], >>> [0.2662, 0.1463, 0.2145], >>> [0.1463, 0.2662, 0.1188], >>> [0.4125, 0.1463, 0.2145], >>> [0.2662, 0.4125, 0.5479]]) >>> >>> xtalGr = rm.CrystalFromCell( >>> 'graphite', (0, 0, 2), a=2.456, c=6.696, gamma=120, >>> atoms=[6]*4, atomsXYZ=[[0., 0., 0.], [0., 0., 0.5], >>> [1./3, 2./3, 0.], [2./3, 1./3, 0.5]]) >>> >>> xtalBe = rm.CrystalFromCell( >>> 'Be', (0, 0, 2), a=2.287, c=3.583, gamma=120, >>> atoms=[4]*2, atomsXYZ=[[1./3, 2./3, 0.25], [2./3, 1./3, 0.75]])

__init__
(name='', hkl=[1, 1, 1], a=5.43071, b=None, c=None, alpha=90, beta=90, gamma=90, atoms=[14, 14, 14, 14, 14, 14, 14, 14], atomsXYZ=[[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.25, 0.25, 0.25], [0.25, 0.75, 0.75], [0.75, 0.25, 0.75], [0.75, 0.75, 0.25]], atomsFraction=None, tK=0, t=None, factDW=1.0, geom='Bragg reflected', table='Chantler total', nuPoisson=0.0, calcBorrmann=None, useTT=False, mosaicity=0)¶  name: str
 Crystal name. Not used by xrt.
 hkl: sequence
 hkl indices.
 a, b, c: float
 Cell parameters in Å. a must be given. b, c, if not given, are equlized to a.
 alpha, beta, gamma: float
 Cell angles in degrees. If not given, are equal to 90.
 atoms: list of str or list of int
 List of atoms in the cell given by element Z’s or element names.
 atomsXYZ: list of 3sequences
List of atomic coordinates in cell units.
Note
atoms and atomsXYZ must contain all the atoms, not only the unique ones for a given symmetry group (we do not consider symmetry here). For example, the unit cell of magnetite (Fe3O4) has 3 unique atomic positions and 56 in total; here, all 56 are needed.
 atomsFraction: a list of float or None
 Atomic fractions. If None, all values are 1.
 nuPoisson: float
 Poisson’s ratio. Used to calculate the properties of bent crystals.
 calcBorrmann: str
 Controls the origin of the ray leaving the crystal. Can be ‘None’,
‘uniform’, ‘Bessel’ or ‘TT’. If ‘None’, the point of reflection
is located on the surface of incidence. In all other cases the
coordinate of the exit point is sampled according to the
corresponding distribution: ‘uniform’ is a fast approximation for
thick crystals, ‘Bessel’ is exact solution for the flat crystals,
‘TT’ is exact solution of TakagiTaupin equations for bent and flat
crystals (‘TT’ requires targetOpenCL in the Optical Element to be
not ‘None’ and useTT in the
Crystal
to be ‘True’. Not recommended for crystals thicker than 100 µm due to heavy computational load).  useTT: bool
 Specifies whether the reflectivity will by calculated by analytical formula or by solution of the TakagiTaupin equations (so far only for the Laue geometry). Must be set to ‘True’ in order to calculate the reflectivity of bent crystals.

class
xrt.backends.raycing.materials.
Powder
(CrystalFromCell)¶ A derivative class from
CrystalFromCell
with randomly distributed atomic plane orientations similar to the real polycrystalline powders. The distribution is uniform in the spherical coordinates, so that the angles of longitudinal and transverse deflection (θ and χ) are both functions of uniformly sampled over [0, 1) variables μ and ν: θ = arccos(μ), χ = 2πν.The class parameter hkl defines the highest reflex, so that reflectivities are calculated for all possible combinations of indices [mnp], where 0 ≤ m ≤ h, 0 ≤ n ≤ k, 0 ≤ p ≤ l. Only one reflection with the highest amplitude is picked for each incident ray.
Warning
Heavy computational load. Requires OpenCL.

__init__
(*args, **kwargs)¶  chi: 2list of floats [min, max]
 Limits of the χ angle distribution. Zero and π/2 angles correspond to the positive directions of x and z axes.


class
xrt.backends.raycing.materials.
CrystalHarmonics
(CrystalFromCell)¶ A derivative class from
CrystalFromCell
, used to calculate multiple orders of the given reflex in one run: n*[hkl], where 1 ≤ n ≤ Nmax i.e. [111], [222], [333] or [220], [440], [660]. Only one harmonic with highest reflectivity is picked for each incident ray. Use this class to estimate the efficiency of higher harmonic rejection schemes.Warning
Heavy computational load. Requires OpenCL.

__init__
(*args, **kwargs)¶  Nmax: int
 Specifies the highest order of reflection to be calculated.


class
xrt.backends.raycing.materials.
MonoCrystal
(CrystalFromCell)¶ A derivative class from
CrystalFromCell
, used for calculation of the single crystal diffraction patterns (so far cubic symettries only). Similar to the parent class, parameter hkl defines the cut orientation, whereas Nmax stands for the highest index to consider, i.e. for every ray the code would calculate the range of reflexes from [Nmax, Nmax, Nmax] to [Nmax, Nmax, Nmax] (required amount of reflectivity calculations is therefore 2*(2*Nmax+1)^3 per every ray), but only return one of them regarding their intensities. Brighter reflexes would be selected with higher probability.Warning
Heavy computational load. Requires OpenCL. Decent GPU highly recommended.

__init__
(*args, **kwargs)¶  Nmax: int
 Specifies the highest order of reflection to be calculated.

Apertures¶
Module apertures
defines rectangular and round apertures and a set of
coplanar rectangular apertures. Rectangular apertures may have one or more
defining edges. For example, a simple obstacle, like a beam stop block would
have one edge, a block of frontend slits would have two edges at 90 degrees to
each other, and a collimator would have all four edges.
The classes have useful methods for getting divergence from the aperture size, for setting divergence (calculating the aperture size given the divergence) and for touching the beam with the aperture, i.e. calculating the minimum aperture size that lets the whole beam through.

class
xrt.backends.raycing.apertures.
RectangularAperture
¶ Implements an aperture or an obstacle as a combination of straight edges.

__init__
(bl=None, name='', center=[0, 0, 0], kind=['left', 'right', 'bottom', 'top'], opening=[10, 10, 10, 10], x='auto', z='auto', alarmLevel=None)¶  bl: instance of
BeamLine
 Container for beamline elements. Optical elements are added to its slits list.
 name: str
 Userspecified name, can be used for diagnostics output.
 center: 3sequence of floats
 3D point in global system.
 kind: sequence
 Any combination of ‘top’, ‘bottom’, ‘left’, ‘right’.
 opening: sequence
 Distances (with sign according to the local coordinate system) from the blade edges to the aperture center with the length corresponding to kind.
 x, z: 3tuples or ‘auto’.
Normalized 3D vectors in the global system which determine the local x and z axes lying in the aperture plane. If x is ‘auto’, it is horizontal and normal to the beam line. If z is ‘auto’, it is vertical.
Warning
If you change x and/or z outside of the constructor, you must invoke the method
set_orientation()
. alarmLevel: float or None.
 Allowed fraction of number of rays absorbed at the aperture relative to the number of incident rays. If exceeded, an alarm output is printed in the console.
 bl: instance of

get_divergence
(source)¶ Gets divergences given the blade openings.

prepare_wave
(prevOE, nrays, rw=None)¶ Creates the beam arrays used in wave diffraction calculations. prevOE is the diffracting element: a descendant from
OE
,RectangularAperture
orRoundAperture
. nrays: if int, specifies the number of randomly distributed samples over the slit area; if 2tuple of ints, specifies (nx, ny) sizes of a uniform mesh of samples.

propagate
(beam=None, needNewGlobal=False)¶ Assigns the “lost” value to beam.state array for the rays intercepted by the aperture. The “lost” value is
self.ordinalNum  1000.

set_divergence
(source, divergence)¶ Gets the blade openings given divergences. divergence is a sequence corresponding to kind

touch_beam
(beam)¶ Adjusts the aperture (i.e. sets self.opening) so that it touches the beam.


class
xrt.backends.raycing.apertures.
RoundAperture
¶ Implements a round aperture meant to represent a pipe or a flange.

__init__
(bl=None, name='', center=[0, 0, 0], r=1, x='auto', z='auto', alarmLevel=None)¶ A round aperture aperture.
r is the radius.
 x, z: 3tuples or ‘auto’.
Normalized 3D vectors in the global system which determine the local x and z axes lying in the aperture plane. If x is ‘auto’, it is horizontal and normal to the beam line. If z is ‘auto’, it is vertical.
Warning
If you change x and/or z outside of the constructor, you must invoke the method
set_orientation()
.

get_divergence
(source)¶ Gets the full divergence given the aperture radius.

propagate
(beam=None, needNewGlobal=False)¶ Assigns the “lost” value to beam.state array for the rays intercepted by the aperture. The “lost” value is
self.ordinalNum  1000.


class
xrt.backends.raycing.apertures.
RoundBeamStop
¶ Implements a round beamstop. Descends from RoundAperture and has the same parameters.

class
xrt.backends.raycing.apertures.
DoubleSlit
¶ Implements an aperture or an obstacle with a combination of horizontal and/or vertical edge(s).

__init__
(*args, **kwargs)¶ Same parameters as in
RectangularAperture
and additionally shadeFraction as a value from 0 to 1.


class
xrt.backends.raycing.apertures.
PolygonalAperture
¶ Implements an aperture or an obstacle defined as a set of polygon vertices.

__init__
(bl=None, name='', center=[0, 0, 0], opening=None, x='auto', z='auto', alarmLevel=None)¶  bl: instance of
BeamLine
 Container for beamline elements. Optical elements are added to its slits list.
 name: str
 Userspecified name, can be used for diagnostics output.
 center: 3sequence of floats
 3D point in global system.
 opening: sequence
 Coordinates [(x0, y0),…(xN, yN)] of the polygon vertices.
 x, z: 3tuples or ‘auto’.
Normalized 3D vectors in the global system which determine the local x and z axes lying in the aperture plane. If x is ‘auto’, it is horizontal and normal to the beam line. If z is ‘auto’, it is vertical.
Warning
If you change x and/or z outside of the constructor, you must invoke the method
set_orientation()
. alarmLevel: float or None.
 Allowed fraction of number of rays absorbed at the aperture relative to the number of incident rays. If exceeded, an alarm output is printed in the console.
 bl: instance of

Screens¶
Module screens
defines a flat screen and a
hemispheric screen that intercept a beam and give its image.

class
xrt.backends.raycing.screens.
Screen
¶ Flat screen for beam visualization.

__init__
(bl=None, name='', center=[0, 0, 0], x='auto', z='auto', compressX=None, compressZ=None)¶  bl: instance of
BeamLine
 Container for beamline elements. Optical elements are added to its screens list.
 name: str
 Userspecified name, can be used for diagnostics output.
 center: tuple of 3 floats
 3D point in the global system.
 x, z: 3sequence or ‘auto’.
Normalized 3D vectors in the global system which determine the local x and z axes lying in the screen plane. If x is ‘auto’, it is horizontal and perpendicular to the beam line. If z is ‘auto’, it is vertical.
Warning
If you change x and/or z outside of the constructor, you must invoke the method
set_orientation()
. compressX, compressZ: float
 Multiplicative compression coefficients for the corresponding axes. Typically are not needed. Can be useful to account for the viewing camera magnification or when the camera sees the screen at an angle.
 bl: instance of

expose
(beam=None, onlyPositivePath=False)¶ Exposes the screen to the beam. beam is in global system, the returned beam is in local system of the screen and represents the desired image.

prepare_wave
(prevOE, dim1, dim2, dy=0, rw=None, condition=None)¶ Creates the beam arrays used in wave diffraction calculations. prevOE is the diffracting element: a descendant from
OE
,RectangularAperture
orRoundAperture
. dim1 and dim2 are x and z arrays for a flat screen or phi and theta arrays for a hemispheric screen. The two arrays are generally of different 1D shapes. They are used to create a 2D mesh bymeshgrid
. condition: a callable defined in the user script with two flattened
meshgrid arrays as inputs and outputs. Can be used to select wave samples. An example:
def condition(d1s, d2s): cond = d1s**2 + d2s**2 <= pinholeDia**2 / 4 # in a pinhole return d1s[cond], d2s[cond]

set_orientation
(x=None, z=None)¶ Determines the local x, y and z in the global system.


class
xrt.backends.raycing.screens.
HemisphericScreen
¶ Hemispheric screen for beam visualization.

__init__
(bl=None, name='', center=[0, 0, 0], R=1000.0, x='auto', z='auto', phiOffset=0, thetaOffset=0)¶  x, z: 3tuples or ‘auto’. Normalized 3D vectors in the global system
which determine the local x and z axes of the hemispheric screen. If x (the origin of azimuthal angle φ) is ‘auto’, it coincides with the beamline’s y; if z (the polar axis) is ‘auto’, it is coincides with the beamline’s x. The equator plane is then vertical. The polar angle θ is counted from π/2 to π/2 with 0 at the equator and π/2 at the polar axis direction.
Warning
If you change x and/or z outside of the constructor, you must invoke the method
set_orientation()
. R: float
 Radius of the hemisphere in mm.

Wave propagation (diffraction)¶
Time dependent diffraction¶
We start from the Kirchhoff integral theorem in the general (timedependent) form [Born & Wolf]:
\[V(r,t)=\frac 1{4\pi }\int _S\left\{[V]\frac{\partial }{\partial n} \left(\frac 1 s\right)\frac 1{\mathit{cs}}\frac{\partial s}{\partial n}\left[\frac{\partial V}{\partial t}\right]\frac 1 s\left[\frac {\partial V}{\partial n}\right]\right\}\mathit{dS},\]
where the integration is performed over the selected surface \(S\), \(s\) is the distance between the point \(r\) and the running point on surface \(S\), \(\frac{\partial }{\partial n}\) denotes differentiation along the normal on the surface and the square brackets on \(V\) terms denote retarded values, i.e. the values at time \(t − s/c\). \(V\) is a scalar wave here but can represent any component of the actual electromagnetic wave provided that the observation point is much further than the wave length (surface currents are neglected here). \(V\) depends on position and time; this is something what we do not have in ray tracing. We obtain it from ray characteristics by:
\[V(s,t)=\frac 1{\sqrt{2\pi }}\int U_{\omega }(s)e^{i\omega t}d\omega,\]
where \(U_{\omega }(s)\) is interpreted as a monochromatic wave field and therefore can be associated with a ray. Here, this is any component of the ray polarization vector times its propagation factor \(e^{ikl(s)}\). Substituting it into the Kirchhoff integral yields \(V(r,t)\). As we visualize the wave fields in space and energy, the obtained \(V(r,t)\) must be backFourier transformed to the frequency domain represented by a new resampled energy grid:
\[U_{\omega '}(r)=\frac 1{\sqrt{2\pi }}\int V(r,t)e^{i\omega 't} \mathit{dt}.\]
Ingredients:
\[ \begin{align}\begin{aligned}[V]\frac{\partial }{\partial n}\left(\frac 1 s\right)=\frac{(\hat {\vec s}\cdot \vec n)}{s^2}[V],\\\frac 1{\mathit{cs}}\frac{\partial s}{\partial n}\left[\frac{\partial V}{\partial t}\right]=\frac{ik} s(\hat{\vec s}\cdot \vec n)[V],\\\frac 1 s\left[\frac{\partial V}{\partial n}\right]=\frac{ik} s(\hat {\vec l}\cdot \vec n)[V],\end{aligned}\end{align} \]
where the hatted vectors are unit vectors: \(\hat {\vec l}\) for the incoming direction, \(\hat {\vec s}\) for the outgoing direction, both being variable over the diffracting surface. As \(1/s\ll k\), the 1st term is negligible as compared to the second one.
Finally,
\[U_{\omega '}(r)=\frac{i}{8\pi ^2\hbar ^2c}\int e^{i(\omega '\omega)t} \mathit{dt}\int \frac E s\left((\hat{\vec s}\cdot \vec n)+(\hat{\vec l} \cdot \vec n)\right)U_{\omega }(s)e^{ik(l(s)+s)}\mathit{dS}\mathit{dE}.\]
The timedependent diffraction integral is not yet implemented in xrt.
Stationary diffraction¶
If the time interval \(t\) is infinite, the forward and back Fourier transforms give unity. The Kirchhoff integral theorem is reduced then to its monochromatic form. In this case the energy of the reconstructed wave is the same as that of the incoming one. We can still use the general equation, where we substitute:
\[\delta (\omega \omega ')=\frac 1{2\pi }\int e^{i(\omega '\omega )t} \mathit{dt},\]
which yields:
\[U_{\omega }(r)=\frac {i k}{4\pi }\int \frac1 s\left((\hat{\vec s}\cdot \vec n)+(\hat{\vec l}\cdot \vec n)\right)U_{\omega }(s)e^{ik(l(s)+s)} \mathit{dS}.\]
How we treat nonmonochromaticity? We repeat the sequence of raytracing from the source down to the diffracting surface for each energy individually. For synchrotron sources, we also assume a single electron trajectory (so called “filament beam”). This single energy contributes fully coherently into the diffraction integral. Different energies contribute incoherently, i.e. we add their intensities, not amplitudes.
The input field amplitudes can, in principle, be taken from raytracing, as it
was done by [Shi_Reininger] as \(U_\omega(s) = \sqrt{I_{ray}(s)}\). This
has, however, a fundamental difficulty. The notion “intensity” in many ray
tracing programs, as in Shadow used in [Shi_Reininger], is different from the
physical meaning of intensity: “intensity” in Shadow is a placeholder for
reflectivity and transmittivity. The real intensity is represented by the
density of rays – this is the way the rays were sampled, while each ray has
\(I_{ray}(x, z) = 1\) at the source [shadowGuide], regardless of the
intensity profile. Therefore the actual intensity must be reconstructed.
We tried to overcome this difficulty by computing the density of rays by (a)
histogramming and (b) kernel density estimation [KDE]. However, the easiest
approach is to sample the source with uniform ray density (and not proportional
to intensity) and to assign to each ray its physical wave amplitudes as s and
p projections. In this case we do not have to reconstruct the physical
intensity. The uniform ray density is an option for the geometric and
synchrotron sources in sources
.
Notice that this formulation does not require paraxial propagation and thus xrt is more general than other wave propagation codes. For instance, it can work with gratings and FZPs where the deflection angles may become large.
[Shi_Reininger]  (1, 2) X. Shi, R. Reininger, M. Sanchez del Rio & L. Assoufid, A hybrid method for Xray optics simulation: combining geometric raytracing and wavefront propagation, J. Synchrotron Rad. 21 (2014) 669–678. 
[shadowGuide] 

[KDE]  Michael G. Lerner (mglerner) (2013) http://www.mglerner.com/blog/?p=28 
Normalization¶
The amplitude factors in the Kirchhoff integral assure that the diffracted wave has correct intensity and flux. This fact appears to be very handy in calculating the efficiency of a grating or an FZP in a particular diffraction order. Without proper amplitude factors one would need to calculate all the significant orders and renormalize their total flux to the incoming one.
The resulting amplitude is correct provided that the amplitudes on the diffracting surface are properly normalized. The latter are normalized as follows. First, the normalization constant \(X\) is found from the flux integral:
\[F = X^2 \int \left(E_s^2 + E_p^2 \right) (\hat{\vec l}\cdot \vec n) \mathit{dS}\]
by means of its MonteCarlo representation:
\[X^2 = \frac{F N }{\sum \left(E_s^2 + E_p^2 \right) (\hat{\vec l}\cdot \vec n) S} \equiv \frac{F N }{\Sigma(J\angle)S}.\]
The area \(S\) can be calculated by the user or it can be calculated automatically by constructing a convex hull over the impact points of the incoming rays. The voids, as in the case of a grating (shadowed areas) or an FZP (the opaque zones) cannot be calculated by the convex hull and such cases must be carefully considered by the user.
With the above normalization factors, the Kirchhoff integral calculated by MonteCarlo sampling gives the polarization components \(E_s\) and \(E_p\) as (\(\gamma = s, p\)):
\[E_\gamma(r) = \frac{\sum K(r, s) E_\gamma(s) X S}{N} = \sum{K(r, s) E_\gamma(s)} \sqrt{\frac{F S} {\Sigma(J\angle) N}}.\]
Finally, the Kirchhoff intensity \(\left(E_s^2 + E_p^2\right)(r)\) must be integrated over the screen area to give the flux.
Note
The above normalization steps are automatically done inside
diffract()
.
Sequential propagation¶
In order to continue the propagation of a diffracted field to the next optical element, not only the field distribution but also local propagation directions are needed, regardless how the radiation is supposed to be propagated further downstream: as rays or as a wave. The directions are given by the gradient applied to the field amplitude. Because \(1/s\ll k\) (validity condition for the Kirchhoff integral), the by far most significant contribution to the gradient is from the exponent function, while the gradient of the preexponent factor is neglected. The new wave directions are thus given by a Kirchhofflike integral:
\[{\vec \nabla } U_{\omega }(r) = \frac {k^2}{4\pi } \int \frac {\hat{\vec s}} s \left((\hat{\vec s}\cdot \vec n)+(\hat{\vec l}\cdot \vec n)\right) U_{\omega }(s)e^{ik(l(s)+s)} \mathit{dS}.\]
The resulted vector is complex valued. Taking the real part is not always correct as the vector components may happen to be (almost) purely imaginary. Our solution is to multiply the vector components by a conjugate phase factor of the largest vector component and only then to take the real part.
The correctness of this approach to calculating the local wave directions can be verified as follows. We first calculate the field distribution on a flat screen, together with the local directions. The wave front, as the surface normal to these directions, is calculated by linear integrals of the corresponding angular projections. The calculated wave front surface is then used as a new screen, where the diffracted wave is supposed to have a constant phase. This is indeed demonstrated in the example applying phase as color axis. The sharp phase distribution is indicative of a true wave front, which in turn justifies the correct local propagation directions.
After having found two polarization components and new directions on the receiving surface, the last step is to reflect or refract the wave samples as rays locally, taking the complex refractive indices for both polarizations. This approach is applicable also to wave propagation regime, as the reflectivity values are purely local to the impact points.
Usage¶
Warning
You need a good graphics card for running these calculations!
Note
OpenCL platforms/devices can be inspected in xrtQook (‘GPU’ button)
Warning
Long calculation on GPU in Windows may result in the system message “Display driver stopped responding and has recovered” or python RuntimeError: out of resources (yes, it’s about the driver response, not the lack of memory). The solution is to change TdrDelay registry key from the default value of 2 seconds to some hundreds or even thousands. Please refer to https://msdn.microsoft.com/enus/library/windows/hardware/ff569918(v=vs.85).aspx
Tip
If you plan to utilize xrt on a multiGPU platform under Linux, we recommend KDE based systems (e.g. Kubuntu). It is enough to install the package
fglrxupdatesdev
without the complete graphics driver.
In contrast to ray propagation, where passing through an optical element requires only one method (“reflect” for a reflective/refractive surface or “propagate” for a slit), wave propagation in our implementation requires two or three methods:
prepare_wave
creates points on the receiving surface where the diffraction integral will be calculated. For an optical element or a slit the points are uniformly randomly distributed (therefore reasonable physical limits must be given); for a screen the points are determined by the screen pixels. The returned wave is an instance of classBeam
, where the arrays x, y, z are local to the receiving surface.diffract
(as imported function from xrt.backends.raycing.waves or as a method of the diffracted beam) takes a beam local to the diffracting surface and calculates the diffraction integrals at the points prepared byprepare_wave
. There are five scalar integrals: two for Es and Ep and three for the components of the diffracted direction (see the previous Section). All five integrals are calculated in the coordinates local to the diffracting surface but at the method’s end these are transformed to the local coordinates of the receiving surface (remained in the wave container) and to the global coordinate system (a Beam object returned bydiffract
).
For slits and screens, the two steps above are sufficient. For an optical element, another step is necessary.
reflect
method of the optical element is meant to take into account its material properties. The intersection points are already known, as provided by the previousprepare_wave
, which can be reported toreflect
bynoIntersectionSearch=True
. Here,reflect
takes the beam right before the surface and propagates it to right after it. As a result of such a zerolength travel, the wave gets no additional propagation phase but only a complexvalued reflectivity coefficient and a new propagation direction.
These three methods are enough to describe wave propagation through the
complete beamline. The first two methods, prepare_wave
and diffract
,
are split from each other because the diffraction calculations may need several
repeats in order to accumulate enough wave samples for attaining dark field at
the image periphery. The second method can reside in a loop that will
accumulate the complex valued field amplitudes in the same beam arrays defined
by the first method. In the supplied examples, prepare_wave
for the last
screen is done before a loop and all the intermediate diffraction steps are
done within that loop.
The quality of the resulting diffraction images is mainly characterized by the
blackness of the dark field – the area of expected zero intensity. If the
statistics is not sufficient, the dark area is not black, and can even be
bright enough to mask the main spot. The contrast depends both on beamline
geometry (distances) and on the number of wave field samples (a parameter for
prepare_wave
). Shorter distances require more samples for the same quality,
and for the distances shorter than a few meters one may have to reduce the
problem dimensionality by cutting in horizontal or vertical, see the
examples of SoftiMAX. In the console output, diffract
reports on samples per zone (meaning per Fresnel zone). As a rule of thumb,
this figure should be greater than ~10^{4} for a good resulting quality.

OE.
prepare_wave
(prevOE, nrays, shape='auto', area='auto', rw=None) Creates the beam arrays used in wave diffraction calculations. prevOE is the diffracting element: a descendant from
OE
,RectangularAperture
orRoundAperture
. nrays: if int, specifies the number of randomly distributed samples the surface withinself.limPhysX
limits; if 2tuple of ints, specifies (nx, ny) sizes of a uniform mesh of samples.

RectangularAperture.
prepare_wave
(prevOE, nrays, rw=None) Creates the beam arrays used in wave diffraction calculations. prevOE is the diffracting element: a descendant from
OE
,RectangularAperture
orRoundAperture
. nrays: if int, specifies the number of randomly distributed samples over the slit area; if 2tuple of ints, specifies (nx, ny) sizes of a uniform mesh of samples.

Screen.
prepare_wave
(prevOE, dim1, dim2, dy=0, rw=None, condition=None) Creates the beam arrays used in wave diffraction calculations. prevOE is the diffracting element: a descendant from
OE
,RectangularAperture
orRoundAperture
. dim1 and dim2 are x and z arrays for a flat screen or phi and theta arrays for a hemispheric screen. The two arrays are generally of different 1D shapes. They are used to create a 2D mesh bymeshgrid
. condition: a callable defined in the user script with two flattened
meshgrid arrays as inputs and outputs. Can be used to select wave samples. An example:
def condition(d1s, d2s): cond = d1s**2 + d2s**2 <= pinholeDia**2 / 4 # in a pinhole return d1s[cond], d2s[cond]

xrt.backends.raycing.waves.
diffract
(oeLocal, wave, targetOpenCL='auto', precisionOpenCL='auto')¶ Calculates the diffracted field – the amplitudes and the local directions – contained in the wave object. The field on the diffracting surface is given by oeLocal. You can explicitly change OpenCL settings targetOpenCL and precisionOpenCL which are initially set to ‘auto’, see their explanation in
xrt.backends.raycing.sources.Undulator
.
Coherence signatures¶
A standard way to define coherence properties is via mutual intensity J and complex degree of coherence j (DoC, normalized J):
\[\begin{split}J(x_1, y_1, x_2, y_2) \equiv J_{12} = \left<E(x_1, y_1)E^{*}(x_2, y_2)\right>\\ j(x_1, y_1, x_2, y_2) \equiv j_{12} = \frac{J_{12}}{\left(J_{11}J_{22}\right)^{1/2}},\end{split}\]
where the averaging \(\left<\ \right>\) in \(J_{12}\) is over different
realizations of filament electron beam (one realization per repeat
) and is
done for the field components \(E_s\) or \(E_p\) of a field diffracted
onto a given \((x, y)\) plane.
Both functions are Hermitian in respect to the exchange of points 1 and 2. There are two common ways of working with them:
The horizontal and vertical directions are considered independently by placing the points 1 and 2 on a horizontal or vertical line symmetrically about the optical axis. DoC thus becomes a 1D function dependent on the distance between the points, e.g. as \(j_{12}^{\rm hor}=j(x_1x_2)\). The intensity distribution is also determined over the same line as a 1D positional function, e.g. as \(I(x)\).
The widths \(\sigma_x\) and \(\xi_x\) of the distributions \(I(x)\) and \(j(x_1x_2)\) give the coherent fraction \(\zeta_x\) [Vartanyants2010]
\[\zeta_x = \left(4\sigma_x^2/\xi_x^2 + 1\right)^{1/2}.\]The transverse field distribution can be analized integrally (not split into the horizontal and vertical projections) by performing the modal analysis consisting of solving the eigenvalue problem for the matrix \(J^{tr=1}_{12}\) – the matrix \(J_{12}\) normalized to its trace – and doing the standard eigendecomposition:
\[J^{tr=1}(x_1, y_1, x_2, y_2) = \sum_i{w_i V_i(x_1, y_1)V_i^{+}(x_2, y_2)},\]with \(w_i, V_i\) being the ith eigenvalue and eigenvector. \(w_0\) is the fraction of the total flux contained in the 0th (coherent) mode or coherent flux fraction.
Note
The matrix \(J_{12}\) is of the size (N_{x}×N_{y})², i.e. squared total pixel size of the image! In the current implementation, we use
eigh()
method fromscipy.linalg
, where a feasible image size should not exceed ~100×100 pixels (i.e. ~10^{8} size of \(J_{12}\)).Note
For a fully coherent field \(j_{12}\equiv1\) and \(w_0=1, w_i=0\ \forall i>0\), \(V_0\) being the coherent field.
We also propose a third method that results in the same figures as the second method above.
It uses Principal Component Analysis (PCA) applied to the filament images \(E(x, y)\). It consists of the following steps.
 Out of r repeats of \(E(x, y)\) build a data matrix \(D\) with N_{x}×N_{y} rows and r columns.
 The matrix \(J_{12}\) is equal to the product \(DD^{+}\). Instead of solving this huge eigenvalue problem of (N_{x}×N_{y})² size, we solve a typically smaller matrix \(D^{+}D\) of the size r².
 The biggest r eigenvalues of \(J_{12}\) are equal to those of \(D^{+}D\) [proof to present in the coming paper]. To find the primary (biggest) eigenvalue is the main objective of the modal analysis (item 2 above); PCA can provide it much easier due to the smaller size of the problem.
 Also the eigen modes of \(J_{12}=DD^{+}\) can be found by PCA via the eigenvectors \(v\)’s of \(D^{+}D\). The matrix \(Dv_iv_i^{+}\) is of the size of \(D\) and has all the columns proportional to each other [proof to present in the coming paper]. These columns are the ith principal components for the corresponding columns of \(D\). Being normalized, all the columns become equal and give the ith eigenvector of \(J_{12}\).
Finally, PCA gives exactly the same information as the direct modal analysis (method No 2 above) but is cheaper to calculate by many orders of magnitude.
One can define another measure of coherence as a single number, termed as degree of transverse coherence (DoTC) [Saldin2008]:
[Saldin2008]  E.L. Saldin, E.A. Schneidmiller, M.V. Yurkov, Coherence properties of the radiation from Xray free electron laser, Opt. Commun. 281 (2008) 1179–88. 
[Vartanyants2010]  I.A. Vartanyants and A. Singer, Coherence properties of hard xray synchrotron sources and xray freeelectron lasers, New Journal of Physics 12 (2010) 035004. 
We propose to calculate DoTC from the matrix traces [derivation to present in the coming paper] as:
 DoTC = Tr(J²)/Tr²(J).
 DoTC = Tr(D^{+}DD^{+}D)/Tr²(D^{+}D), with the matrix D defined above. The exactly same result as in (a) but obtained with smaller matrices.
Note
A good test for the correctness of the obtained coherent fraction is to find it at various positions on propagating in free space, where the result is expected to be invariant. As appears in the examples of SoftiMAX, the analysis based on DoC never gives an invariant coherent fraction at the scanned positions around the focus. The primary reason for this is the difficulty in the determination of the width of DoC, for the latter typically being a complexshaped oscillatory curve. In contrast, the modal analysis (the PCA implementation is recommended) and the DoTC give the expected invariance.
Coherence analysis and related plotting¶
The module coherence
has functions for 1D and 2D
analysis of coherence and functions for 1D plotting of degree of coherence and
and 2D plotting of eigen modes.
The input for the analysis functions is a 3D stack of field images. It can be
obtained directly from the undulator class, or from a plot object after several
repeats of wave propagation of a filament beam through a beamline. Examples can
be found in ...\tests\raycing\test_coherent_fraction_stack.py
and in
SoftiMAX at MAX IV.

xrt.backends.raycing.coherence.
calc_1D_coherent_fraction
(U, axisName, axis, p=0)¶ Calculates 1D degree of coherence (DoC). From its width in respect to the width of intensity distribution also infers the coherent fraction. Both widths are rms. The one of intensity is calculated over the whole axis, the one of DoC is calculated between the first local minima around the center provided that these minima are lower than 0.5.
 U: complex valued ndarray, shape(repeats, nx, ny)
 3D stack of field images. For a 1D cut along axis, the middle of the other dimension is sliced out.
 axis: str, one of ‘x’ or (‘y’ or ‘z’)
 Specifies the 1D axis of interest.
 p: float, distance to screen
 If nonzero, the calculated mutual intensity will be divided by p². This is useful to get proper physical units of the returned intensity if the function is applied directly to the field stacks given by Undulator.multi_electron_stack() that is calculated in angular units.
Returns a tuple of mutual intensity, 1D intensity, 1D DoC, rms width of intensity, rms width of DoC (between the local minima, see above), the position of the minima (only the positive side) and the coherent fraction. This tuple can be fed to the plotting function
plot_1D_degree_of_coherence()
.

xrt.backends.raycing.coherence.
plot_1D_degree_of_coherence
(data1D, axisName, axis, unit='mm', fig2=None, isIntensityNormalized=False, locLegend='auto')¶ Provides two plots: a 2D plot of mutual intensity and a 1D plot of intensity and DoC. The latter plot can be shared between the two 1D axes if this function is invoked two times.
data1D: tuple returned by
calc_1D_coherent_fraction()
.axisName: str, used in labels.
axis: 1D array of abscissa.
unit: str, used in labels.
fig2: matplotlib figure object, if needed for shared plotting of two 1D axes.
isIntensityNormalized: bool, controls the intensity axis label.
locLegend: str, legend location in matplotlib style.
Returns the two figure objects for the user to add suptitles and to export to image files.
Plot examples for one and two 1D curves:

xrt.backends.raycing.coherence.
calc_degree_of_transverse_coherence_4D
(J)¶ Calculates DoTC from the mutual intensity J as Tr(J²)/Tr²(J). This function should only be used for demonstration purpose. There is a faster alternative:
calc_degree_of_transverse_coherence_PCA()
.

xrt.backends.raycing.coherence.
calc_degree_of_transverse_coherence_PCA
(U)¶ Calculates DoTC from the field stack U. The field images of U are flattened to form the matrix D shaped as (repeats, nx×ny). DoTC = Tr(D^{+}DD^{+}D)/Tr²(D^{+}D), which is equal to the original DoTC for the mutual intensity J: DoTC = Tr(J²)/Tr²(J).

xrt.backends.raycing.coherence.
calc_eigen_modes_4D
(J, eigenN=4)¶ Solves the eigenvalue problem for the mutual intensity J. This function should only be used for demonstration purpose. There is a much faster alternative:
calc_eigen_modes_PCA()
.

xrt.backends.raycing.coherence.
calc_eigen_modes_PCA
(U, eigenN=4, maxRepeats=None, normalize=False)¶ Solves the PCA problem for the field stack U shaped as (repeats, nx, ny). The field images are flattened to form the matrix D shaped as (repeats, nx×ny). The eigenvalue problem is solved for the matrix D^{+}*D*.
Returns a tuple of two arrays: eigenvalues in a 1D array and eigenvectors as columns of a 2D array. This is a much faster and exact replacement of the full eigen mode decomposition by
calc_eigen_modes_4D()
.eigenN sets the number of returned eigen modes. If None, the number of modes is inferred from the shape of the field stack U and is equal to the number of macroelectrons (repeats). If given as a 2tuple, eigenN refers to the number of eigenvalues and eigenvectors separately.
if maxRepeats are given, the stack is sliced up to that number. This option is introduced in order not to make the eigenvalue problem too big.
If normalize is True, the eigenvectors are normalized, otherwise they are the PCs of the field stack in the original field units.

xrt.backends.raycing.coherence.
plot_eigen_modes
(x, y, w, v, xlabel='', ylabel='')¶ Provides 2D plots of the first 4 eigen modes in the given x and y coordinates. The eigen modes are specified by eigenvalues and eigenvectors w and v returned by
calc_eigen_modes_PCA()
orcalc_eigen_modes_4D()
.Plot example:
Typical logic for a wave propagation study¶
 According to [Wolf], the visibility is not affected by a small bandwidth. It is therefore recommended to first work with strictly monochromatic radiation and check nonmonochromaticity at the end.
 Do a hybrid propagation (waves start somewhere closer to the end) at zero emittance. Put the screen at various positions around the main focus for seeing the depth of focus.
 Examine the focus and the footprints and decide if vertical and horizontal cuts are necessary (when the dark field does not become black enough at a reasonably high number of wave samples).
 Run nonzero electron emittance, which spoils both the focal size and the degree of coherence.
 Try various ways to improve the focus and the degree of coherence: e.g. decrease the exit slit or elongate interoptics distances.
 Do hybrid propagation for a finite energy band to study the chromaticity in focus position.
[Wolf]  E. Wolf. Introduction to the theory of coherence and polarization of light. Cambridge University Press, Cambridge, 2007. 
Coherent mode decomposition and propagation¶
The most straightforward way for studying the coherent portion of synchrotron light is first to propagate a large collection of filament beams (or macroelectrons) onto a surface of interest (typically, the focal plane at the sample position) and then perform the modal analysis. The propagation itself may become expensive when optical elements are situated close to each other or when the analysis plane is out of focus. These cases require a large number of wave samples which would allow only a small number of filament beams processed in a reasonable time. In such cases, an inverse approach should be preferred: first do the modal analysis of the source radiation and then propagate the 0th or a few main modes corresponding to the biggest eigenvalues (flux fractions).
In xrt, we calculate a collection of undulator field realizations at the first beamline slit (a frontend slit), transform these fields into modes and save them for the further propagation along the beamline. A selected number of modes, also optionally a selected number of original fields, are prepared for the following three ways of propagation:
 As waves. The wave samples are sequentially diffracted from the first slit by calculating the Kirchhoff diffraction integral, see Sequential propagation.
 As hybrid waves. The wave samples at the first slit are treated as rays with the directions projected from the source center (for the modes) or from the filament beam position (for the fields). The beams of these rays can be propagated in ray tracing and then diffracted closer to the end of the beamline. Note that if the beams were propagated as rays down to the image plane, the individual fields would be seen as filament beam dots and the individual modes would be seen as centered dots, i.e. the source would have zero size.
 As rays. The wave samples at the first slit are propagated backward to the source plane, which gives the field distribution in real space. The field at the first slit is sampled by its intensity; these samples give the ray directions. These beams are suitable for ray traycing down to the image plane. They are not suitable for wave propagation, as the propagation phase of each ray is destroyed by the above resampling procedure.

xrt.backends.raycing.modes.
make_and_save_modes
(bl, nsamples, nElectrons, nElectronsSave, nModes, fixedEnergy, phaseEsEp=0, output='all', basename='local', limitsOrigin=[])¶ Produces pickled files of nModes wave modes and nElectronsSave wave fields. The beamline object bl must have at least one aperture; the first of them fill be used to generate nElectrons fields for the eigenmode decomposition. The aperture will be sampled by nsamples wave samples. The fields are normalized such that the intensity sum over the sumples and over nElectrons gives the total flux returned as totalFlux in use_saved(). Note that the total flux is made independent (in average) of nElectrons.
fixedEnergy is photon energy.
phaseEsEp is phase difference between Es and Ep components.
output can be ‘all’ or any combination of words ‘wave’, ‘hybr’, ‘rays’ in one string. The case ‘rays’ can take quite long time for field resampling.
basename is the output file name that will be prepended by the selected ‘wave’, ‘hybr’, ‘rays’ output modes and appended by .pickle.
limitsOrigin as [xmin, xmax, zmin, zmax] must be given for generating ‘rays’.

xrt.backends.raycing.modes.
use_saved
(what, basename)¶ Loads the saved modes and fields produced by make_and_save_modes().
what is a string starting with one of ‘wave’, ‘hybr’, ‘rays’ and ending with one of ‘modes’, ‘fields’.
basename is the same parameter used in make_and_save_modes().
Returnes a tuple (savedBeams, wAll, totalFlux), where savedBeams is a list of beam objects corresponding to what. wAll is a list of all eigenvalues (flux weights) of the eigenmode decomposition. totalFlux is described in make_and_save_modes().
See a usage example in
/examples/withRaycing/11_Waves/coherentModePropagation.py
.
A model case of 5 macroelectrons¶
This example was made for the following beamline scheme, where the mirror M1 is an ideal ellipsoid providing 1:1 image of the source.
The eigenmode decomposition at the first slit is accurate within the machine precision, i.e. the total flux summed over all individual fields is equal to the total flux summed over all modes, here down to ~1e16 relative error. This fact can be proven by comparing the last animation frames of the ‘field’ images vs. ‘mode’ images. This comparison “all fields” vs. “all modes” was the main objective of this “a few macroelectrons” test; in usual studies one would typically need a few thousand macroelectrons and just a small number of modes.
Individual fields
focus/source  FE slit  

rays  
hybrid  
wave 
Coherent modes
focus/source  FE slit  

rays  
hybrid  
wave 
A model case of 5000 macroelectrons¶
Coherent radiation
The coherent radiation in the 0th mode has a weight (coherent fraction) ≈3.66%. Here it was individually propagated in three different ways:
focus/source  FE slit  

rays  
hybrid  
wave 