xrt (XRayTracer)¶

DOI:
GitHub:
PyPI:: version history, plans
License:
Authors:: Konstantin Klementiev (MAX IV Laboratory), Roman Chernikov (Canadian Light Source)

Package xrt is a python software library for ray tracing and wave propagation in x-ray regime. It is primarily meant for modeling synchrotron sources, beamlines and beamline elements. Includes a GUI for creating a beamline and interactively viewing it in 3D.

Intensity of undulator radiation on a transverse flat screen. The long axis is energy. One can observe two harmonics crossing the central line: the one at 2400 eV is odd and the one at 4800 eV is even. Intensity of undulator radiation on a transverse flat screen. The long axis is energy. One can observe two harmonics crossing the central line: the one at 2400 eV is odd and the one at 4800 eV is even.

Absorbed power and power density on a mirror at varying pitch angle. Absorbed power and power density on a mirror at varying pitch angle.

Vortex beam with l=1 and p=1. Colored by phase. Vortex beam with l=1 and p=1. Colored by phase.

GUIs¶

xrtQook – a GUI for creating a beamline¶

The main interface to xrt is through a python script. Many examples of such scripts can be found in the supplied folders examples and tests. The script imports the modules of xrt, instantiates beamline parts, such as synchrotron or geometric sources, various optical elements, apertures and screens, specifies required materials for reflection, refraction or diffraction, defines plots and sets job parameters.

The Qt tool xrtQook takes these ingredients as GUI elements and prepares a ready to use script that can be run within the tool itself or in an external Python context. xrtQook has a parallelly updated help panel that provides a complete list of parameters for the used objects. xrtQook writes/reads the recipes of beamlines into/from xml files.

In the present version, xrtQook does not provide automated generation of scans and does not create wave propagation sequences. For these two tasks, the corresponding script parts have to be written manually based on the supplied examples and the present documentation.

See a brief tutorial for xrtQook.

A view of xrtQook with an empty beamline tree on the left and a help panel on the right.

Features of xrt¶

Rays and waves. Classical ray tracing and wave propagation via Kirchhoff integrals, also freely intermixed. No further approximations, such as thin lens or paraxial. The optical surfaces may have figure errors, analytical or measured. In wave propagation, partially coherent radiation is treated by incoherent addition of coherently diffracted fields generated per electron. Propagation of individual coherent source modes is possible as waves, hybrid waves (i.e. partially as rays and then as waves) and only rays.
Publication quality graphics. 1D and 2D position histograms are simultaneously coded by hue and brightness. Typically, colors represent energy and brightness represents beam intensity. The user may select other quantities to be encoded by colors: angular and positional distributions, various polarization properties, beam categories, number of reflections, incidence angle etc. Brightness can also encode partial flux for a selected polarization and incident or absorbed power. Publication quality plots are provided by matplotlib with image formats PNG, PostScript, PDF and SVG.
Unlimited number of rays. The colored histograms are cumulative. The accumulation can be stopped and resumed.
Parallel execution. xrt can be run in parallel in several threads or processes (can be opted), which accelerates the execution on multi-core computers. Alternatively, xrt can use the power of GPUs via OpenCL for running special tasks such as the calculation of an undulator source or performing wave propagation.
Scripting in Python. xrt can be run within Python scripts to generate a series of images under changing geometrical or physical parameters. The image brightness and 1D histograms can be normalized to the global maximum throughout the series.
Synchrotron sources. Bending magnet, wiggler, undulator and elliptic undulator are calculated internally within xrt. Please look the section Comparison of synchrotron source codes for the comparison with other popular codes. If the photon source is one of the synchrotron sources, the total flux in the beam is reported not just in number of rays but in physical units of ph/s. The total power or absorbed power can be opted instead of flux and is reported in W. The power density can be visualized by isolines. The magnetic gap of undulators can be tapered. Undulators can be calculated in near field. Custom magnetic field is also possible. Undulators can be calculated on GPU, with a high gain in computation speed, which is important for tapering and near field calculations.
Shapes. There are several predefined shapes of optical elements implemented as python classes. The python inheritance mechanism simplifies creation of other shapes: the user specifies methods for surface height and surface normal. The surface and the normal are defined either in local Cartesian coordinates or in user-defined parametric coordinates. Parametric representation enables closed shapes such as capillaries or wave guides. It also enables exact solutions for complex shapes (e.g. a logarithmic spiral or an ellipsoid) without any expansion. The methods of finding the intersections of rays with the surface are very robust and can cope with pathological cases such as sharp surface kinks. Notice that the search for intersection points does not involve any approximation and has only numerical inaccuracy which is set by default as 1 fm. Any surface can be combined with a (differently and variably oriented) crystal structure and/or (variable) grating vector. Surfaces can be faceted.
Energy dispersive elements. Implemented are crystals in dynamical diffraction, gratings (also with efficiency calculations), Fresnel zone plates, Bragg-Fresnel optics and multilayers in dynamical diffraction. Crystals can work in Bragg or Laue cases, in reflection or in transmission. The two-field polarization phenomena are fully preserved, also within the Darwin diffraction plateau, thus enabling the ray tracing of crystal-based phase retarders.
Materials. The material properties are incorporated using three different tabulations of the scattering factors, with differently wide and differently dense energy meshes. Refractive index and absorption coefficient are calculated from the scattering factors. Two-surface bodies, such as plates or refractive lenses, are treated with both refraction and absorption.
Multiple reflections. xrt can trace multiple reflections in a single optical element. This is useful, for example in ‘whispering gallery’ optics or in Montel or Wolter mirrors.
Non-sequential optics. xrt can trace non-sequential optics where different parts of the incoming beam meet different surfaces. Examples of such optics are poly-capillaries and Wolter mirrors.
Singular optics. xrt correctly propagates vortex beams, which can be used for studying the creation of vortex beams by transmissive or reflective optics.
Global coordinate system. The optical elements are positioned in a global coordinate system. This is convenient for modeling a real synchrotron beamline. The coordinates in this system can be directly taken from a CAD library. The optical surfaces are defined in their local systems for the user’s convenience.
Beam categories. xrt discriminates rays by several categories: good, out, over and dead. This distinction simplifies the adjustment of entrance and exit slits. An alarm is triggered if the fraction of dead rays exceeds a specified level.
Portability. xrt runs on Windows and Unix-like platforms, wherever you can run python.
Examples. xrt comes with many examples. See the galleries, the links are at the top bar.
Tests. xrt comes with the tests of its several main aspects. See the folder tests and also this documentation for the tests of undulators, optical elements, materials and wave propagation.

Dependencies¶

numpy, scipy and matplotlib are required. If you use OpenCL for calculations on GPU or CPU, you need AMD/NVIDIA drivers, Intel CPU only OpenCL runtime (these are search key words), pytools and pyopencl. PyQt4 or PyQt5 are needed for xrtQook. Spyder (as library of Spyder IDE) is highly recommended for nicer view of xrtQook. OpenGL is required for xrtGlow.

See detailed instructions for installing dependencies.

Get xrt¶

xrt is available as source distribution from pypi.python.org and from GitHub. The distribution archive also includes tests and examples. The complete documentation is available online at Read the Docs and offline as zip file at GitHub.

Installation¶

Unzip the .zip file into a suitable directory and use sys.path.append(path-to-xrt) in your script. You can also install xrt to the standard location by running python setup.py install from the directory where you have unzipped the archive, which is less convenient if you try different versions of xrt and/or different versions of python. Note that python-64-bit is by ~20% faster than the 32-bit version (tested with WinPython). Also consider Speed tests for making choice between Windows and Linux and between Python 2 and Python 3.

Get help¶

For getting help and/or reporting a bug please use GitHub xrt Issues.

Citing xrt¶

Please cite xrt as: K. Klementiev and R. Chernikov, “Powerful scriptable ray tracing package xrt”, Proc. SPIE 9209, Advances in Computational Methods for X-Ray Optics III, 92090A; doi:10.1117/12.2061400. Online documentation at xrt.readthedocs.io; doi:10.5281/zenodo.1252468.

Acknowledgments¶

Josep Nicolás and Jordi Juanhuix (synchrotron Alba) are acknowledged for discussion and for their Matlab codes used as examples at early stages of the project.
Andrew Geondzhian and Victoria Kabanova (summer students at DESY) are acknowledged for their help in coding the classes of synchrotron sources.
Rami Sankari and Alexei Preobrajenski (MAX IV Laboratory) are thanked for discussion, testing and comparing with external codes.
Hasan Yavaş, Jozef Bednarčik, Dmitri Novikov and Alexey Zozulya (DESY Photon Science) are acknowledged for supplied cases.
Hamed Tarawneh (MAX IV Laboratory) has initiated the custom field undulator project and supplied tables of magnetic field.
Bernard Kozioziemski (LLNL) has pointed to the importance of considering the total absorption cross-section (not just photoelectric absorption) in the cases of light materials at high energy. He also did deep tests of the mosaic crystal model.
Emilio Heredia (CLS) is thanked for valuable suggestions on GUI.
Tim May (CLS) is acknowledged for inspiring the custom magnetic field calculation development, providing field tables and comparing the results with other codes.
Louisa Pickworth (MAX IV Laboratory) has provided curves from other codes for testing reflectivity/transmittivity of multilayers.

User’s Guide¶

Using xrt¶

Scripting in python¶

You need to prepare a script that gives instructions on how to get the wanted ray properties and prepare the graphs. The scripting is different for different backends (backend is a module or an external program that supplies ray distributions). Currently, xrt supports two backends: raycing – an internal backend – and shadow.

Please consider the supplied examples.

Running a python script¶

You can run it from your shell as python your-script.py or from within an IDE, e.g. Spyder.

Interacting with the plots¶

matplotlib provides an interactive navigation toolbar that allows zoom, pan and save functions.

Saving the results¶

You can save a plot using its navigation toolbar or by specifying the parameter saveName of the plot in your Python script.

You can also save any of the plot attributes through a script. For this, use the parameter afterScript of xrt.runner.run_ray_tracing to point to a function that will be executed after all the repeats have been completed. A typical usage of such a function is to pickle 1D and 2D histograms, fwhm values etc. For example:

plot.xaxis.total1D is the 1D histogram of xaxis of plot. Use an appropriate axis, there are three available: xaxis, yaxis and caxis.
plot.xaxis.binCenters and axis.binEdges are the edges and centers of xaxis bins
plot.total2D is the 2D histogram
plot.total4D is the 4D mutual intensity, see fluxKind for the appropriate options
plot.dx, plot.dy and plot.dE store fwhm’s of the three 1D histograms
flux or power are accessed as plot.flux and plot.power

Using xrtQook for script generation¶

Start xrtQook: type python xrtQookStart.pyw from xrt/gui or, if you have installed xrt by running setup.py, type xrtQookStart.pyw from any location.

Note

The scaling of GUI may behave differently in different systems and Qt versions. If it looks wrong, try the scaling options in xrtQookStart.pyw.

Note

If you want to start xrtQook from Spyder, select the run option “Execute in an external system terminal”.
Rename beamLine to myTestBeamline by double clicking on it (you do not have to, only for demonstration).
Right-click on myTestBeamline and Add Source → BendingMagnet. The same can be done from the icon buttons on the left.
In its properties change eMin to 10000-10 and eMax to 10000+10. The middle of this range will be used to automatically align crystals (one crystal in this example) unless the parameter myTestBeamline.alignE sets another value. Blue color indicates non-default values. These will be included into the generated script. All the default-valued parameters do not propagate into the script.
Create a crystalline material CrystalSi. This will create a Si111 crystal at room temperature.
Add a generic OE -> OE. This will add an optical element with a flat surface.

Note

The sequence of the inserted optical elements does matter! This sequence determines the order of beam propagation.
In its properties select the created crystal as ‘material’, put [0, 20000, 0] as ‘center’ (i.e. 20 m from source) and “auto” (with or without quotes) as ‘pitch’.
Add a screen to the beamline.
Give it [0, 21000, auto] as ‘center’. Its height – the last coordinate – will be automatically calculated from the previous elements.
Check the beamline layout with xrtGlow.
Add a plot and select the local screen beam.
Define an offset to the color (energy) axis.
Save the beamline layout as xml.
Generate python script (the button with a code page and the python logo), save the script and run it.
In the console output you can read the actual pitch (Bragg angle) for the crystal and the screen position.

Using custom optical elements in xrtQook¶

Create custom classes of optical elements in a module located in a writable folder. Copy xrt/gui/xrtQookStart.py starter script to that folder, edit it and run it. Editing is very simple: uncomment a few lines at the top of the script that import the custom module and enable integration with xrtQook.

Notes on using xrtGlow¶

3D glasses button is a two-state button. When it is pressed, xrtGlow will update its view whenever changes are made to the beamline in xrtQook. If you close the window of xrtGlow and the button is pressed, xrtGlow will pop up again after any change in xrtQook. To really close xrtGlow, deactivate the button.
You can use the CTRL-F1 shortcut to open xrtGlow, F4 to dock/undock it into/from the xrtQook window.
The Navigation panel of xrtGlow has several columns. The last columns may be hidden in the initial view. You can access them by enlarging the window.
The element will appear on the Navigation panel only if interacts with the beam, i.e. has an assigned method returning Beams.
Export to Image is available under the context menu -> File. You can save and load the scene settings (camera position, model orientation, rays opacity and so on) as well.
Load the example …/examples/withRaycing/_QookBeamlines/lens3.xml and follow the instructions in Description tab in order to understand the visualization precision vs the swiftness of the 3D manipulations.
From xrtGlow, press F1 to see the available keyboard shortcuts. Also observe the available pop-up menu by right mouse click.
Movements of the model are separated for the transverse plane and longitudinal direction. Use SHIFT-MouseLeft and ALT-MouseLeft for corresponding movements.
The color histogram without Virtual Screen shows the color map – the correspondence between the selected physical parameter (e.g. energy) and the colors. With Virtual Screen active (by F3), the plot shows a histogram of the selected parameter as distributed on Virtual Screen. In both cases the user may select a sub-band on the color plot by the mouse. The vertical extent in that selection is irrelevant.
Virtual Screen is instantiated on the Beam as close as possible to the center of the window. There are several ways to move it:
1. Holding CTRL-MouseLeft: moves the Virtual Screen along the Beam.
2. CTRL-SHIFT-MouseLeft and CTRL-ALT-MouseLeft: moves the whole beamline through the fixed Virtual Screen in transverse or longitudinal directions correspondingly.
If color gradients overlap on the Virtual Screen it can be useful to expand the color axis in real space by enabling the Color Bump. Do not forget that the resulting height distribution is artificial, does not present the real intersections of rays and is only used for convenience.

Rays or footprints visualisation can be enabled/disabled either by setting corresponding checkboxes in the Navigation Panel for individual elements or globally by changing the opacity of the lines and points in the Color Panel. The same applies for the Projections.
Intensity cut-off allows to omit the visualisation of the darkest/weakest rays. It is especially important if Intensity defines the Value key in HSV color space when dark rays can shadow the whole beam.
Convenient way to inspect the detailed beam footprint on the coordinate grid is to use Projections: disable the Perspective, select only the footprint of interest on the Navigation Panel (or disable all and just leave the Virtual Screen on), enable the projection, set to zero the Projection Line Opacity (or Line Width, it will do the job too), increase the Projection Point Opacity to improve the visibility, enable the Fine Grid. Increase the number of rays in the source if necessary.

If you have any doubts regarding the orientation of the optical element or trying to identify the directions, you can plot local coordinate axes by checking the corresponding option on the Scene panel or in the context menu. Make sure that the surface rendering is enabled for this element on the Navigation panel. Orientation of the diffraction planes will be represented by the yellow arrow in case of the crystals with asymmetric cut.

Depth test is disabled by default for the Points. Enable it if you do not want the footprints to shine through solid surfaces of the optical elements. Be aware that the Points may be obscured by rays in this case.
Scene checkbox ‘Virtual Screen for Indexing’ can be used to filter the rays hitting the Virtual Screen. This is convenient for retrospective analysis, to highlight the rays of the initial beam that reach the final point.

Antialiasing can improve the visual quality of the scene, but it seriously affects the performance (depending on the number of rays / elements in the model), only enable it after all modifications to the scene are applied, prior the Export to file. Nevertheless antialiasing is always enabled for the coordinate grid.
Default Zoom does not involve the coordinate grid, if you want to Zoom In/Out the whole scene, use CTRL-MouseWheel.

Using xrt as x-ray calculator¶

xrt can be used as a library for calculations of synchrotron sources and material properties related to x-ray scattering, diffraction and propagation: reflectivity, transmittivity, refractive index, absorption coefficient etc.

See the scripts in \examples\withRaycing\00_xRayCalculator\. Each script consists of:

imports, with possibly defining a path to xrt if it is installed in a non-standard location,
a definition of a source or a material object,
invocation of a corresponding method of interest for a specified range of energy, angles etc. and
plotting.

Example 1a: Undulator

import numpy as np
import matplotlib.pyplot as plt
import sys; sys.path.append(r"c:\Ray-tracing")
import xrt.backends.raycing.sources as rs

source = rs.Undulator(eE=3.0, eI=0.5, eEpsilonX=0.263, eEpsilonZ=0.008,
    betaX=9.539, betaZ=1.982, period=18.5, n=108, K=0.52, distE='BW')

energy = np.linspace(3850, 4150, 601)
theta = np.linspace(-1, 1, 51) * 30e-6
psi = np.linspace(-1, 1, 51) * 30e-6
I0, l1, l2, l3 = source.intensities_on_mesh(energy, theta, psi)
dtheta, dpsi = theta[1] - theta[0], psi[1] - psi[0]
flux = I0.sum(axis=(1, 2)) * dtheta * dpsi

plt.plot(energy, flux)
plt.show()

Example 1b: Undulator, tuning curves

Use the method:

tunesE, tunesF = source.tuning_curves(energy, theta, psi, harmonics, Ks)

See the script \examples\withRaycing\00_xRayCalculator\calc_undulator_tune.py.

Example 2a: Crystal reflectivity

import numpy as np
import matplotlib.pyplot as plt
import sys; sys.path.append(r"c:\Ray-tracing")
import xrt.backends.raycing.materials as rm

crystal = rm.CrystalSi(hkl=(1, 1, 1))

E = 9000
dtheta = np.linspace(-20, 80, 501)
theta = crystal.get_Bragg_angle(E) + dtheta*1e-6
curS, curP = crystal.get_amplitude(E, np.sin(theta))

plt.plot(dtheta, abs(curS)**2, 'r', dtheta, abs(curP)**2, 'b')
plt.show()

Example 2b: Crystal reflectivity: Single and double crystal

See the script \examples\withRaycing\00_xRayCalculator\calc_crystal_rocking_curve.py.

Example 3: Multilayer reflectivity

mSi = rm.Material('Si', rho=2.33)
mW = rm.Material('W', rho=19.3)
mL = rm.Multilayer(mSi, 27, mW, 18, 40, mSi)

E = 10000
theta = np.linspace(0, 2.0, 1001)  # degrees
rs, rp = mL.get_amplitude(E, np.sin(np.deg2rad(theta)))[0:2]

plt.plot(theta, abs(rs)**2, 'r', theta, abs(rp)**2, 'b')
plt.show()

Example 4: Mirror reflectivity

stripe = rm.Material(('Si', 'O'), quantities=(1, 2), rho=2.2)

E = np.logspace(1 + np.log10(3), 4 + np.log10(5), 501)
theta = 7e-3
rs, rp = stripe.get_amplitude(E, np.sin(theta))[0:2]

plt.semilogx(E, abs(rs)**2, 'r', E, abs(rp)**2, 'b')
plt.gca().set_xlim(E[0], E[-1])
plt.show()

Example 5: Material absorption

mat = rm.Material('Be', rho=1.848)

E = np.logspace(1, 4 + np.log10(3), 501)
mu = mat.get_absorption_coefficient(E)

plt.loglog(E, mu)
plt.gca().set_xlim(E[0], E[-1])
plt.show()

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 right-handed):

The global coordinate system. It is arbitrary (user-defined) with one requirement driven by code simplification: Z-axis is vertical. For example, the system origin of Alba synchrotron is in the center of the ring at the ground level with Y-axis 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.
1. 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.
2. of an optical element. The origin is on the optical surface. Z is out-of-surface. 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 user-supplied functions for the surface height (z) and the normal as functions of (x, y) are defined in the local oe system.
3. 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 module-level 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¶

Sampling strategies¶

There are three main strategies for creating samples of optical sources, depending on study purpose, implemented in xrt.

1. Monte Carlo sampling by intensity, or ray generation. This is the default mode of all synchrotron sources in xrt and is the best sampling strategy for ray tracing or hybrid wave propagation (when the wave part starts at a downstream optical element). Technically it is done by the acceptance-rejection method. In this method, the optical field is calculated in a 3D space – energy and two transverse angles – for a set of uniformly random values of these three variables within the predefined ranges. The transverse angles also get normally distributed angular shifts within the emittance distribution, independently for each ray. The acceptance random variable is uniformly sampled in the range from zero to the intensity maximum for each field sample, and the sample is accepted if the acceptance variable is below the sample intensity, otherwise it is rejected. The sampling continues until the requested number of samples is reached. The resulted space density of the samples (rays) is proportional to the field intensity. Simultaneously, this sampling gives the value of total flux as the product of the acceptance ratio times the 3D sample volume times the intensity maximum.

2. Uniform Monte Carlo sampling, or wave generation. This mode is necessary for wave propagation when it starts right from the source. The 3D calculation space – energy and two transverse angles – is sampled uniformly. This way of sampling is still possible to use for ray tracing, and it is even faster at the source as no samples are rejected, but it is less efficient at the downstream part of the beamline as a significant part of the samples (or even a majority of them) is of low intensity. The main usage pattern of this sampling is, however, for single electron (or macro electrons) field propagation (enabled by the source option filamentBeam=True), when all wave samples are attributed to the same electron with a single displacement within emittance distribution and a single shift of gamma (relative electron energy) within energy spread distribution. The flux is a sum of all intensities times the 3D sample volume over the number of samples.

3. Grid sampling. This is frequently a quick method for studying synchrotron sources in reciprocal space and energy that must be used with care as it is full of pitfalls. Grid samples can be obtained by intensities_on_mesh() method of the synchrotron source. Three aspects must be considered when making a grid. (i) For a sharp field distribution, e.g. for an undulator at zero emittance and zero electron beam energy spread, the sharp field rings may interfere with the grid (form moiré patterns) and result in false fringes in spectral flux density. In such cases, the grid has to be very finely spaced. (ii) When the electron beam emittance is non-zero, the calculated field values have to be convolved with the electron beam divergence; this convolution is done inside intensities_on_mesh(). Due to the edge effects of the convolution, the grid borders of the thickness of the doubled electron beam divergence must be discarded from the field arrays returned by intensities_on_mesh() and therefore these grid margins have to be added beforehand. (iii) The number of energy spread samples must be set sufficient for the used energy spread value, see the end of the next paragraph.

Electron beam energy spread is treated differently in the above three sampling strategies. In ray sampling, every ray gets its gamma sample (relative electron energy) of the normal energy spread distribution. In wave sampling, energy spread is sampled once for all wave samples when in filamentBeam regime and individually per wave sample otherwise. In grid sampling, a new grid dimension is introduced internally that contains samples of gamma within energy spread distribution with subsequent averaging over that dimension. Notice that a wide energy spread distribution (that can be a case for linacs) may require more energy spread samples (the parameter eSpreadNSamples of intensities_on_mesh()) than the default number 36 that is usually good for energy spread sigma of 0.1%. Because of this new dimension, the grid may become too large to fit into RAM, which might force the user to invoke the method intensities_on_mesh() in a loop that scans over individual photon energies.

Electron beam real space size is also treated differently in the above three sampling strategies. In ray sampling, every ray gets its transverse origin shift in the source plane within emittance distribution. In wave sampling, the wave source position is sampled once for all wave samples when in filamentBeam regime and individually per wave sample otherwise. In grid sampling, real space electron beam size is totally ignored.

See the example Undulator radiation through rectangular aperture, where several calculation methods are compared.

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.

concatenate(beam)¶: Adds beam to self. Useful when more than one source is presented.

export_beam(fileName, fformat='npy')¶: Saves the beam to a binary file. File format can be Numpy ‘npy’, Matlab ‘mat’ or python ‘pickle’. Matlab format should not be used for future imports in xrt as it does not allow correct load.

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, roll=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:: Linear (distx, disty, distz) and angular (distxprime, distzprime) source distributions. Accepted values: ‘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: Linear (dx, dy, dz) and angular (dxprime, dzprime) source sizes. 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: 1-D array-like: 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 ray-tracing. 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, roll, yaw: float: rotation angles around x, y and z axes. 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: 1-D array-like: 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, roll=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, roll=0, yaw=0)¶

bl: instance of BeamLine

name: str

center: tuple of 3 floats: 3D point in global system
w0: float or 2-sequence: Gaussian beam waist size. If a 2-sequence, 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: 1-D array-like: 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, roll, yaw: float: rotation angles around x, y and z axes. Useful for canted sources.

class xrt.backends.raycing.sources.LaguerreGaussianBeam(*args, **kwargs)¶

Implements Laguerre-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__(*args, **kwargs)¶

vortex: None or tuple(l, p): specifies Laguerre-Gaussian beam with l the azimuthal index and p >= 0 the radial index.

class xrt.backends.raycing.sources.HermiteGaussianBeam(*args, **kwargs)¶

Implements Hermite-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__(*args, **kwargs)¶

TEM: None or tuple(m, n): specifies Hermite-Gaussian 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}\]

user-defined (polarization is 4-sequence 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.copy_beam(beamTo, beamFrom, indarr, includeState=False, includeJspEsp=True)¶: Copies arrays of beamFrom to arrays of beamTo. The slicing of the arrays is given by indarr.

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 of run_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 z-axis to be directed along the beamline in order to be compatible with the cited works. Elsewhere in xrt z-axis 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, ST-M-91-12B, July 1991, Presented at 4th Int. Conf. on Synchrotron Radiation Instrumentation, Chester, England, Jul 15-19, 1991

[WS]

R. J. Dejus, Computer simulations of the wiggler spectrum, Nucl. Instrum. Methods Phys. Res., Sect. A, 347 (1994) 56-60

[XOP] (1,2,3,4,5,6)

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/data-analysis/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 3-dimensional 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 far-field 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 non-equivalent 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 Runge-Kutta fourth-order method.

For the Undulator and custom field models we directly calculate the integral using the Clenshaw-Curtis quadrature, it proves to converge as quickly as the previously used Gauss-Legendre 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 on-axis/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 re-evaluated 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 n-th 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, 22-27 Apr 1996: proceedings (CERN. Geneva, 1998) 129-190.

[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: User-specified 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.

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: User-specified 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 center-to-screen 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 ray-tracing.
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.

intensities_on_mesh(energy='auto', theta='auto', psi='auto', harmonic=None, eSpreadSigmas=3.5, eSpreadNSamples=36, mode='constant', resultKind='Stokes')¶: Returns Stokes parameters (as a 4-list of arrays, when resultKind == ‘Stokes’) or intensities and OAM (Orbital Angular Momentum) matrix elements (as [Is, Ip, OAMs, OAMp, Es, Ep], when resultKind == ‘vortex’)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 arrays is a 1D array of an individually selectable length. Energy spread is sampled by a normal distribution and the resulting field values are averaged over it. eSpreadSigmas is sigma value of the distribution; eSpreadNSamples sets the number of samples. The resulted transverse field is convolved with angular spread by means of scipy.ndimage.filters.gaussian_filter. mode is a synonymous parameter of that filter that controls its behaviour at the borders.

Note

This method provides incoherent averaging over angular and energy spread of electron beam. The photon beam phase is lost here!

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 “macro-electrons” 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 macro-electron 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 and Undulator. 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 1e-6 and below. Custom field calculation may require setting the precision of 1e-3.

gNodes: int

Number of integration nodes in each of the integration intervals. If not provided at init, will be defined automatically.

targetOpenCL: None, str, 2-tuple or tuple of 2-tuples

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() and platform.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 vendor-specific 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.

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 2-tuple 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 2-tuple 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 5-dimensional. 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 by numpy.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 non-zero 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 mesh-based codes interpolate the source intensity pre-calculated 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) 1221-8.

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 off-axis 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.

Single-electron and multi-electron undulator radiation¶

Here we compare single-electron and multi-electron (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, 22-26 June 1998, p.1177-1179.

	SRW	xrt
single electron
single electron	execution time 974 s	execution time 17.4 s
non-zero emittance
non-zero emittance	execution time 65501 s (sic)	execution time 18.6 s
non-zero emittance, non-zero energy spread
non-zero emittance, non-zero energy spread	execution time 66180 s (sic)	execution time 216 s

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) 596-603.

[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 off-axis radiation (greenish color) but not to the on-axis radiation (bluish color).

Undulator radiation in near field¶

Notice that on the following pictures the p-polarized 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, p-polarized
near field at 5 m, full flux
near field at 5 m, p-polarized

at 25 m	SPECTRA	xrt
far field at 25 m, full flux
far field at 25 m, p-polarized
near field at 25 m full flux
near field at 25 m, p-polarized

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 non-zero 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 s-p 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

User-specified name, occasionally used for diagnostics output.

center: 3-sequence 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() or get_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() and xyz_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 self-contained in the above mentioned user-supplied functions. For example, these can be viewed as cylindrical-like 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. Class SurfaceOfRevolution gives an example of the transformation functions and represents a useful kind of parametric surface. The methods local_n() (surface normal) and local_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, ρ₀, P₀, P₁, …], where ρ₀ is the constant line density in inverse mm, P₀ – 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, 1e-6, 3.1e-7] 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 Bragg-Fresnel 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, 2-tuple or tuple of 2-tuples

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.

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 3-tuple 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 6-sequence: 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 3-tuple or a 6-tuple. 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 or RoundAperture. nrays: if int, specifies the number of randomly distributed samples the surface within self.limPhysX limits; if 2-tuple 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)¶: Ground-bent (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)¶: Ground-2D-bent (Johansson) reflective optical element.

class xrt.backends.raycing.oes.GeneralBraggToroid(JohannToroid)¶: Ground-2D-bent 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). This element supports volumetric diffraction model, if corresponding parameter is enabled in the assigned material.

__init__(*args, **kwargs)¶

R: float or 2-tuple.: 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.BentLaue2D(OE)¶

Parabolically bent reflective optical element in Laue geometry. Meridional and sagittal radii (Rm, Rs) can be defined independently and have same or opposite sign, representing concave (+, +), convex (-, -) or saddle (+, -) shaped profile. This element supports volumetric diffraction model, if corresponding parameter is enabled in the assigned material.

__init__(*args, **kwargs)¶

Rm: float or 2-tuple.: Meridional bending radius.
Rs: float or 2-tuple.: Sagittal radius.

class xrt.backends.raycing.oes.GroundBentLaueCylinder(BentLaueCylinder)¶: Ground-bent 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 (3-sequences). 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.

Values of the ellipse’s semi-major and semi-minor axes lengths can be accessed after init at ellipseA and ellipseB respectively.

Note

Any of p, q, f1, f2 or pAxis can be set as instance attributes of this mirror object; the ellipsoid parameters parameters will be recalculated automatically.

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 (3-sequence) 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.

Note

Any of p, q, f1, f2 or parabolaAxis can be set as instance attributes of this mirror object; the ellipsoid parameters parameters will be recalculated automatically.

The usage is exemplified in test_param_mirror.py.

class xrt.backends.raycing.oes.ConicalMirror(OE)¶

Conical mirror with its base parallel to the side of the cone.

__init__(*args, **kwargs)¶

L0: float: Distance from the center of the mirror to the vertex of the cone. This distance is measured along the surface, NOT along the axis.
theta: float: Opening angle of the cone (axis to surface) in radians.

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 N-th local beam returns the absorbed intensity on N-th 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. It 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: Thickness in mm.
wedgeAngle: float: Relative angular misorientation of the back plane.

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 multi-surface 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 curved-flat 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.ParaboloidCapillaryMirror(SurfaceOfRevolution)¶

Paraboloid of revolution a.k.a. Mirror Lens. By default will be oriented for focusing. Set yaw to 180deg for collimation.

__init__(*args, **kwargs)¶

q: float: Distance from the center of the element to focus.
r0: float: Radius at the center of the element.

class xrt.backends.raycing.oes.EllipsoidCapillaryMirror(SurfaceOfRevolution)¶

Ellipsoid of revolution a.k.a. Mirror Lens. Do not forget to set reasonable limPhysY.

__init__(*args, **kwargs)¶

ellipseA: float: Semi-major axis, half of source-to-sample distance.
ellipseB: float: Semi-minor axis. Do not confuse with the size of the actual capillary!
workingDistance: float: Distance between the end face of the capillary and focus. Mind the length of the optical element for proper positioning.

class xrt.backends.raycing.oes.NormalFZP(OE)¶

Implements a circular Fresnel Zone Plate, as it is described in X-Ray 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 4-sequence or str: The two foci given by 3-sequences 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 property kind='grating' but rather kind='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 class BlazedGrating implements an ad hoc method find_intersection() for selecting the 1st intersection point among the several possible ones. The left picture below illustrates the behavior of OE (the footprint shown by circles is partially in the shadowed area). The right picture demonstrates the correct behavior of BlazedGrating 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.

OE

BlazedGrating

__init__(*args, **kwargs)¶

blaze, antiblaze: float: Angles in radians
rho: float: Constant line density in inverse mm. If the density is variable, use gratingDensity from the parental class OE with the 1st argument ‘y’ (i.e. along y-axis).

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.

__init__(*args, **kwargs)¶

aspect: float: Top-to-period ratio of the groove.
depth: float: Depth of the groove in mm.

For the VLS density, use gratingDensity of the parental class OE with the 1st argument ‘y’ (i.e. along y-axis).

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

Tests of Optical elements¶

See the tests here.

Materials¶

Module materials defines atomic and material properties related to x-ray 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[1E-12cm] IncoherentX-section[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) 185-200, http://www-phys.llnl.gov/Research/scattering/RTAB.html

[Waasmaier]

D. Waasmaier & A. Kirfel, Acta Cryst. A51 (1995) 416-413

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. See also predefined materials in modules materials_compounds and materials_elemental.

__init__(elements=None, quantities=None, kind='auto', rho=0, t=None, table='Chantler total', efficiency=None, efficiencyFile=None, name='', refractiveIndex=None)¶

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

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 cross-sections. The tabulation by Chantler can optionally have total absorption cross-sections. 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, X-ray interactions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, Z=1-92, Atomic Data and Nuclear Data Tables 54 (no.2) (1993) 181-342.

[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) 71-643.

[BrCo]

http://www.bmsc.washington.edu/scatter/periodic-table.html ftp://ftpa.aps.anl.gov/pub/cross-section_codes/ S. Brennan and P.L. Cowan, A suite of programs for calculating x-ray absorption, reflection and diffraction performance for a variety of materials at arbitrary wavelengths, Rev. Sci. Instrum. 63 (1992) 850-853.

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 (zero-based) or a column number (also zero-based, 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.
refractiveIndex: float or complex or numpy array or str: Explicitly defines the refractive index of the material. Can be used for the energy ranges not covered by the tables of scattering factors (IR, visible). Can be set as: a) float or complex value, for a constant, energy-independent refractive index. b) a 3-column numpy array containing Energy in eV, real and imaginary parts of the complex refractive index c) filename for an *.xls or CSV table containig same columns as a numpy array in b).

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], [Als-Nielsen]. 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 .

[Als-Nielsen] (1,2)

Jens Als-Nielsen, Des McMorrow, Elements of Modern X-ray 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,j-1}$ of the interdiffusion/roughness interface.
substThickness: float in Å: Is only relevant in transmission if substrate is present.
geom: str: Either ‘transmitted’ or ‘reflected’.

get_amplitude(E, beamInDotNormal, x=None, y=None, ucl=None)¶

Calculates amplitude of reflectivity [Als-Nielsen]. 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 ($N-1$) is:

\[R_N = \frac{r_{N-1, N} + r_{N, N+1} p_N^2} {1 + r_{N-1, 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_{j-1, j} + R_{j+1} p_j^2} {1 + r_{j-1, 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_{j-1,z}\sigma^{2}_{j,j-1})$, where $k_{j,z}$ is longitudinal component of the wave vector in j-th layer [Nevot-Croce].

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_{j-1, j}p_j}{1 + r_{j-1, 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}}}\]

[Nevot-Croce]

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 period-averaged 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 Å.

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='', volumetricDiffraction=False, useTT=False, nu=None, mosaicity=0)¶

hkl: sequence: hkl indices.
d: float: Interatomic spacing in Å.
V: float: Unit cell volume in Å³. If not given, is calculated from d assuming a cubic symmetry.
factDW: float: Debye-Waller 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.

volumetricDiffraction: bool: By default the diffracted ray originates in the point of incidence on the surface of the crystal in both Bragg and Laue case. When volumetricDiffraction is enabled, the point of diffraction is generated randomly along the transmitted beam path, effectively broadening the meridional beam profile in plain Laue crystal. If the crystal is bent, local deformation of the diffracting plane is taken into account, creating the polychromatic focusing effect.

useTT: bool

Specifies whether the reflectivity is calculated by analytical formula (useTT is False) or by solution of the Takagi-Taupin equations (when useTT is True). The latter case is based on PyTTE code [PyTTE1] [PyTTE2] that was adapted to running the calculations on GPUs.

[PyTTE1] (1,2)

https://github.com/aripekka/pyTTE

[PyTTE2]

A.-P. Honkanen, S. Huotari, IUCrJ 8 (2021) 102-115. doi:10.1107/S2052252520014165

Warning

You need a good graphics card to run these calculations! The corresponding optical element, that utilizes the present crystal material class, must specify targetOpenCL (typically, ‘auto’) and precisionOpenCL (in Bragg cases ‘float32’ is typically sufficient and ‘float64’ is typically needed in Laue cases).

nu: float

Poisson’s ratio. Can be used for calculation of reflectivity in bent isotropic crystals with [PyTTE1]. Not required for plain crystals or for crystals with predefined compliance matrix, see materials_crystals. If provided, overrides existing compliance matrix.

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 of the new ray origins and (iii) the reflectivity is calculated following the work [BaconLowde].

In the above stage (ii), the impact points are sampled according to the secondary extinction distribution. For a thin crystal (when t is specified) those rays that go over the crystal thickness retain the original incoming direction and their x, y, z coordinates (the ray heads) are put on the back crystal surface. These rays are also attenuated by the mosaic crystal. Note again that the impact points do not lie on the front crystal surface, so to plot them in a 2D XY plot becomes useless. You can still plot them as YZ or XZ to study the secondary extinction depth. The remaining rays (those sampled within the crystal depth) get reflected and also attenuated depending on the penetration depth.

Note

The amplitude calculation in the mosaic case is implemented only in the reflection geometry. The transmitted beam can still be studied by ray-tracing as we split the beam in our modeling of secondary extinction, see the previous paragraph.

Note

The mosaicity is assumed large compared with the Darwin width. Therefore, there is no continuous transition mosaic-to-perfect 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 p-polarizations ($\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(1-b)}{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 x-ray 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 cross-section) but with flux.

[BD]

V. A. Belyakov and V. E. Dmitrienko, Polarization phenomena in x-ray 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_amplitude_pytte(E, beamInDotNormal, beamOutDotNormal=None, beamInDotHNormal=None, xd=None, yd=None, alphaAsym=None, Ry=None, Rx=None, ucl=None, tolerance=1e-06, maxSteps=10000000.0, autoLimits=True, signal=None)¶

Calculates complex amplitude reflectivity for s- and p-polarizations ($\gamma = s, p$) in Bragg and Laue cases, based on modified PyTTE code

alphaAsymm: float: Angle of asymmetry in radians.
Ry: float: Meridional radius of curvature in mm. Positive for concave bend.
Rx: float: Sagittal radius of curvature in mm. Positive for concave bend.
ucl:: instance of XRT_CL class, defines the OpenCL device and precision of calculation. Calculations should run fine in single precision, float32. See XRT_CL.
tolerance: float: Precision tolerance for RK adaptive step algorithm.
maxSteps: int: Emergency exit to avoid kernel freezing if the step gets too small.
autoLimits: bool: If True, the algorithm will try to automatically determine the angular range where reflectivity will be calculated by numeric integration. Useful for ray-tracing applications where angle of incidence can be too far from Bragg condition, and integration might take unnesessarily long time.

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, X-Ray Optics High-Energy-Resolution Applications, Springer-Verlag 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 diamond-like 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 d-spacing 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. See also predefined crystals in module materials_crystals.

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 3-sequences: 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 Takagi-Taupin 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 Takagi-Taupin 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: 2-list 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.

Predefined Materials¶

The module materials_crystals contains predefined classes for most commonly used crystals. Lattice parameters, atomic positions and references have been semi-automatically parsed from XOP/DABAX [XOP] Crystals.dat. To use a crystal in a script simply import the module and instantiate its class:

import xrt.backends.raycing.materials_crystals as xcryst
myInSbXtal = xcryst.InSb(hkl=(3, 1, 1))  # default hkl=(1, 1, 1)

The crystals inherit from Crystal and can use its methods to calculate diffraction amplitudes, the Darwin width, extinction depth etc.

The following crystal classes are defined in this module (sorted by cell volume in Å³), marked in bold are those with available elastic constants:

Be (16.23), Graphite (34.98), Titanium (35.32), Diamond (45.38), Iron (46.31), Copper (47.24), Platinum (60.43), LiF (65.27), Aluminum (66.41), Gold (67.83), LaB6 (71.61), LaB6NIST (71.83), SiC (82.2), AlphaQuartz (113), Si (160.2), GaP (161.9), NaCl (179.4), GaAs (180.7), Ge (181.1), InP (202.1), CsF (216.9), InAs (219.9), GaSb (226.5), KCl (249.2), Sapphire (254.7), AlphaAlumina (254.8), InSb (271.9), LiNbO3 (318.2), PET (324.8), CsCl (345.9), Beryl (657.3), TlAP (855.9), KAP (858.9), RbAP (862.9), KTP (871.3), Muscovite (934.2)

The module materials_compounds contains predefined classes for compound materials found at the CXRO Table of Densities of Common Materials. Air composition and density are given as in [Cox].

[Cox]

Cox, Arthur N., ed. (2000), Allen’s Astrophysical Quantities (Fourth ed.), AIP Press, pp. 258–259, ISBN 0-387-98746-0

To use a compound material in a script, simply import the module and instantiate its class:

import xrt.backends.raycing.materials_compounds as xcomp
kapton = xcomp.Polyimide()

The compound materials inherit from Material and can use its methods to calculate reflection or transmission amplitudes, absorption coefficient, refractive index etc.

Note

The compound materials do not provide crystal diffraction amplitudes even if they occur naturally as crystals. To calculate diffraction on crystals please use materials_crystals.

The following compound classes are defined in this module (sorted by chemical formula):

Vacuum ( ), SilverBromide (AgBr), Sapphire (Al2O3), AluminumArsenide (AlAs), AluminumPhosphide (AlP), BoronOxide (B2O3), BoronCarbide (B4C), BoronNitride (BN), BerylliumOxide (BeO), CVDDiamond (C), Mylar (C10H8O4), Polycarbonate (C16H14O3), Kimfol (C16H14O3), Polyimide (C22H10N2O5), Teflon (C2F4), Polypropylene (C3H6), PMMA (C5H8O2), ParyleneC (C8H7Cl), ParyleneN (C8H8), Fluorite (CaF2), CadmiumSulfide (CdS), CadmiumTelluride (CdTe), CadmiumTungstate (CdWO4), CobaltSilicide (CoSi2), Cromium3Oxide (Cr2O3), CesiumIodide (CsI), CopperIodide (CuI), GalliumArsenide (GaAs), GalliumNitride (GaN), GalliumPhosphide (GaP), Water (H2O), HafniumOxide (HfO2), Indium3Oxide (In2O3), IndiumNitride (InN), IndiumAntimonide (InSb), IridiumOxide (IrO2), Mica (KAl3Si3O12H2), LithiumFluoride (LiF), LithiumHydride (LiH), LithiumHydroxide (LiOH), MagnesiumSilicide (Mg2Si), MagnesiumFluoride (MgF2), MagnesiumOxide (MgO), Manganese2Oxide (MnO), Manganese4Oxide (MnO2), Molybdenum4Oxide (MoO2), Molybdenum6Oxide (MoO3), MolybdenumSilicide (MoSi2), Air (N0.781O0.209Ar0.009), RockSalt (NaCl), NiobiumNitride (NbN), NiobiumSilicide (NbSi2), NickelSilicide (Ni2Si), NickelOxide (NiO), RutheniumSilicide (Ru2Si3), Ruthenium4Oxide (RuO2), Zerodur (Si.56Al.5P.16Li.04Ti.02Zr.02Zn.03O2.46), ULEGlass (Si.925Ti.075O2), SiliconNitride (Si3N4), SiliconCarbide (SiC), Silica (SiO2), Quartz (SiO2), TantalumOxide (Ta2O5), TantalumSilicide (Ta2Si), TantalumNitride (TaN), TitaniumNitride (TiN), Rutile (TiO2), TitaniumSilicide (TiSi2), Uranium4Oxide (UO2), VanadiumNitride (VN), TungstenCarbide (WC), ZinсOxide (ZnO), ZinсSulfide (ZnS), ZirconiumNitride (ZrN), Zirconia (ZrO2), ZirconiumSilicide (ZrSi2)

The module materials_elemental contains predefined classes for elemental materials. Most atomic densities have been adopted from the NIST table of X-Ray Mass Attenuation Coefficients. Densities for Fr and At given according to [Lavrukhina_Pozdnyakov].

[Lavrukhina_Pozdnyakov]

Lavrukhina, A. K. and Pozdnyakov, A. A. (1970). Analytical Chemistry of Technetium, Promethium, Astatine, and Francium. Translated by R. Kondor. Ann Arbor–Humphrey Science Publishers. p. 269. ISBN 978-0-250-39923-9

Note

Densities are given for the phase state at ambient conditions, e.g. Nitrogen as N₂ gas.

To use an elemental material in a script simply import the module and instantiate its class:

import xrt.backends.raycing.materials_elemental as xmat
nitrogenGas = xmat.N()

The elemental materials inherit from Material and can use its methods to calculate reflection or transmission amplitudes, absorption coefficient, refractive index etc.

Note

The elemental materials do not provide crystal diffraction amplitudes even if they occur naturally as crystals. To calculate diffraction on crystals please use materials_crystals.

The following elemental classes are defined in this module:

H, He, Li, Be, B, C, N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ga, Ge, As, Se, Br, Kr, Rb, Sr, Y, Zr, Nb, Mo, Tc, Ru, Rh, Pd, Ag, Cd, In, Sn, Sb, Te, I, Xe, Cs, Ba, La, Ce, Pr, Nd, Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu, Hf, Ta, W, Re, Os, Ir, Pt, Au, Hg, Tl, Pb, Bi, Po, At, Rn, Fr, Ra, Ac, Th, Pa, U

Tests of Materials¶

See the tests here.

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 front-end 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: User-specified name, can be used for diagnostics output.
center: 3-sequence 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: 3-tuples 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. Both x and z can also be set as instance attributes.
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.

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 or RoundAperture. nrays of samples are randomly distributed over the slit area.

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: 3-tuples 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. Both x and z can also be set as instance attributes.

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: User-specified name, can be used for diagnostics output.
center: 3-sequence of floats: 3D point in global system.
opening: sequence: Coordinates [(x0, y0),…(xN, yN)] of the polygon vertices.
x, z: 3-tuples 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. Both x and z can also be set as instance attributes.
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.

class xrt.backends.raycing.apertures.GridAperture¶

Implements a grid of rectangular apertures. See tests/raycing/test_polygonal_aperture.py

__init__(bl=None, name='', center=[0, 0, 0], x='auto', z='auto', alarmLevel=None, dx=0.5, dz=0.5, px=1.0, pz=1.0, nx=7, nz=7)¶

dx and dz: float: Opening sizes (horizontal and vertical).
px and pz: float: Pitch steps (periods).
nx and nz: int: The number of grid openings, counted from the central hole in -ve and +ve directions, making (2*nx+1)*(2*nz+1) rectangular holes in total.

class xrt.backends.raycing.apertures.SiemensStar¶

Implements a Siemens Star pattern. See tests/raycing/test_polygonal_aperture.py

__init__(bl=None, name='', center=[0, 0, 0], x='auto', z='auto', alarmLevel=None, nSpokes=9, r=0, rx=0, rz=0, phi0=0, vortex=0, vortexNradial=7)¶

nSpokes: int: The number of spoke openings.
r or (rx and rz): float: The radius of spokes or ellipse semiaxes.
phi0: float: The angle of rotation of the star. At 0 (default), the center of the top spoke is vertical.
vortex: float: Lets the spokes bend. At =1, they bend by one spoke position.
vortexNradial: int: Used with non-zero vortex. The number of segments in the bent spokes.

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: User-specified name, can be used for diagnostics output.
center: tuple of 3 floats: 3D point in the global system.
x, z: 3-sequence 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. Both x and z can also be set as instance attributes.
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.

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 or RoundAperture. 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 by meshgrid.

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]

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: 3-tuples 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. Both x and z can also be set as instance attributes.
R: float: Radius of the hemisphere in mm.

Wave propagation (diffraction)¶

Time dependent diffraction¶

We start from the Kirchhoff integral theorem in the general (time-dependent) 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 back-Fourier transformed to the frequency domain represented by a new re-sampled 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 time-dependent 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 non-monochromaticity? We repeat the sequence of ray-tracing 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 ray-tracing, 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 X-ray optics simulation: combining geometric ray-tracing and wavefront propagation, J. Synchrotron Rad. 21 (2014) 669–678.

[shadowGuide]

Cerrina, “SHADOW User’s Guide” (1998).

[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 Monte-Carlo 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 Monte-Carlo 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 pre-exponent factor is neglected. The new wave directions are thus given by a Kirchhoff-like 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/en-us/library/windows/hardware/ff569918(v=vs.85).aspx

Tip

If you plan to utilize xrt on a multi-GPU platform under Linux, we recommend KDE based systems (e.g. Kubuntu). It is enough to install the package fglrx-updates-dev 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 class Beam, 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 by prepare_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 by diffract).

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 previous prepare_wave, which can be reported to reflect by noIntersectionSearch=True. Here, reflect takes the beam right before the surface and propagates it to right after it. As a result of such a zero-length travel, the wave gets no additional propagation phase but only a complex-valued 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⁴ 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 or RoundAperture. nrays: if int, specifies the number of randomly distributed samples the surface within self.limPhysX limits; if 2-tuple 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 or RoundAperture. nrays of samples are randomly distributed over the slit area.

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 or RoundAperture. 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 by meshgrid.

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 (another common name is cross-spectral density) 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_1-x_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_1-x_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 from scipy.linalg, where a feasible image size should not exceed ~100×100 pixels (i.e. ~10⁸ 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.
1. Out of r repeats of $E(x, y)$ build a stacked data matrix $D$ with N_x×N_y rows and r columns.
2. 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 covariance matrix $D^{+}D$ of the size r².
3. The spectra of eigenvalues of matrices $DD^{+}$ and $D^{+}D$ are equal, plus zeroes for the bigger matrix.
4. Their eigenvectors (being eigenmodes $V_i$ of $DD^{+}$ and principal component axes $v_i$ of $D^{+}D$) corresponding to the same eigenvalue are transformed to each other with a factor $D$ or $D^{+}$: $\tilde{V}_i = Dv_i$ and $\tilde{v}_i = D^{+}V_i$, after which transformation they must be additionally normalized ($\tilde{v}_i$ and $\tilde{V}_i$ are non-normalized eigenvectors) [the proof is a one line matrix equation].
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]:

\[{\rm DoTC} = \frac{\iiiint |J_{12}|^2 dx_1 dy_1 dx_2 dy_2} {\left[\iint J_{11} dx_1 dy_1\right]^2}\]

[Saldin2008]

E.L. Saldin, E.A. Schneidmiller, M.V. Yurkov, Coherence properties of the radiation from X-ray free electron laser, Opt. Commun. 281 (2008) 1179–88.

[Vartanyants2010]

I.A. Vartanyants and A. Singer, Coherence properties of hard x-ray synchrotron sources and x-ray free-electron 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:

1. DoTC = Tr(J²)/Tr²(J).
2. 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 complex-shaped 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 non-zero, 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 2-tuple, 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() or calc_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 non-monochromaticity 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 non-zero 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 inter-optics 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 macro-electrons) 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 front-end 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 macro-electrons¶

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 ~1e-16 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 macro-electrons” test; in usual studies one would typically need a few thousand macro-electrons 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 macro-electrons¶

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

Tests of wave propagation¶

See the tests here.

Shadow backend¶

This is a deprecated backend. raycing is much more functional. Module shadow works with shadow input files, starts the ray-tracing and gets its output.

Warning

Shadow3 is not supported.

Warning

xrtQook and xrtGlow do not work with this backend. A beamline created in Shadow cannot be visualized by xrtGlow.

Description of shadow¶

… can be found in the manual pages of your shadow distribution. In connection with using shadow in xrt, it is important to understand the naming of the output beam files and the naming of parameters in the setup files.

Preparation for a shadow run¶

Note

on shadow under Windows Vista and 7:

Under Windows Vista and 7 shadow does not work out of the box because of epath (a part of shadow) reporting an error. There is a workaround consisting of simply stopping the Windows’ Error Reporting Service.

Create a folder where you will store your ray-tracing script and the output images. Make it as Python’s working directory. Create there a sub-folder tmp0. Put there your shadow project file along with all necessary data files (reflectivities, crystal parameters etc.). Run shadow and make sure it correctly produces output files you want to accumulate (like star.01).

Now you need to generate two command files that run shadow source and shadow trace. These are system-specific and also differ for different shadow sources. Under Windows, this can be done as follows: set the working directory of shadowVUI as your tmp0, run Source in shadowVUI and rename the produced shadowvui.bat to shadow-source.bat; then run Trace and rename the produced shadowvui.bat to shadow-trace.bat.

Try to run the generated command files in order to check their validity.

If you want to use multi-threading then copy tmp0 to tmp1, tmp2 etc. (in total as many directories as you have threads).

Scripting in python¶

The simplest script consists of 4 lines:

import xrt.runner as xrtr
import xrt.plotter as xrtp
plot1 = xrtp.XYCPlot('star.01')
xrtr.run_ray_tracing(plot1, repeats=40, updateEvery=2, threads=1)

Modifying input files¶

There are two types of input files in shadow:

Of ‘Namelist’ or ‘GFile’ type (both are in terminology of shadow). These are parameter files which consist of lines field = value. Examples of this type are: start.xx and end.xx. Such files describe optical elements and two sources: geometrical and bending magnet.
Of a non-named type consisting of lines of values, one value per line. Examples are xsh_input_source_tmp.inp and xsh_nphoton_tmp.inp. Such files describe two other sources: wiggler and undulator.

If you want to run a series of ray tracing studies for variable physical or geometrical parameters (e.g. for a variable meridional radius of a focusing mirror), you have to find out which parameter in the shadow’s input files controls the desired variable. The only way for this is to play with the parameter in a GUI (I use shadowVUI) and look for changes in the corresponding start.xx text file. Once you have discovered the needed parameter, you can change it in your Python script. There are two functions for this:

xrt.backends.shadow.modify_input(fileNameList, *editlines)¶

modifies a shadow text input file (like start.NN) which consists of lines field = value.

Parameters:

fileNameList: str or list of str: A list of file names is useful when executing shadow in parallel in several directories.

editlines: list of tuples of strings (field, value).

Returns:

0 if successful, otherwise -1.

Example:

>>> modify_input('start.00',('istar1',str(seed)))  #change seed

xrt.backends.shadow.modify_xsh_input(fileNameList, *editlines)¶

modifies a shadow xsh text input file (like xsh_nphoton_tmp.inp) which consist of lines of values, one value per line.

Parameters:

fileNameList: str or list of str: A list of file names is useful when executing shadow in parallel in several directories.
editlines: list of tuples (fieldNo: int, value: str): fieldNo is zero-based index of the modified parameter.

Returns:

0 if successful, otherwise -1.

Example:

>>> modify_xsh_input('xsh_nphoton_tmp.inp', (2, energyRange[0]),
                     (3, energyRange[1]))

The 1st parameter in the above functions is a simple file name or a list of file names. If you run with several threads or processes then you must modify all the versions of the shadow input file in the directories tmp0, tmp1 etc. The following function helps in doing this:

xrt.backends.shadow.files_in_tmp_subdirs(fileName, processes=1)¶

Creates and returns a list of full file names of copies of a given file located in the process directories. This list is needed for reading and writing to several versions of one file (one for each process) in one go. Useful in user’s scripts.

Example:

>>> start01 = shadow.files_in_tmp_subdirs('start.01', processes=4)
>>> shadow.modify_input(start01, ('THICK(1)', str(thick * 1e-4)))

Writing a loop generator¶

A sequence of ray tracing runs is controlled by a generator (a function that returns by yield) which modifies the shadow input files, optionally specifies information text panels, define the output file names etc. in a loop. The same generator is used for normalization, if requested, when it is (quickly) run in the second pass.

Consider an example:

import xrt.runner as xrtr
import xrt.plotter as xrtp
import xrt.backends.shadow as shadow
plot1 = xrtp.XYCPlot('star.01') #create a plot
textPanel = plot1.fig.text(0.88, 0.8, '', transform=plot1.fig.transFigure,
  size=14, color='r', ha='center') #create a text field, see matplotlib help
threads = 2
start01 = shadow.files_in_tmp_subdirs('start.01', threads)
def plot_generator():
    for thick in [0, 60, 400, 1000, 1500]: #thickness in um
        shadow.modify_input(start01, ('THICK(1)', str(thick * 1e-4)))
        filename = 'filt%04imum' %thick #output file name without extension
        plot1.saveName = [filename + '.pdf', filename + '.png']
        textPanel.set_text('filter\nthickness\n%s $\mu$m' %thick)
        yield

xrtr.run_ray_tracing(plot1, repeats=40, generator=plot_generator,
  threads=threads, globalNorm=True)

Customizing your plots¶

Module plotter provides classes describing axes and plots, as well as containers for the accumulated arrays (histograms) for subsequent pickling/unpickling or for global flux normalization. The module defines several constants for default plot positions and sizes. The user may want to modify them in the module or externally as in the xrt_logo.py example.

Note

Each plot has a 2D positional histogram, two 1D positional histograms and, typically, a 1D color histogram (e.g. energy).

Warning

The two 1D positional histograms are not calculated from the 2D one!

In other words, the 1D histograms only respect their corresponding limits and not the other dimension’s limits. There can be situations when the 2D image is black because the screen is misplaced but one 1D histogram may still show a beam distribution if in that direction the screen is positioned correctly. This was the reason why the 1D histograms were designed not to be directly dependent on the 2D one – this feature facilitates the troubleshooting of misalignments. On the other hand, this behavior may lead to confusion if a part of the 2D distribution is outside of the visible 2D area. In such cases one or two 1D histograms may show a wider distribution than the one visible on the 2D image. For correcting this behavior, one can mask the beam by apertures or by selecting the physical or optical limits of an optical element.

Tip

If you do not want to create plot windows (e.g. when they are too many or when you run xrt on a remote machine) but only want to save plots, you can use a non-interactive matplotlib backend such as Agg (for PNGs), PDF, SVG or PS:

matplotlib.use('agg')

Importantly, this must be done at the very top of your script, right after import matplotlib and before importing anything else.

The plots can be customized at the time of their creation using many optional parameters of xrt.plotter.XYCPlot.__init__() constructor.

class xrt.plotter.XYCPlot(beam=None, rayFlag=(1,), xaxis=None, yaxis=None, caxis=None, aspect='equal', xPos=1, yPos=1, ePos=1, title='', invertColorMap=False, negative=False, fluxKind='total', fluxUnit='auto', fluxFormatStr='auto', contourLevels=None, contourColors=None, contourFmt='%.1f', contourFactor=1.0, saveName=None, persistentName=None, oe=None, raycingParam=0, beamState=None, beamC=None, useQtWidget=False)¶

Container for the accumulated histograms. Besides giving the beam images, this class provides with useful fields like dx, dy, dE (FWHM), cx, cy, cE (centers) and intensity which can be used in scripts for producing scan-like results.

__init__(beam=None, rayFlag=(1,), xaxis=None, yaxis=None, caxis=None, aspect='equal', xPos=1, yPos=1, ePos=1, title='', invertColorMap=False, negative=False, fluxKind='total', fluxUnit='auto', fluxFormatStr='auto', contourLevels=None, contourColors=None, contourFmt='%.1f', contourFactor=1.0, saveName=None, persistentName=None, oe=None, raycingParam=0, beamState=None, beamC=None, useQtWidget=False)¶

beam: str

The beam to be visualized.

In raycing backend:: The key in the dictionary returned by run_process(). The values of that dictionary are beams (instances of Beam).
In shadow backend:: The Shadow output file (star.NN, mirr.NN` or screen.NNMM). It will also appear in the window caption unless title parameter overrides it.

This parameter is used for the automatic determination of the backend in use with the corresponding meaning of the next two parameters. If beam contains a dot, shadow backend is assumed. Otherwise raycing backend is assumed.

rayFlag: int or tuple of ints

shadow: 0=lost rays, 1=good rays, 2=all rays. raycing: a tuple of integer ray states: 1=good, 2=out, 3=over, 4=alive (good + out), -NN = dead at oe number NN (numbering starts with 1).

xaxis, yaxis, caxis: instance of XYCAxis or None.

If None, a default axis is created. If caxis=’category’ and the backend is raycing, then the coloring is given by ray category, the color axis histogram is not displayed and ePos is ignored.

Warning

The axes contain arrays for the accumulation of histograms. If you create the axes outside of the plot constructor then make sure that these are not used for another plot. Otherwise the histograms will be overwritten!

aspect: str or float

Aspect ratio of the 2D histogram, = ‘equal’, ‘auto’ or numeric value (=x/y). aspect =1 is the same as aspect =’equal’.

xPos, yPos: int

If non-zero, the corresponding 1D histograms are visible.

ePos: int

Flag for specifying the positioning of the color axis histogram:

ePos =1: at the right (default, as usually the diffraction plane is vertical)
ePos =2: at the top (for horizontal diffraction plane)
ePos =0: no color axis histogram

title: str

If non-empty, this string will appear in the window caption, otherwise the beam will be used for this.

invertColorMap: bool

Inverts colors in the HSV color map; seen differently, this is a 0.5 circular shift in the color map space. This inversion is useful in combination with negative in order to keep the same energy coloring both for black and for white images.

negative: bool

Useful for printing in order to save black inks. See also invertColorMap.

=False: black bknd for on-screen presentation
=True: white bknd for paper printing

The following table demonstrates the combinations of invertColorMap and negative:

	invertColorMap =False	invertColorMap =True
negative =False
negative =True

Note that negative inverts only the colors of the graphs, not the white global background. Use a common graphical editor to invert the whole picture after doing negative=True:

(such a picture would nicely look on a black journal cover, e.g. on that of Journal of Synchrotron Radiation ;) )

fluxKind: str

Can begin with ‘s’, ‘p’, ‘+-45’, ‘left-right’, ‘total’, ‘power’, ‘Es’, ‘Ep’ and ‘E’. Specifies what kind of flux to use for the brightness of 2D and for the height of 1D histograms. If it ends with ‘log’, the flux scale is logarithmic.

If starts with ‘E’ then the field amplitude or mutual intensity is considered, not the usual intensity, and accumulated in the 2D histogram or in a 3D stack:

If ends with ‘xx’ or ‘zz’, the corresponding 2D cuts of mutual intensity are accumulated in the main 2D array (the one visible as a 2D histogram). The plot must have equal axes.
If ends with ‘4D’, the complete mutual intensity is calculated and stored in plot.total4D with the shape (xaxis.bins*yaxis.bins, xaxis.bins*yaxis.bins).

Warning

Be cautious with the size of the mutual intensity object, it is four-dimensional!

If ends with ‘PCA’, the field images are stored in plot.field3D with the shape (repeats, xaxis.bins, yaxis.bins) for further Principal Component Analysis.
If without these endings, the field aplitudes are simply summed in the 2D histogram.

fluxUnit: ‘auto’ or None

If a synchrotron source is used and fluxUnit is ‘auto’, the flux will be displayed as ‘ph/s’ or ‘W’ (if fluxKind == ‘power’). Otherwise the flux is a unitless number of rays times transmittivity | reflectivity.

fluxFormatStr: str

Format string for representing the flux or power. You can use a representation with powers of ten by utilizing ‘p’ as format specifier, e.g. ‘%.2p’.

contourLevels: sequence

A sequence of levels on the 2D image for drawing the contours, in [0, 1] range. If None, the contours are not drawn.

contourColors: sequence or color

A sequence of colors corresponding to contourLevels. A single color value is applied to all the contours. If None, the colors are automatic.

contourFmt: str

Python format string for contour values.

contourFactor: float

Is applied to the levels and is useful in combination with contourFmt, e.g. contourFmt = r’%.1f mW/mm$^2$’, contourFactor = 1e3.

saveName: str or list of str or None

Save file name(s). The file type(s) are given by extensions: png, ps, svg, pdf. Typically, saveName is set outside of the constructor. For example:

filename = 'filt%04imum' %thick #without extension
plot1.saveName = [filename + '.pdf', filename + '.png']

persistentName: str or None

File name for reading and storing the accumulated histograms and other ancillary data. Ray tracing will resume the histogramming from the state when the persistent file was written. If the file does not exist yet, the histograms are initialized to zeros. The persistent file is rewritten when ray tracing is completed and the number of repeats > 0.

Warning

Be careful when you use it: if you intend to start from zeros, make sure that this option is switched off or the pickle files do not exist! Otherwise you do resume, not really start anew.

if persistentName ends with ‘.mat’, a Matlab file is generated.

oe: instance of an optical element or None

If supplied, the rectangular or circular areas of the optical surfaces or physical surfaces, if the optical surfaces are not specified, will be overdrawn. Useful with raycing backend for footprint images.

raycingParam: int

Used together with the oe parameter above for drawing footprint envelopes. If =2, the limits of the second crystal of DCM are taken for drawing the envelope; if =1000, all facets of a diced crystal are displayed.

beamState: str

Used in raycing backend. If not None, gives another beam that determines the state (good, lost etc.) instead of the state given by beam. This may be used to visualize the incoming beam but use the states of the outgoing beam, so that you see how the beam upstream of the optical element will be masked by it. See the examples for capillaries.

beamC: str

The same as beamState but refers to colors (when not of ‘category’ type).

Each of the 3 axes: xaxis, yaxis and caxis can be customized using many optional parameters of xrt.plotter.XYCAxis.__init__() constructor.

class xrt.plotter.XYCAxis(label='', unit='mm', factor=None, data='auto', limits=None, offset=0, bins=128, ppb=2, density='histogram', invertAxis=False, outline=0.5, fwhmFormatStr='%.1f')¶

Contains a generic record structure describing each of the 3 axes: X, Y and Color (typ. Energy).

__init__(label='', unit='mm', factor=None, data='auto', limits=None, offset=0, bins=128, ppb=2, density='histogram', invertAxis=False, outline=0.5, fwhmFormatStr='%.1f')¶

label: str

The label of the axis without unit. This label will appear in the axis caption and in the FWHM label.

unit: str

The unit of the axis which will follow the label in parentheses and appear in the FWHM value

factor: float

Useful in order to match your axis units with the units of the ray tracing backend. For instance, the shadow length unit is cm. If you want to display the positions as mm: factor=10; if you want to display energy as keV: factor=1e-3. Another usage of factor is to bring the coordinates of the ray tracing backend to the world coordinates. For instance, z-axis in shadow is directed off the OE surface. If the OE is faced upside down, z is directed downwards. In order to display it upside, set minus to factor. if not specified, factor will default to a value that depends on unit. See def auto_assign_factor().

data: int for shadow, otherwise array-like or function object

shadow:

zero-based index of columns in the shadow binary files:

0	x
1	y
2	z
3	x’
4	y’
5	z’
6	Ex s polariz
7	Ey s polariz
8	Ez s polariz
9	lost ray flag
10	photon energy
11	ray index
12	optical path
13	phase (s polarization)
14	phase (p polarization)
15	x component of the electromagnetic vector (p polar)
16	y component of the electromagnetic vector (p polar)
17	z component of the electromagnetic vector (p polar)
18	empty

raycing:

use the following functions (in the table below) or pass your own one. See raycing for more functions, e.g. for the polarization properties. Alternatively, you may pass an array of the length of the beam arrays.

x	raycing.get_x
y	raycing.get_y
z	raycing.get_z
x’	raycing.get_xprime
z’	raycing.get_zprime
energy	raycing.get_energy

If data = ‘auto’ then label is searched for “x”, “y”, “z”, “x’”, “z’”, “energy” and if one of them is found, data is assigned to the listed above index or function. In raycing backend the automatic assignment is additionally implemented for label containing ‘degree (for degree of polarization)’, ‘circular’ (for circular polarization rate), ‘path’, ‘incid’ or ‘theta’ (for incident angle), ‘order’ (for grating diffraction order), ‘s’, ‘phi’, ‘r’ or ‘s’ (for parametric representation of OE).

limits: 2-list of floats [min, max]

Axis limits. If None, the limits are taken as np.min and np.max for the corresponding array acquired after the 1st ray tracing run. If limits == ‘symmetric’, the limits are forced to be symmetric about the origin. Can also be set outside of the constructor as, e.g.:

plot1.xaxis.limits = [-15, 15]

offset: float

An offset value subtracted from the axis tick labels to be displayed separately. It is useful for the energy axis, where the band width is most frequently much smaller than the central value. Ignored for x and y axes.

no offset	non-zero offset

bins: int

Number of bins in the corresponding 1D and 2D histograms. See also ppb parameter.

ppb: int

Screen-pixel-per-bin value. The graph arrangement was optimized for bins * ppb = 256. If your bins and ppb give a very different product, the graphs may look ugly (disproportional) with overlapping tick labels.

density: ‘histogram’ or ‘kde’

The way the sample density is calculated: by histogram or by kde [KDE].

invertAxis: bool

Inverts the axis direction. Useful for energy axis in energy- dispersive images in order to match the colors of the energy histogram with the colors of the 2D histogram.

outline: float within [0, 1]

Specifies the minimum brightness of the outline drawn over the 1D histogram. The maximum brightness equals 1 at the maximum of the 1D histogram.

=0	=0.5	=1

fwhmFormatStr: str

Python format string for the FWHM value, e.g. ‘%.2f’. if None, the FWHM value is not displayed.

Customizing the look of the graphs¶

The elements of the graphs (fonts, line thickness etc.) can be modified by editing the matplotlibrc file.

Alternatively, you can do this “dynamically”, e.g.:

import matplotlib as mpl
mpl.rcParams['font.family'] = 'serif'
mpl.rcParams['font.size'] = 12
mpl.rcParams['axes.linewidth'] = 0.5

The sizes and positions of the graphs are controlled by several constants specified in xrt. You can either modify them directly in the module or in your script, like it was done for creating the logo image (see also logo_xrt.py script in the examples):

import xrt.plotter as xrtp
xrtp.xOrigin2d = 4
xrtp.yOrigin2d = 4
xrtp.height1d = 64
xrtp.heightE1d = 64
xrtp.xspace1dtoE1d = 4

Saving the histogram arrays¶

The histograms can be saved to and later restored from a persistent file. Doing so is always recommended as this can save you time afterwards if you decide to change a label or a font size or to make the graphs negative: you can do this without re-doing ray tracing. Remember, however, that a new run of the same script will initialize the histograms not from zero but from the last saved state. You will notice this by the displayed number of rays. If you want to initialize from zero, just delete the persistent file. Note that you can run run_ray_tracing() with repeats=0 in order to re-display the saved plots without starting ray tracing.

Controlling ray tracing execution¶

The function run_ray_tracing() has many options on the program flow.

xrt.runner.run_ray_tracing(plots=[], repeats=1, updateEvery=1, pickleEvery=None, energyRange=None, backend='raycing', beamLine=None, threads=1, processes=1, generator=None, generatorArgs=[], generatorKWargs='auto', globalNorm=0, afterScript=None, afterScriptArgs=[], afterScriptKWargs={})¶

This function is the entry point of xrt. Parameters are all optional except the 1st one. Please use them as keyword arguments because the list of parameters may change in future versions.

plots: instance of XYCPlot or a sequence of
instances or an empty sequence if no graphical output is wanted.

repeats: int
The number of ray tracing runs. It should be stressed that accumulated are not rays, which would be limited by the physical memory, but rather the histograms from each run are summed up. In this way the number of rays is unlimited.

updateEvery: int
Redrawing rate. Redrawing happens when the current iteration index is divisible by updateEvery.

pickleEvery: int
Saving rate. Applicable to plots with a defined persistentName. If None, the pickling will happen once at the end.

energyRange: [eMin: float, eMax: float]
Only in shadow backend: If not None, sets the energy range of shadow source. Alternatively, this can be done directly inside the generator.

backend: str
so far supported: {‘shadow’ | ‘raycing’ | ‘dummy’}

beamLine: instance of BeamLine, used
with raycing backend.

threads, processes: int or str
The number of parallel threads or processes, should not be greater than the number of cores in your computer, otherwise it gives no gain. The bigger of the two will be used as a signal for using either threading or multiprocessing. If they are equal, threading is used. See also performance tests. If ‘all’ is given then the number returned by multiprocessing.cpu_count() will be used.

Warning

You cannot use multiprocessing in combination with OpenCL because the resources (CPU or GPU) are already shared by OpenCL. You will get an error if processes > 1. You can still use threads > 1 but with a little gain.

Note

For the shadow backend you must create tmp0, tmp1 etc. directories (counted by threads or processes) in your working directory. Even if the execution is not parallelized, there must be tmp0 with the shadow files prepared in it.

generator: generator object
A generator for running complex ray-tracing studies. It must modify the optics, specify the graph limits, define the output file names etc. in a loop and return to xrt by yield. See the supplied examples.

generatorArgs, generatorKWargs: list and (dictionary or ‘auto’)
If generatorKWargs is ‘auto’, the following keyword dictionary will be used for the generator: kwargs = {} if generator is defined within the caller of run_ray_tracing() or if generatorArgs is not empty, otherwise kwargs = {‘plots’=pots, ‘beamLine’=beamLine}.

globalNorm: bool
If True, the intensity of the histograms will be normalized to the global maximum throughout the series of graphs. There are two flavors of normalization:

only the heights of 1D histograms are globally normalized while the brightness is kept with the normalization to the local maximum (i.e. the maximum in the given graph).

both the heights of 1D histograms and the brightness of 1D and 2D histograms are globally normalized.

The second way is physically more correct but sometimes is less visual: some of the normalized pictures may become too dark, e.g. when you compare focused and strongly unfocused images. Both normalizations are saved with suffixes _norm1 and _norm2 for you to select the better one.

Here is a normalization example where the intensity maximum was found throughout a series of images for filters of different thickness. The brightest image was for the case of no filter (not shown here) and the normalization shown below was done relative to that image:

normalized to local maximum

global normalization, type 1

global normalization, type 2

afterScript: function object
This function is executed at the end of the current script. For example, it may run the next ray-tracing script.

afterScriptArgs, afterScriptKWargs: list and dictionary
args and kwargs for afterScript.

Using multiprocessing and multithreading¶

Module multipro defines BackendProcess as a subclass of multiprocessing.Process or threading.Thread. You can opt between deriving from multiprocessing or threading by selecting the corresponding parameter in run_ray_tracing(). The multiprocessing is normally faster than multithreading but has an inconvenience when the user aborts the execution: the processes have to be killed manually.

See the performance tests.

Calculations on GPU¶

GPU can be used for several types of calculations, e.g. of an undulator source. The gain is enormous. You should try it. Even on an “ordinary” Intel processor the execution with OpenCL becomes significantly faster.

Here are some benchmarks on a system with Intel Core i7-3930K 3.20 GHz CPU, ASUS Radeon R9 290 GPU, Python 2.7.6 64-bit, Windows 7 64-bit. Script examples/withRaycing/01_SynchrotronSources/synchrotronSources.py, 1 million rays, execution times in seconds:

CPU 1 process	5172, 1 CPU process loaded
CPU 10 processes	1245, with heavily loaded system
openCL with CPU	163, with highly loaded system
openCL with GPU	132, with almost idle CPU

You will need AMD/NVIDIA drivers (if you have a GPU, however this is not a must), a CPU only OpenCL runtime, pytools and pyopencl.

Note

When using OpenCL, no further parallelization is possible by means of multithreading or multiprocessing. You should turn them off by using the default values processes=1 and threads=1 in the run properties.

Please run the script tests/raycing/info_opencl.py for getting information about your OpenCL platforms and devices. You will pass then the proper indices in the lists of the platforms and devices as parameters to pyopencl methods or, alternatively, pass ‘auto’ to targetOpenCL.

Important

Consider the warnings and tips on using xrt with GPUs.

Hint

Consider also Speed tests for a few selected cases.

OpenCL Server¶

This example contains a Python script designed to run on a GPU server, leveraging ZeroMQ (ZMQ) for efficient data transfer. The script acts as a remote accelerator device, receiving data from a client Python script, performing calculations on the GPU, and returning the results for plotting on a local computer.

Users can seamlessly execute their scripts in their favorite IDE, offloading resource-intensive calculations to a remote server over the network. The only trade-off is the potential delay due to data transfer, which is outweighed by the benefits when local computations take longer than data transfer time. Furthermore, the local graphical user interface (GUI) remains responsive without freezes or issues caused by high GPU/CPU loads. This script now supports all acceleration scenarios:

synchrotron sources,
wave propagation,
bent crystals,
multilayer reflectivity.

Script Components¶

The GPU accelerator script is comprised of two files located at tests/raycing/RemoteOpenCLCalculation:

zmq_server.py: The server script is the main component, responsible for receiving data and getting kernel names from the client. It listens on a predefined port, processes the received package, executes the specified kernel on the GPU and sends the computed data back to the client. This server script can be executed independently or in conjunction with the queue manager.
queue_device.py: The queue manager script facilitates the handling of multiple user requests and the distribution of computational tasks across multiple servers. It provides scalability and load balancing capabilities. The queue manager can be executed on the same machine as the server or on a dedicated node. However, when running the queue manager on a separate node, data transfer times may increase.

Running the Script¶

To execute the GPU accelerator script, follow these steps:

Set up the GPU server environment with the necessary dependencies, including pyzmq, xrt and underlying dependencies (numpy, scipy, matplotlib, pyopencl). Start the server script, either as a standalone process or in conjunction with the queue manager, based on your specific requirements.

Ensure that the client Python script is configured to connect to the correct server (or queue manager) address and port:

targetOpenCL="GPU_SERVER_ADDRESS:15559"

Gallery of plots and scripts 1. Synchrotron sources¶

Synchrotron sources¶

The images below are produced by \tests\raycing\test_sources.py and by \examples\withRaycing\01_SynchrotronSources\synchrotronSources.py.

Bending magnet¶

On a transversal screen the image is unlimited horizontally (practically limited by the front end). The energy distribution is red-shifted for off-plane photons. The polarization is primarily horizontal. The off-plane radiation has non-zero projection to the vertical polarization plane.

source	total flux	horiz. pol. flux	vert. pol. flux
using WS
internal xrt

The off-plane radiation is in fact left and right polarized:

source	circular polarization rate
using WS
internal xrt

The horizontal phase space projected to a transversal plane at the origin is parabolic:

zero electron beam size	σ_x = 49 µm

Multipole wiggler¶

The horizontal image size is determined by the parameter K. The energy distribution is red-shifted for off-plane photons. The polarization is primarily horizontal. The off-plane radiation has non-zero projection to the vertical polarization plane.

source	total flux	horiz. pol. flux	vert. pol. flux
using WS
internal xrt

The horizontal longitudinal cross-section reveals a sinusoidal shape of the source. The horizontal phase space projected to the transversal plane at the origin has individual branches for each pole.

zero electron beam size	σ_x = 49 µm

Undulator¶

The module test_sources has functions for visualization of the angular and energy distributions of the implemented sources in 2D and 3D. This is especially useful for undulators because they have sharp peaks, which requires a proper selection of angular and energy meshes.

The ray traced images of an undulator source (produced by \examples\withRaycing\01_SynchrotronSources\synchrotronSources.py) are feature-rich. The polarization is primarily horizontal. The off-plane radiation has non-zero projection to the vertical polarization plane.

source	total flux	hor. pol. flux	ver. pol. flux	deg. of pol.
using Urgent
internal xrt

Elliptical undulator¶

An elliptical undulator gives circular images with a higher circular polarization rate in the inner rings:

source	total flux	hor. pol. flux	ver. pol. flux
using Urgent
internal xrt

source	deg. of pol.	circular polarization rate
using Urgent
internal xrt

Custom field undulator¶

A custom magnetic field can be specified by an Excel file or a column text file. The example below is based on a table supplied by Hamed Tarawneh [Tarawneh]. The idea of introducing quasi-periodicity is to shift the n-th harmonics down in energy relative to the exact n-fold multiple of the 1st harmonic energy. This trick eliminates higher monochromator harmonics that are situated at the exact n-fold energies, which is a safer solution compared to a gas absorption filter.

Compare the harmonic energies (half-maximum position at the higher energy side) of the 3rd harmonic with the triple energy of the 1st harmonic.

[Tarawneh]

Quasi-periodic undulator field for ARPES beamline at MAX IV 1.5 GeV ring, (2016) unpublished.

Note

The definition of xyz coordinate system differs for the tabulated field and for xrt screens: z is along the beam direction in the tabulation and as a vertical axis in xrt.

	periodic	quasi-periodic
tabulated field
trajectory top view
wide band image and spectrum
1st harmonic image and spectrum
3rd harmonic image and spectrum

For validation, our calculations are compared here with those by Spectra for a particular case — quasi-periodic undulator defined by the same tabulated field, the 3rd harmonic, at E=20.5 eV. Notice again that Spectra provides either a spectrum or a transverse image while xrt can combine both by using colors and brightness. Notice also that on the following pictures the p-polarized flux is only ~3% of the total flux.

	SPECTRA	xrt
total flux
p-pol flux

Undulator radiation through rectangular aperture¶

The images below are produced by \examples\withRaycing\01_SynchrotronSources\fluxThroughAperture.py.

This example illustrates the sampling concepts explained in Section Sampling strategies.

Grid sampling comes in two flavors: with a grid exactly matching the 100×100 µrad² aperture and with the same grid plus extra 30% margins on each side. Adding margins is crucial for doing convolution with electron beam angular spread. Notice that for the sharp field distributions at zero emittance and zero energy spread, the grid has to be set with fine angular steps, otherwise the flux density spectrum gets unreal ripples. One may play with ntheta and npsi values to observe the effect.

In this particular example (MAX IV Linac), eEspread value is four times bigger than a typical value for a storage ring and so the grid here needs a bigger number of gamma (relative electron beam energy) samples to obtain a smooth spectrum. One may play with eSpreadNSamples parameter to observe the effect.

Ray tracing examples illustrate two different ways of field sampling: uniform reciprocal space sampling and sampling by intensity distribution. Energy can also be sampled variously: uniformly within a given range and on an equidistant mesh (energy scan). The latter case also delivers single energy field images, which also has some calculation overhead.

The following cases are sorted in increased complexity (and calculation time) from top to bottom.

Zero emittance, zero energy spread. Here, angular mesh has to be very dense due to very sharp field features.
Non-zero emittance, zero energy spread. Here, extra margins have to be added to the angular mesh in order to be able to convolve with electron beam divergence. In this example, adding 30% margins (the orange curve) was not enough.
Non-zero emittance, non-zero energy spread. Here, the number of energy spread samples has to be adjusted for a large energy spread sigma otherwise the spectrum may have false sharp peaks.

Finally, the standard ray tracing approach (“rays-b”) should be the method of choice as it does not need any parameter adjustment (in contrast to the grid methods). For not extreme cases (not too small aperture or not to big emittance and not too big energy spread) grid methods can be most efficient but anyways need a sanity check by ray tracing.

Orbital Angular Momentum of helical undulator radiation¶

The images below are produced by \examples\withRaycing\01_SynchrotronSources\undulatorVortex.py.

This example calculates flux and Orbital Angular Momentum (OAM) of helical undulator radiation. The calculation is done on a 3D (energy, theta, psi) or 4D (energy, theta, psi, gamma) rectangular mesh for zero and true electron beam energy spread, respectively. To calculate OAM projection along the propagation direction, we first calculate OAM intensity $I^l_s$ for the field $E_s$

\[I^l_s = E_s^* i\left(\psi\frac{\partial}{\partial\theta}- \theta\frac{\partial}{\partial\psi}\right) E_s\]

and similarly for the field $E_p$. OAM intensity is an incoherent distribution, similarly to the “usual” intensity, as the phase information is lost in it, and therefore it is suitable for incoherent accumulation of fields originated from different electrons distributed within emittance and energy spread distributions. After averaging intensity and OAM intensity over electron beam energy spread distribution and convolving with electron beam angular distributions, we obtain vorticity – the averaged normalized OAM value:

\[l_s = \int I^l_s d\theta d\psi \bigg/ \int I_s d\theta d\psi\]

and similarly for the field $E_p$.

The images below were calculated at the maximum flux energies at the 1st to 4th harmonics. As a reminder, the flux is maximized at energies slightly below the formal harmonic energies; at these energies the transverse distribution is of a donut shape. As a second reminder, the radiation from a helical undulator has only the first harmonic on the undulator axis but also has higher harmonics at finite observation angles.

The coloring of the transverse images is done by field phase. For this purpose, the radiation field (here, $E_s$) was coherently averaged over electron beam energy spread distribution and electron beam angular distributions. This operation is physically illegal but was used here only for the visual effect. The measurable values – intensity, flux and vorticity – were determined by the correct incoherent averaging. The image brightness represents intensity.

	zero emittance, zero energy spread	true emittance, true energy spread
Flux and vorticity
1st harmonic
2nd harmonic
3rd harmonic
4th harmonic

As expected, vorticity approximately equals the harmonic number minus one in the perfect case of zero emittance and zero energy spread. This pattern is quite strongly affected by electron beam energy spread. Emittance plays a much lesser role as the used angular acceptance is much bigger than the electron beam angular distributions (360 µrad vs 5.8 and 2.0 µrad rms).

Gallery of plots and scripts 2. X-ray optics¶

Beamline optics¶

The images below are produced by the scripts in \examples\withRaycing\02_Balder_BL\. The examples show the scans of various optical elements at Balder@MaxIV beamline. The source is a multipole conventional wiggler.

Diamond filter of varying thickness¶

Shown are i) intensity downstream of the filter and its energy spectrum and ii) absorbed power in the filter with its energy spectrum and power density isolines at 85% and 95% of the maximum.

White-beam collimating mirror at varying pitch angle¶

	bare Si	Ir coating
flux downstream of the mirror with its energy spectrum
absorbed power in the mirror with its energy spectrum and power density isolines at 90% of the maximum

Bending of collimating mirror¶

Shown are i) image downstream of the DCM and ii) image at the sample.

Bending of focusing mirror¶

Shown is the image at the sample position.

Both mirrors (collimating + focusing) at varying pitch angle¶

The sagittal radius of the focusing toroid mirror is optimal at 2 mrad pitch angle.

Scanning of Double Crystal Monochromator¶

Scanning of Double Multilayer Monochromator¶

This example shows a scan of a double multilayer monochromator – analog of a double crystal monochromator. On the left is the footprint on the 1st multilayer. On the right is a transversal image of the exit beam. The two multilayers are equal and have 40 pairs of 27-Å-thick silicon on 18-Å-thick tungsten layers on top of a silicon substrate.

Laue Monochromator¶

Bending of a single crystal Laue Monochromator¶

Files in \examples\withRaycing\03_LaueMono

This example shows the reflectivity of a bent 200-µm-thick Si111 Laue crystal at various bending radii and energies. Watch how the band width is growing and the flux is lowering in going to smaller radii.

E = 9 keV	E = 16 keV	E = 25 keV	E = 36 keV

Double bent-crystal Laue monochromator (beam cleaner)¶

This example shows the beam images and rocking curves of a bent 200-µm-thick Si111 double Laue crystal monochromator (similar to the beam cleaner [Karanfil2004]) at various bending radii and energies.

[Karanfil2004]

C. Karanfil, D. Chapman, C. U. Segre and G. Bunker, A device for selecting and rejecting X-ray harmonics in synchrotron radiation beams, J. Synchrotron Rad. 11 (2004) 393-8.

Beam images at various detuning angles of the second crystal, R = 25 m, E ~ 9 keV. Watch how the energy band becomes split and the flux goes down in going away from the parallel positioning (dθ = 0).

Flux vs. detuning angle of the second crystal (rocking curves)

E = 9 keV	E = 16 keV	E = 25 keV	E = 36 keV

Compound Refractive Lenses¶

Files in \examples\withRaycing\04_Lenses

This example demonstrates refraction in x-ray regime. Locus that refracts a collimated beam into a point focus is a paraboloid. The focal distance of such a vacuum-to-solid interface is, as in the usual optics, 2p/δ where p is the focal parameter of the lens paraboloid and δ = 1 - Re(n), n is the refractive index [snigirev]. As for the usual lenses, the diopters of several consecutive lenses are summed up to give the total diopter: $\frac{1}{f} = \frac{1}{f_1} + \frac{1}{f_2} + \ldots$.

This example considers focusing of collimated x-rays of 9 keV at a distance q = 5 m from the lenses. The lenses are double-sided paraboloids (then f = p/δ) with p = 1 mm and zero spacing between the apices of the paraboloids. The thickness of each lens is 2 mm. Their number is an integer number N = round(p/qδ).

The following images demonstrate the focusing along the optical axis close to the nominal focal position for Be and Al CRL’s. The real focal position deviates from the nominal one, where dq = 0, due to the rounding of p/qδ:

This graph shows the relative flux in the focused beam at 9 keV after the given number of double-sided lenses which give approximately equal focal distance of q = 5 m. As seen, low absorbing materials are preferred:

This graph shows the depth of focus as a function of the on-axis coordinate around the nominal focal position. For heavy materials the depth of focus is larger due to the higher absorption of the peripherical rays of the incoming beam. Such lenses act effectively also as apertures thus reducing the focal spot at the expense of flux:

[snigirev]

A. Snigirev, V. Kohn, I. Snigireva, A. Souvorov, and B. Lengeler, Focusing High-Energy X Rays by Compound Refractive Lenses, Applied Optics 37 (1998), 653-62.

Quarter wave plates¶

Files in \examples\withRaycing\05_QWP

Collimated beam, Bragg transmission case¶

This example shows the polarization properties of a collimated linearly polarized beam passed through a diamond plate at various departure from the nominal Bragg angle. The polarization plate is put at 45º to the diffraction plane. Notice that the phase difference between the s- and p-polarized components was calculated here not in the 1-field approximation as elsewhere [Malgrange] that has a pole at $\Delta\theta=0$ but in the general 2-field approximation, see materials.

[Malgrange]

C. Giles, C. Malgrange, J. Goulon, F. de Bergevin, C. Vettier, E. Dartyge, A. Fontaine, C. Giorgetti and S. Pizzini, J. Appl. Cryst. 27 (1994) 232; C. Giles, C. Vettier, F. de Bergevin, C. Malgrange, G. Grübel, and F. Grossi, Rev. Sci. Instrum. 66 (1995) 1518; J. Goulon, C. Malgrange, C. Giles, C. Neumann, A. Rogalev, E. Moguiline, F. De Bergevin and C. Vettier, J. Synchrotron Rad. 3 (1996) 272.

Beam images after the QWP with color axis as 1) energy, 2) circular polarization rate, 3) phase shift between s- and p-components and 4) ratio of axes of the polarization ellipse. Watch how the circular polarization rate becomes close to 1 or -1 at certain departure angles; here, between 16 and 32 arcsec (plus or minus). At the same time also the ratio of axes of the polarization ellipse becomes close to 1 or -1 and with narrow distribution.

E ~ 9 keV, crystal thickness = 200 µm.

Collimated beam, Laue transmission case¶

Same as in the revious subsection but for the Laue case. The thickness of the crystal was selected as to give the path length similar to that in the Bragg case.

E ~ 9 keV, crystal thickness = 500 µm.

Convergent beam, Laue transmission case¶

In this example the beam is focused by a toroidal mirror. Due to the large angular variation in the beam (3 mrad at the position of the QWP), the resulting circular polarization rate is low. The initial polarization is horizontal and the diffraction plane of the QWP is turned by 45º from vertical.

E ~ 9 keV, crystal thickness = 500 µm.

Convergent beam, bent-Laue transmission case¶

One way to improve the low circular polarization rate in a divergent or convergent beam is to aperture it, with the obvious disadvantage of lowering the beam flux. Another, less obvious, way is to bend the crystal with the radius equal to the distance from the sample (focal point) to the QWP. In this particular example the Laue case is considered. The bending is done cylindrically in the diffraction plane which is at 45º to the initial polarization plane (horizontal). Watch the circular polarization rate at $\pm$ 32 arcsec departure angle.

E ~ 9 keV, crystal thickness = 500 µm.

Comparison of 1D-bent crystal analyzers¶

Files in /examples/withRaycing/06_AnalyzerBent1D

Rowland circle based analyzers¶

This study compares simply bent and ground-bent spectrometers utilizing Bragg and Laue analyzing crystals. The bending is cylindrical (one-dimensional).

Conditions:: Rowland circle diameter = 1 m, 70v × 200h µm² unpolarized fluorescence source, crystal size = 100meridional × 20saggittal mm².

Energy resolution was calculated as described in the CDR of a diced Johansson-like spectrometer at Alba/CLÆSS beamline. This calculation requires two images: 1) of a flat energy distribution source and 2) of a monochromatic source. The images are dispersive in the diffraction plane, which can be used in practice with a position sensitive detector or with a slit scan in front of a bulk detector. From these two images energy resolution δE was calculated and then 3) a verifying image was ray-traced for a source of 7 energy lines evenly spaced with the found step δE. Such images are shown for the four crystal geometries at a particular Bragg angle:

geometry	flat source	line source	7 lines
Bragg simply bent
Bragg ground bent
Laue simply bent
Laue ground bent

The energy distribution over the crystal surface is hyperbolic for Bragg and ellipsoidal for Laue crystals. Therefore, Laue crystals have limited acceptance in the sagittal direction whereas Bragg crystals have the hyperbola branches even for large sagittal sizes. Notice the full crystal coverage in the meridional direction for the two ground-bent cases.

geometry	flat source	line source	7 lines
Bragg simply bent
Bragg ground bent
Laue simply bent
Laue ground bent

As a matter of principles checking, let us consider how the initially unpolarized beam becomes partially polarized after being diffracted by the crystal analyzer. As expected, the beam is fully polarized at 45° Bragg angle (Brewster angle in x-ray regime). CAxis here is degree of polarization:

Bragg	Laue

Comments

The ground-bent crystals are more efficient as the whole their surface works for a single energy, as opposed to simply bent crystals which have different parts reflecting the rays of different energies.
When the crystal is close to the source (small θ for Bragg and large θ for Laue spectrometers), the images are distorted, even for the ground-bent crystals.
The Bragg case requires small pixel size in the meridional direction (~10 µm for 1-m-diameter Rowland circle) for a good spatial resolution but can profit from its compactness. The Laue case requires a big detector of a size comparable to that of the crystal but the pixel size is not required to be small.
The comparison of energy resolution in Bragg and Laue cases is not strictly correct here. While the former case can use the small beam size at the detector for utilizing energy dispersive property of the spectrometer, the latter one has a big image at the detector which is restricted by the size of the crystal. The size of the ‘white’ beam image is therefore correct only for the crystal size selected here. The Laue case can still be used in energy dispersive regime if 2D image analysis is utilized. At the present conditions, the energy resolution of Bragg crystals is better than that of Laue crystals except at small Bragg angles and low diffraction orders.
The energy resolution in ground-bent cases is not always better than that in simply bent cases because of strongly curved images. If the sagittal size of the crystal is smaller or sagittal bending is used, the advantage of ground-bent crystals is clearly visible not only in terms of efficiency but also in terms of energy resolution.

Johansson analyzer with elastically deformed crystal¶

Example in /examples/withRaycing/06_AnalyzerBent1D/01B_BentTT.py

A Johansson 1D ground-bent analyzer is studied here by ray tracing with two ways of reflectivity calculations: for a perfect crystal and for an elastically deformed crystal by means of the Takagi-Taupin equations. For higher reflection indices, the deviation from perfect crystal reflectivity becomes increasingly pronounced. Notice here an increased flux value and decreased energy resolution for the latter case.

study

flat source

line source

7 lines

perfect crystal

deformed crystal

(Takagi- Taupin)

Von Hamos analyzer¶

A von Hamos spectrometer has axial symmetry around the axis connecting the source and the detector. The analyzing crystal is cylindrically bent with the radius equal to the crystal-to-axis distance. In this scheme, the emission escape direction depends on Bragg angle (energy). In practice, the spectrometer axis is adapted such that the escape direction is appropriate for a given sample setup. If the emission escape direction is kept fixed relatively to the sample, the mechanical model becomes more complex and includes three translations and two rotations. In the figure below, the crystal is sagittally curved around the source-detector line. The detector plane is perpendicular to the sketch. Left: the classical setup [vH] with 2 translations. Right: the setup with an invariant escape direction.

The geometrical parameters for the von Hamos spectrometer were taken from [vH_SLS]: a diced 100 (sagittal) × 50 (meridional) mm² Si(444) crystal is curved with Rs = 250 mm. The width of segmented facets was taken equal to 5 mm (as in [vH_SLS]) and 1 mm together with a continuously bent case. The crystal thickness = 100 µm.

[vH] (1,2)

L. von Hámos, Röntgenspektroskopie und Abbildung mittels gekrümmter Kristallreflektoren II. Beschreibung eines fokussierenden Spektrographen mit punktgetreuer Spaltabbildung, Annalen der Physik 411 (1934) 252-260

[vH_SLS] (1,2)

J. Szlachetko, M. Nachtegaal, E. de Boni, M. Willimann, O. Safonova, J. Sa, G. Smolentsev, M. Szlachetko, J. A. van Bokhoven, J.-Cl. Dousse, J. Hoszowska, Y. Kayser, P. Jagodzinski, A. Bergamaschi, B. Schmitt, C. David, and A. Lücke, A von Hamos x-ray spectrometer based on a segmented-type diffraction crystal for single-shot x-ray emission spectroscopy and time-resolved resonant inelastic x-ray scattering studies, Rev. Sci. Instrum. 83 (2012) 103105.

The calculation of energy resolution requires two detector images: 1) of a flat energy distribution source and 2) of a monochromatic source. From these two images, energy resolution δE was calculated and then 3) a verifying image was ray-traced for a source of 7 energy lines evenly spaced with the found step δE. Such images are shown for different dicing sizes at a particular Bragg angle. In addition to perfect crystal reflectivity calculations, elastically deformed crystal reflectivity was calculated by means of the Takagi-Taupin equations (labelled TT).

crystal	flat source	line source	7 lines
diced 5 mm
diced 1 mm
not diced
not diced TT

With the coloring by stripe (crystal facet) number, the image below explains why energy resolution is worse when stripes are wider and the crystal is sagittally larger. The peripheral stripes contribute to aberrations which increase the detector image.

crystal	line source colored by stripe number
diced 5 mm
diced 1 mm

The efficiency of a von Hamos spectrometer is significantly lower as compared to Johann and Johansson crystals. The reason for the lower efficiency can be understood from the figure below, where the magnified footprint on the crystal is shown: only a narrow part of the crystal surface contributes to a given energy band. Here, in the 5-mm-stripe case a bandwidth of ~12 eV uses less than 1 mm of the crystal!

Comparison of Rowland circle based and von Hamos analyzers¶

An additional case was also included here: when a Johann crystal is rotated by 90⁰ around the sample-to-crystal line, it becomes a von Hamos crystal that has to be put at a correct distance corresponding to the 1 m sagittal radius. This case is labelled as “Johann as von Hamos”.

In comparing with a von Hamos spectrometer, one should realize its strongest advantage – inherent energy dispersive operation without a need for energy scan. This advantage is especially important for broad emission lines. Below, the comparison is made for two cases: (1) a narrow energy band (left figure), which is more interesting for valence band RIXS and which assumes a high resolution monochromator in the primary beam and (2) a wide energy band (right figure), which is more interesting for core state RIXS and normal fluorescence detection. The desired position on the charts is in the upper left corner. As seen in the figures, the efficiency of the von Hamos crystals (i) is independent of the energy band (equal for the left and right charts), which demonstrates truly energy-dispersive behavior of the crystals but (ii) is significantly lower as compared to the Johann and Johansson crystals. A way to increase efficiency is to place the crystal closer to the source, which obviously worsens energy resolution because of the increased angular source size. Inversely, if the crystal is put at a further distance, the energy resolution is improved (square symbols) but the efficiency is low because of a smaller solid angle collected. The left figure is with a narrow energy band equal to the 6-fold energy resolution. The right figure is with a wide energy band equal to 8·10 ^-4·E (approximate width of K β lines [Henke]).

Finally, among the compared 1D-bent spectrometers the Johansson type is the best in the combination of good energy resolution and high efficiency. It is the only one that can function both as a high resolution spectrometer and a fluorescence detector. One should bear in mind, however, two very strong advantages of von Hamos spectrometers: (1) they do not need alignment – a crystal and a detector positioned approximately will most probably immediately work and (2) the image is inherently energy dispersive with a flat (energy independent) detector response. The low efficiency and mediocre energy resolution are a price for the commissioning-free energy dispersive operation. Rowland circle based spectrometers will always require good alignment, and among them only the Johansson-type spectrometer can be made energy dispersive with a flat detector response.

Circular and elliptical von Hamos analyzers¶

The axial symmetry of the classical von Hamos spectrometer [vH] results in a close detector-to-sample position at large Bragg angles. A single detector may find enough space there but when the spectrometer has several branches, the corresponding detectors come close to each other, which restricts both the space around the sample and the accessible Bragg angle range. A solution to this problem could be an increased magnification of the crystal from the classical 1:1. The axis of the circular cilynder is then split into two axes representing the two foci of an ellipsoid, see the scheme below. The lower axis holds the source (sample) and the upper one holds the detector. The crystal in the figure has the magnification 1:1 for the circular crystal (left part) and 1.5:1 for the elliptical one (right part).

The crystal is diced along the cylinder axis with 1 mm pitch. The difference in the circular and elliptical figures is shown below.

The elliptical figure results in some aberrations, as seen by the monochromatic images below, which worsens energy resolution.

crystal	flat source	line source	7 lines
bent as circular cylinder
bent as elliptical cylinder

Comparison of 2D-bent Bragg crystal analyzers¶

Files in \examples\withRaycing\07_AnalyzerBent2D

This study compares simply bent (Johann) and ground-bent (Johansson) analyzers in diced and non-diced versions. This time the bending is two-dimensional with the sagittal radius equal to the distance from the crystal to the source-detector line. Such bending gives a family of Rowland circles all going through the source and the detector.

The conditions are equal to those in the previous section. The crystal size is 100 × 100 mm². The crystal facets of the diced version are 1.4meridional × 2.1sagittal mm² with 50 µm gaps. The non-diced crystals are 350 µm thick.

The images for flat and line sources were used to calculate energy resolution δE. After this, a 7-line source is created with the energy spacing between the lines equal to δE.

In addition to perfect crystal reflectivity calculations, elastically deformed crystal reflectivity was calculated by means of the Takagi-Taupin equations (labelled TT).

crystal	flat source	line source	7 lines
Johann
Johann TT
Johansson
Johansson TT
Johann diced
Johansson diced

Notice the energy distribution over the crystal area: the sagittal bending makes it uniform in the sagittal direction (here horizontal) and the ground-bent technology makes it uniform also in the meridional direction (here vertical):

crystal	footprint image	zoomed in footprint
Johann
Johann TT
Johansson
Johansson TT
Johann diced
Johansson diced

ALBA CLÆSS beamline¶

Files in \examples\withRaycing\08_CLAESS_BL

See the optical scheme of the beamline here.

This script produces images at various positions along the beamline.

The following 13 images are:

FSM image after the front end with the projected absorbed rays (red) at

the fixed front end mask,

upstream half and

downstream half of the movable front end mask

footprint on VCM,
footprint on the 1st crystal of DCM,
footprint on the 2nd crystal of DCM,
beam at the Bremsstrahlung block,
image at the foil holder of 4-diode XBPM,
footprint on VFM,
front collimator of the photon shutter,
image at the reducer flange 100CF-to-40CF,
image at the EH 4-blade slit,
image at the focal (sample) point.

The script also exemplifies the usage of touch_beam() for finding the optimal size of slits.

Gratings, FZPs, Bragg-Fresnel optics, cPGM beamline¶

Files in \examples\withRaycing\09_Gratings

Simple gratings¶

The following pictures exemplify simple gratings with the dispersion vector a) in the meridional plane and b) orthogonal to the meridional plane. Coloring is done by energy and by diffraction order.

Fresnel Zone Plate¶

This example shows focusing of a quasi-monochromatic collimated beam by a normal (orthogonal to the beam) FZP. The energy distribution is uniform within 400 ± 5 eV. The focal length is 2 mm. The phase shift in the zones is variable and is relative to the central ray. As expected, the phase shift does not influence the focusing properties and can be selected at will.

zoomed footprint on FZP	focal spot

Bragg-Fresnel optics¶

One can combine an arbitrarily curved crystal surface, also (variably) asymmetrically cut, with a grating or zone structure on top of it. The following example shows a Fresnel zone structure that focuses a collimated beam at q = 20 m, whereas the Bragg crystal provides good energy resolution. One can easily study how the band width affects the focusing properties (not shown here).

Generic cPGM beamline¶

This example shows a generic cPGM beamline aligned for a fixed focus regime. The angles at the mirrors equal 2 degrees, c_ff = 2.25, the line density is 1221 mm^-1.

An energy scan at a given vertical slit (here, 30 µm) between M3 and M4. Shown are images at the slit and at the final focus ‘Exp2’:

A vertical slit scan at a given energy (here, 40 eV) with a final dependency of energy resolution and flux on the slit size:

Multiple reflections¶

Files in \examples\withRaycing\10_MultipleReflect

Montel mirror¶

Montel mirror consists of two orthogonal mirrors positioned side-by-side. It is very similar to a KB mirror but more compact in the longitudinal direction. In a Montel mirror a part of the incoming beam is first reflected by one side of the mirror and then by the other side. Another part of the beam has the opposite reflection sequence. There are also single reflections on either side of the mirror. The non-sequential way of reflections makes such ray tracing impossible in most of ray-tracing programs.

The images below show the result of reflection by a Montel mirror consisting of a pair of parabolic mirrors. The mirrors can be selected by the user to have another shape. The coloring is by categories and the number of reflections. Notice a gap between the mirrors (here 0.2 mm) that transforms into a diagonal gap in the final image.

In the present example one can visualize the local footprints on either of the mirrors. The footprints are colored by the number of reflections.

Polycapillary¶

This example also demonstrates the technique of propagating the beam through an array of non-sequential optical elements. These are 397 glass capillaries, close packed into a hexagonal bunch. The capillaries here serve for collimating a divergent beam (fluorescence) born at the origin. Each capillary here is parabolically (user-supplied shape) bent such that the left end tangents are directed towards the source and the right ends are parallel to the y axis. The capillary radius here is constant but can also be given by a user-function. The images below show the geometry of the polycapillary: the screen at the entrance that is colored by the categories of the exit beam (the green rays are reflected by the capillaries, the orange ones are transmitted without any reflection and the red ones are absorbed) and a longitudinal cross-section (note very different scales of x and y axes).

Each capillary reflects the rays many times until they are out. The local footprints over the capillary surface are vs. the polar angle and the longitudinal coordinate. The coloring is histogrammed over incidence angle and number of reflections. The longitudinal coordinate s has its zero at the exit from the capillary.

	layer 1	layer 6	layer 12
n refl
θ

At the exit from the polycapillary the beam is expanded and attenuated at the periphery. The attenuation is due to the losses at each reflection; the number of reflections increases at the periphery, as shown by the colored histogram below. The phase space of the exit beam (shown is the horizontal one) shows the quality of collimation, see below. The divergence of ~1 mrad is large and cannot be efficiently used with a flat crystal analyzer.

One may attempt to add a second-stage to collimate the beam with ~1 mrad divergence. However, this would not work because the rays collimated by a polycapillary have a very large distribution of the ray origins over the longitudinal direction y, as shown below.

Powder Diffraction¶

Simulation of the real powder diffraction experiment on PETRA-III High Resolution Powder Diffraction Beamline P02.1. Uses Undulator source and double Laue Plate monochromator. Cerium Dioxide powder as the sample.

Warning

Heavy computational load. Requires OpenCL.

Gallery of plots and scripts 3. Wave propagation¶

These examples are based on wave diffraction via Kirchhoff integral.

Warning

You need a good graphics card for running these calculations!

Note

Consider the warnings and tips on using xrt with GPUs.

Diffraction from mirror surface¶

This example shows wave diffraction from a geometric source onto a flat or hemispheric screen. The source is rectangular with no divergence. The mirror has no material properties (reflectivity is 1) for simplicity. Notice the difference in the calculated flux between the rays and the waves.

rays	wave

The flux losses are not due to the integration errors, as was proven by variously dense meshes. The losses are solely caused by cutting the tails, as proven by a wider image shown below.

Slit diffraction with undulator source¶

This example shows wave diffraction of an undulator beam on a square slit. The calculations were done for (1) a monochromatic beam with zero electron beam emittance and (2) a finite energy band beam with the real electron emittance of Petra III ring. The resulted images are compared with experimentally measured ones [Zozulya_Sprung].

Notice that vertical stripes are less pronounced in the complete calculations because the horizontal emittance is much larger than the vertical one.

[Zozulya_Sprung]

A. Zozulya and M. Sprung, measured at P10 beamline, DESY Photon Science (2014) unpublished.

slit size (µm²)	ray tracing, zero emittance	diffraction, zero emittance	diffraction, real emittance	exp P10
100×100
150×150
200×200
300×300

Young’s experiment with undulator source¶

This example shows double slit diffraction of an undulator beam. The single slit width is 10 µm, the slit separation is variable (displayed is edge-to-edge distance), the slit position is 90 m from the source and the screen is at 110 m.

Diffraction from grating¶

Various gratings described in [Boots] have been tested with xrt for diffraction efficiency. The efficiency curves in [Boots] were calculated by means of the code peg which provides almost identical results to those by REFLEC [REFLEC] but with reportedly better convergence. In order to have comparison curves, we got the REFLEC results calculated by R. Sankari [Sankari] which were basically equal to those in [Boots].

[Boots] (1,2,3,4)

M. Boots, D. Muir and A. Moewes, Optimizing and characterizing grating efficiency for a soft X-ray emission spectrometer, J. Synchrotron Rad. 20 (2013) 272–285.

[REFLEC]

F. Schäfers and M. Krumrey, Technischer Bericht, BESSY TB 201 (1996).

[Sankari]

Sankari, private communication (2015).

The diffraction orders are shown below as transverse images and meridional (“polar exit angle”) and sagittal (“azimuthal exit angle”) cuts. Notice that the diffraction orders are positioned on the screen “by themselves”, i. e. with only the use of the Kirchhoff diffraction integral. Also before it was possible to work with grating diffraction orders in xrt within the geometrical ray tracing approach. In that approach the rays were deflected according to the grating equation. Here, in wave propagation, the grating equation was only used to position the screen.

Notice that in contrast to the conventional grating theories (also used in REFLEC), the diffraction orders here have the sagittal dimension. And that dimension has diffraction fringes and a variable width, too!

The resulting efficiency was obtained as ratio of the flux into the given order over the incoming flux. The incoming radiation was considered as uniform, parallel and fully coherent.

For the LEG (Low Energy Grating, see its properties in the figure below), the efficiency curves are pretty similar to those by REFLEC. The main difference is the low-energy part. Our 1st order does not decrease so rapidly. If we consider not the 2D exit angle but only the central azimuthal cut, the resulted low-energy efficiencies are very similar to that of REFLEC (not shown). Ref. [Boots] also provides experimental measurements which also have a rapid low-energy decrease. It seems that the detector had a pinhole that might cut the beam at low energies as the diffracted beam becomes wider there (see the transverse pictures above), which may explain lower measured efficiency.

For the IMP grating (“impurity”, see its properties in the figure below), the difference is bigger.

We believe that REFLEC is essentially wrong at high energies. If we mentally translate the working terraces of a blazed grating to form a continuous plane, we get a mirror at the pitch + blaze angle. By energy conservation, the overall grating efficiency (the sum into all orders) cannot be higher than the reflectivity of such a mirror. The REFLEC curves can violate this limit even for a single order, compare with the blue curve in the figure above. The reason for such behavior seems to be the artificially shadowless illumination by the incoming wave. Indeed, REFLEC assumes the complete saw profile to work in the diffraction, whereas the back side and a portion of the front side behind it stay in the shadow. We compare the two gratings shown below, one is with 90 degree anti-blaze angle and the other is with pitch as anti-blaze angle.

REFLEC gives different efficiencies for these two cases (see above) whereas xrt cannot distinguish them. We tried to artificially remove the shadows by making the surface “emit” a coherent wave. The result was an increase in the high-energy efficiency, similarly to the REFLEC’s behaviour.

The factors which definitely will affect the efficiency are (1) restricted coherence radius and (2) roughness. Both will be added into this example in a later release of xrt.

We are open for further discussion on the above results with interested scientists.

Diffraction from FZP¶

This examples demonstrates diffraction from a Fresnel Zone Plate with variously thick outer zone and at variable energy. The radial intensity distribution is shown in the figure below for a 70-nm-outer-zone FZP. Notice that the 2nd order was also calculated and together with other even orders indeed results in vanishing intensity.

The energy dependence of efficiency for 3 different FZPs is shown below. The horizontal bars mark the expected $1/m^2\pi^2$ levels for odd orders and 25% transmission for the 0th order. Watch how a zone plate becomes a band pass filter as the outer zone size approaches the wavelength, here ~10 nm.

SoftiMAX at MAX IV¶

The images below are produced by scripts in \examples\withRaycing\14_SoftiMAX.

The beamline will have two branches: - STXM (Scanning Transmission X-ray Microscopy) and - CXI (Coherent X-ray Imaging),

see the scheme provided by K. Thånell.

STXM branch¶

Rays vs. hybrid

The propagation through the first optical elements – from undulator to front end (FE) slit, to M1, to M2 and to plane grating (PG) – is done with rays:

FE	M1	M2	PG

Starting from PG – to M3, to exit slit, to Fresnel zone plate (FZP) and to variously positioned sample screen – the propagation is done by rays or waves, as compared below. Despite the M3 footprint looks not perfect (not black at periphery), the field at normal surfaces (exit slit, FZP (not shown) and sample screen) is of perfect quality. At the best focus, rays and waves result in a similar image. Notice a micron-sized depth of focus.

	rays	wave
M3
exit slit
sample

Influence of emittance

Non-zero emittance radiation is treated in xrt by incoherent addition of single electron intensities. The single electron (filament) fields are considered as fully coherent and are resulted from filament trajectories (one per repeat) that attain positional and angular shifts within the given emittance distribution. The following images are calculated for the exit slit and the focus screen for zero and non-zero emittance (for MAX IV 3 GeV ring: ε_x=263 pm·rad, β_x=9 m, ε_z=8 pm·rad, β_z=2 m). At the real emittance, the horizontal focal size increases by ~75%. A finite energy band, as determined by vertical size of the exit slit, results in somewhat bigger broadening due to a chromatic dependence of the focal length.

	0 emittance	real emittance	real emittance, finite energy band
exit slit
sample

Correction of emittance effects

The increased focal size can be amended by closing the exit slit. With flux loss of about 2/3, the focal size is almost restored.

	80 µm exit slit	20 µm exit slit
exit slit
sample

Coherence signatures

The beam improvement can also be viewed via the coherence properties by the four available methods (see Coherence signatures). As the horizontal exit slit becomes smaller, one can observe the increase of the coherent fraction ζ and the increase of the primary (coherent) mode weight. The width of degree of coherence (DoC) relative to the width of the intensity distribution determines the coherent beam fraction. Both widths vary with varying screen position around the focal point such that their ratio is not invariant, so that the coherent fraction also varies, which is counter-intuitive. An important advantage of the eigen-mode or PCA methods is a simple definition of the coherent fraction as the eigenvalue of the zeroth mode (component); this eigenvalue appears to be invariant around the focal point, see below. Note that the methods 2 and 3 give equal results. The method 4 that gives the degree of transverse coherence (DoTC) is also invariant around the focal point, see DoTC values on the pictures of Principal Components.

	80 µm exit slit	20 µm exit slit
method 1
method 2
method 3, method 4b

CXI branch¶

2D vs 1D

Although the sample screen images are of good quality (the dark field is almost black), the mirror footprints may be noisy and not well convergent in the periphery. Compare the M3 footprint with that in the previous section (STXM branch) where the difference is in the mirror area and thus in the sample density. The used 10⁶ wave samples (i.e. 10¹² possible paths) are not enough for the slightly enlarged area in the present example. The propagation is therefore performed in separated horizontal and vertical directions, which dramatically improves the quality of the footprints. Disadvantages of the cuts are losses in visual representation and incorrect evaluation of the flux.

	2D	1D horizontal cut	1D vertical cut
M3 footprint
sample screen

Flat screen vs normal-to-k screen (wave front)

The following images demonstrate the correctness of the directional Kirchhoff-like integral (see Sequential propagation). Five diffraction integrals are calculated on flat screens around the focus position: for two polarizations and for three directional components. The latter ones define the wave fronts at every flat screen position; these wave fronts are further used as new curved screens. The calculated diffraction fields on these curved screens have narrow phase distributions, as shown by the color histograms, which is indeed expected for a wave front by its definition. In contrast, the flat screens at the same positions have rapid phase variation over several Fresnel zones.

Note

In the process of wave propagation, wave fronts – surfaces of constant phase – are not used in any way. We therefore call it “wave propagation”, not “wave front propagation” as frequently called by others. The wave fronts in this example were calculated to solely demonstrate the correctness of the local propagation directions after having calculated the diffracted field.

flat screen	curved screen (wave front)

The curvature of the calculated wave fronts varies across the focus position. The wave fronts become more flat as one approaches the focus, see the figure below. This is in contrast to ray propagation, where the angular ray distribution is invariant at any position between two optical elements.

Rays, waves and hybrid

The following images are horizontal cuts at the footprints and sample screens calculated by

rays,
rays + waves hybrid (rays up to PG and wave from PG) and
purely by waves.

	rays	hybrid	waves
front end slit		same as rays
footprint on M1		same as rays
footprint on M2		same as rays
footprint on PG		same as rays
footprint on M3
exit slit
footprint on M4
footprint on M5
sample screen

Coherence signatures

This section demonstrates the methods 1 and 3 from Coherence signatures. Notice again the difficulty in determining the width of DoC owing to its complex shape (at real emittance) or the restricted field of view (the 0 emittance case). In contrast, the eigen mode analysis yields an almost invariant well defined coherent fraction.

	0 emittance	real emittance
method 1
method 3

Defocusing by a distorted mirror¶

The images below are produced by \examples\withRaycing\13_Warping\warp.py.

This example has two objectives:

to demonstrate how one can add a functional or measured figure error to an ideal optical element and
to study the influence of various figure errors onto image non-uniformity in focused and defocused cases. The study will be done in ray tracing and wave propagation, the latter being calculated in partial coherence with the actual emittance of the MAX IV 3 GeV ring.

Here, a toroidal mirror focuses an undulator source in 1:1 magnification. The sagittal radius of the torus was determined for p = q = 25 m and pitch = 4 mrad. Defocusing in horizontal is done by going to a smaller pitch angle, here 2.2 mrad, and in vertical by unbending the meridional figure.

Three distorted surfaces are of Gaussian, waviness and as measured shapes, see below. They were normalized such that the meridional slope error be 1 µrad rms. The surfaces are determined on a 2D mesh. Interpolation splines for the height and the normal vector are found at the time of mirror instantiation and used in two special methods: local_z_distorted and local_n_distorted, see Section Distorted surfaces. If the distorted shape is known analytically, as for waviness, the two methods may directly invoke the corresponding functions without interpolation. The scattered circles in the figures are random samples where the height is calculated by interpolation (cf. the color (height) of the circles with the color of the surface) together with the interpolated normals (white arrows as projected onto the xy plane).

Gaussian	waviness	mock NOM measurement

Defocused images reveal horizontal stripes seen both by ray tracing and wave propagation. Notice that wave propagation ‘sees’ less distortion in the best focusing case.

	ray tracing	wave propagation
ideal
Gaussian
waviness
mock NOM

Logo¶

Uses dummy backend. File: \examples\withDummy\logo_xrt.py

xrt logo (on the right, in two versions) created from a “flat” python logo (on the left).

————–>

Speed tests¶

The scripts used for these tests can be found in \tests\speed\

The following computers have been used:

Note

The tests here were reduced in the number of rays/samples as compared to real calculations to let them run reasonably quickly. Longer calculations would demonstrate yet bigger difference between the slowest and the fastest cases, as the overheads (job distribution, collecting of histograms and plotting) would become relatively less important.

The tables below show execution times in seconds. Some cells have two values: for Python 2 and for Python 3.

Multithreading and multiprocessing in ray tracing¶

Script: \tests\speed\1_SourceZCrystalThetaAlpha_speed.py.

The test is based on the example \examples\withRaycing\07_AnalyzerBent2D\01BD_SourceZCrystalThetaAlpha.py.

This script calculates diffraction of a 2D-bent diced crystal anayzer from three types of geometric source.

system		1	multithreading		multiprocessing
system		1	2	4	2	4
[1]	Windows 10 64 bit	510.2, 539.2	306.0, 318.5	220.1, 220.6	458.8, 573.4	338.8, 441.2
[1]	Ubuntu 16.04 64 bit	436.1, 446.7	257.2, 263.3	183.3, 188.6	246.5, 250.0	157.1, 158.7
[2]	Windows 10 64 bit	359.9, 546.9	290.8, 292.0	218.7, 219.7	362.2, 436.4	290.8, 363.6
[2]	Ubuntu 16.10 64 bit	293.3, 293.1	172.3, 173.6	169.5, 168.3	173.3, 172.4	156.1, 154.6

OpenCL performance with Undulator source¶

Script: \tests\speed\2_synchrotronSources_speed.py.

The test is based on the example \examples\withRaycing\01_SynchrotronSources\synchrotronSources.py.

This script calculates characteristics of an undulator source at energies around one harmonic.

system		no OpenCL	OpenCL on CPU	OpenCL on GPU
[1]	Windows 10 64 bit	1471 1385	36.0 34.1	25.7 23.9
[1]	Ubuntu 16.04 64 bit	950 950	34.6 35.4	20.6 21.0
[2]	Windows 10 64 bit	1801 1909	61.0 60.3	126 123
[2]	Ubuntu 16.10 64 bit	1245 1255	57.6 60.2	122 127
[9]	local			30.0
[9]	ZMQ 1Gb			182.9

OpenCL performance with wave propagation¶

Script: \tests\speed\3_Softi_CXIw2D_speed.py.

The test is based on the example \examples\withRaycing\14_SoftiMAX\Softi_CXIw2D.py.

This script calculates several consecutive wave propagation integrals from the source down to the focus at the sample. Here, each wave is represented by 2·10⁵ samples and thus each integral considers 4·10¹⁰ scattering paths. Such calculations are impossible in numpy and have to be carried out with the help of OpenCL. Even with OpenCL, these calculations are not feasible on low end graphics cards and therefore are not exemplified here on a laptop. Notice also that for the real example a larger number of samples, > 10⁶ (i.e. > 10¹² scattering paths), should be opted for better convergence.

system		OpenCL on CPU	OpenCL on GPU
[1]	Windows 10 64 bit	637 635	76.8 76.4
[1]	Ubuntu 16.04 64 bit	602 605	71.1 71.0
[3]	1×K80		196
[3]	4×K80		76.5
[4]	1×P100		53.0
[4]	4×P100		25.8
[5]	Xeon E5-2650 v3	321
[6]	Xeon E5-2650 v4	251
[7]	Xeon Gold 6130	162
[8]	1×A100		17.5
[8]	2×A100		11.5
[9]	local		40.4
[9]	ZMQ 1Gb		48.4

Summary¶

Except for the case of computing with OpenCL on GPU, calculations in Linux are usually significantly faster than in Windows. Especially when using multithreading or multiprocessing, the execution in Linux is dramatically faster.
There is no significant difference in speed between Python 2 and Python 3, except for multiprocessing in Windows, where Python 2 performs better.
For geometric ray tracing a decent laptop can be a reasonable choice.

The MIT License¶

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

normalized to local maximum
global normalization, type 1
global normalization, type 2

xrt (XRayTracer)¶

GUIs¶

xrtQook – a GUI for creating a beamline¶

xrtGlow – an interactive 3D beamline viewer¶

xrtBentXtal – a GUI for bent crystal calculations¶

Features of xrt¶

Dependencies¶

Get xrt¶

Installation¶

Get help¶

Citing xrt¶

Acknowledgments¶

User’s Guide¶

Using xrt¶

Scripting in python¶

Running a python script¶

Interacting with the plots¶

Saving the results¶

Using xrtQook for script generation¶

Using custom optical elements in xrtQook¶

Notes on using xrtGlow¶

Using xrt as x-ray calculator¶

Raycing backend¶

Coordinate systems¶

Units¶

Beam categories¶

Scripting in python¶

Sources¶

Sampling strategies¶

Geometric sources¶

Synchrotron sources¶

Comparison of synchrotron source codes¶

Using xrt synchrotron sources on a mesh¶

Undulator spectrum across a harmonic¶

Undulator spectrum at very high energies¶

Single-electron and multi-electron undulator radiation¶

Undulator spectrum with tapering¶

Undulator spectrum in transverse plane¶

Undulator radiation in near field¶

Field phase in near field¶

Undulator source size dependent on energy spread and detuning¶

Optical elements¶

Distorted surfaces¶

Tests of Optical elements¶

Materials¶

Predefined Materials¶

Tests of Materials¶

Apertures¶

Screens¶

Wave propagation (diffraction)¶

Time dependent diffraction¶

Stationary diffraction¶

Normalization¶

Sequential propagation¶

Usage¶

Coherence signatures¶

Coherence analysis and related plotting¶

Typical logic for a wave propagation study¶

Coherent mode decomposition and propagation¶

A model case of 5 macro-electrons¶

A model case of 5000 macro-electrons¶

Tests of wave propagation¶

Shadow backend¶

Description of shadow¶

Preparation for a shadow run¶

Scripting in python¶

Modifying input files¶

Writing a loop generator¶

Customizing your plots¶

Customizing the look of the graphs¶

Saving the histogram arrays¶

Controlling ray tracing execution¶

Using multiprocessing and multithreading¶

Calculations on GPU¶

OpenCL Server¶

Script Components¶

Running the Script¶

Gallery of plots and scripts 1. Synchrotron sources¶

Synchrotron sources¶

Bending magnet¶