SimOptDevice : a library for virtual optical experiments

Virtual experiments have become an indispensable tool for the design and the accuracy assessment of novel measurement procedures and instruments. Virtual experiments are particularly relevant in modern optics due to its challenging demands for highly accurate measurements. This paper introduces SimOptDevice, a flexible library for opto-mechanical virtual experiments. After describing the scope and general structure of the library, its underlying mathematical tools used for solving the related numerical tasks are described. Finally, the application of SimOptDevice to a recent interferometric measurement procedure is presented. Copyright statement. The author’s copyright for this publication is transferred to the Physikalisch-Technische Bundesanstalt (PTB).


Introduction
Following the advance of technology and the demand for highly accurate measurements, optical instruments and experiments have become very complex in recent years. In addition, sophisticated data analysis has become an important part of modern optical measurement devices. To ensure that a measurement principle is fit for its purpose, it is beneficial to first test it in a virtual environment prior to building the physical setup. That way, experimenters can save time and costs in the development of novel procedures. Furthermore, virtual experiments are often essential for the assessment of accuracies that can be reached.
The Physikalisch-Technische Bundesanstalt (PTB) has developed the SimOptDevice software library for optical virtual experiments. SimOptDevice is a flexible library implemented in MATLAB ® (MATLAB, 2018) that covers a large range of applications. It facilitates the design of experiments and allows one to develop and test measurement procedures. Furthermore, SimOptDevice is used to optimize existing measurement procedures with respect to measurement time and uncertainty. The latter is particularly relevant for PTB, which aims at measurements at the highest level of accuracy. Successful examples using SimOptDevice include accuracy evaluation of interferometric measurements of a synchrotron mirror (Wiegmann et al., 2011), development of a deflectometric flatness reference at PTB (Schulz et al., 2010;Ehret et al., 2012), and accuracy tests for multispectral imaging systems (Dierl et al., 2018).
In this paper we will describe the structure of SimOptDevice and its underlying mathematical methods. We will address the key issues of optical virtual experiments and refer to the corresponding solutions implemented in SimOptDevice. The purpose of this paper is to describe the mathematical methods necessary for implementing a complex optical simulation environment, enabling readers to set up a software solution of their own.

R. Schachtschneider et al.: SimOptDevice
The general structure of SimOptDevice and its basic principles are described in Sect. 2. In Sect. 3 we introduce the mathematical methods used for the solution of related numerical problems. A prominent example is then presented in Sect. 4. Finally, some conclusions are given.

Scope and general description
The SimOptDevice library can be used to perform virtual optical measurements with complex beam paths running through a series of optical elements. It considers nested scanning stages accounting for translations and rotations, and supports the use of various sensors such as cameras. SimOpt-Device is based on the application of ray optics.
The basic principle of SimOptDevice is a system of hierarchical coordinate systems, combined with ray tracing routines. Within each coordinate system, optical elements can be placed. A related local coordinate system is assigned to each considered optical element, along with a superordinate coordinate system relative to which the local coordinate system is defined. In this way, a tree structure of coordinate systems is built. Each local coordinate system can undergo individual rotations and translations with respect to its superordinate system. The coordinates of each element can be transformed into any of the other employed coordinate systems. Those transformations are made simple by using homogeneous coordinates which are introduced in Sect. 3.1.
The power of SimOptDevice lies in tracing rays and performing ray aiming accurately and efficiently according to the laws of refraction and reflection while being in control of all optimization parameters and the applied algorithms. Using the library for our experiments, we view and verify all intermediate results, which is very helpful for tuning the algorithms for each specific problem. Whereas ray tracing follows a ray through the optical system when start point and direction are given, ray aiming seeks the path through the system for given start and end points. The latter is a highly nonlinear optimization problem. Details are described in Sect. 3.
The accurate modelling of all elements of a measurement setup is another advantage of SimOptDevice. This includes not only optical elements like lenses, mirrors and sensors but also linear stages and rotary tables. Ensembles of elements can be saved and reused in other virtual experiments.
During development, the software has been successfully checked against ray tracing results obtained by ZEMAX. This included comparisons of optical path lengths and of points reached by ray tracing. The differences were in the sub-nanometre range.

Coordinate transformations
One key feature of SimOptDevice is its easy way to transform coordinates from any local system to any other coordi-nate system within the defined structure (see Fig. 1). Coordinates are defined as homogeneous coordinates (see e.g. Cox et al., 2007, p. 357 ff.). Rotations and translations can be represented as matrices, and it is simple to connect transformations and compute a transformation's inverse. Homogeneous coordinates for describing points in R 3 have four components. A vector v = (x, y, z) T is represented in homogeneous coordinates as v = (x , y z , w) T , where w = 0 is arbitrary and x = x/w, y = y/w, and z = z/w. Usually, w is set to one.
With this coordinate definition, a translation or rotation is represented by a 4 × 4 matrix: Translation : Rotations about the other axes are defined analogously. The composition of several homogeneous transformations H 1 , H 2 , . . . , H N is equal to a single homogeneous transformation Using Eq. (1), the translation and rotation of a coordinate system with respect to its superordinate system are represented by a single matrix. The inverse of a homogeneous transformation is represented by the inverse of the corresponding transformation matrix. The aforementioned transformation from one local source system S a to another destination system S b is performed with the help of a common superordinate system S c : For a more comprehensive introduction to homogeneous coordinates, see e.g. Cox et al. (2007) or Hartley and Zisserman (2003).

Ray tracing
Ray propagation through the defined object is performed by sequential ray tracing. For this tracing method the order of elements passed by a ray is known in advance and optical paths are computed element by element. SimOptDevice computes the ray paths locally in each optical element system. At the boundary T i between two elements (e.g. a surface of a lens), the current intersection point p i of the ray with that boundary along with the new ray direction e i are calculated in the new local system (cf. Sect. 3.1). The ray is then traced until it intersects with the next boundary, and so on. The intersection point with T i and the new direction are a function of the previous intersection point and ray direction: The corresponding instrument is illustrated. The position and orientation of each of the instrument's subsystems are defined with respect to the superordinate system through a transformation of homogeneous coordinates. A transformation matrix H can be computed for transformations from one system to another (see Eq. 2). It is a function of the source system and the destination system (in this example Sensor and Topo, respectively). The common superordinate system Ts_Frame is needed for the computation of the transformation matrix (cf. Eq. 2). Subsequently, coordinates of a point in one system can be transformed to any other system by multiplication by the corresponding transformation matrix H.
Function f i entails performing the following steps (cf. Fig. 2 for an illustration): 1. transformation of p i−1 and e i−1 from system T i−1 to T i , 2. calculation of the ray's geometrical path length l i between T i−1 and T i , 3. determination of next intersection point p i = p i−1 + l i · e i−1 , 4. calculation of normal vector n i on T i in p i , and 5. calculation of the new ray direction e i according to the laws of refraction and reflection. This is classical ray tracing with the particular feature that at each boundary the intersection point and the new direction are calculated in the new local coordinate system. Step 3 is calculated analytically for planes and spherical surfaces and has to be calculated numerically for more complex surfaces like Zernike surfaces, aspheres, or surfaces described by a Gauss function. If the ray passes N surfaces, there are N different functions f i . The last intersection point and ray direction can be expressed as a function of the start point and direction by concatenation of functions f 1 to f N : Furthermore, for each crossing of a boundary T i a Jacobian matrix containing the partial derivatives of the intersection point p i and the ray direction e i with respect to the previous Figure 2. Schematics of a ray tracing step. Each topography T has its own coordinate system. The intersection point p i with topography T i is computed from the previous intersection point p i−1 and the previous ray direction e i−1 . Subsequently, the new ray direction e i at T i is computed according to the laws of reflection and refraction.
intersection point p i−1 and the previous ray direction e i−1 is calculated. It is associated with Eq. (3) and has the following form: where x i and y i are the x and y components of the point p i , andx i andŷ i are the x and y components of the direction e i . We only need the derivatives with respect to the x and y components, since the intersected topographies are known and the z component is fixed at the intersection point.
Since the direction vector has unit length, its z component can be calculated byẑ = 1 −x 2 −ŷ 2 . The Jacobian matrix J corresponding to the concatenation of N − 1 interfaces i = 1 . . . N − 1 associated with Eq. (4) is obtained by multiplication of the corresponding matrices: The inverse Jacobian matrix J −1 describes partial derivatives of the ray's behaviour along the same path in the opposite direction. These Jacobian matrices are used extensively when performing the ray aiming. Their usage is explained in more detail in Sect. 3.3. A detailed description of the ray tracing and ray aiming procedures can also be found in Fortmeier (2016, in German).

Ray aiming
Given a start point p 0 and an end point p dest N the task of ray aiming is to find the start direction e 0 from p 0 in order to minimize the distance between p N and the desired end point p dest N (see Fig. 3): where || · || denotes the Euclidean norm and p N is reached through application of ray tracing for a chosen start direction e 0 . For successful ray aiming the resulting norm in Eq. (6) is close to zero. This is usually a highly nonlinear problem. It can be solved using one of MATLAB ® 's parallel nonlinear solver routines, e.g. lsqnonlin. The latter is a solver for nonlinear least-squares problems in which the trust-region-reflective or Levenberg-Marquard algorithm can be applied. For better convergence of the solver, the ray tracing Jacobian matrix from Eq. (5) is utilized in this step. The Jacobian matrices calculated during ray tracing are used to update the start direction e 0 when minimizing the distance between point p N and point p dest N . Furthermore, SimOptDevice delivers the Jacobian matrices for the change in total optical path lengths with respect to changes in a point p i or a ray direction e i for each surface T i (Fortmeier et al., 2014). It is calculated analytically, thereby omitting the computationally expensive extra ray tracing steps for the numerical differentiation. The Jacobian matrices for ray aiming are important in many interferometric applications, e.g. the tiltedwave interferometer (cf. Sect. 4), where the change in the optical path length of a ray is a significant quantity of the experiment.
A requirement for successful ray aiming is that the destination point p dest N can be reached. Therefore, the valid area on the CCD (or any other final surface) is determined in a preceding step. A ray aiming between the light source and some characteristic points at the smallest aperture in the optical system is performed. The final ray directions at the aperture are used to trace the rays to the CCD, thereby defining Figure 3. Ray aiming principle: the start and end points p 0 and p dest N of a ray are given. The task is to find the correct start angle e 0 such that the desired end point is hit. This is achieved by minimizing the distance between points p N and p dest N , when p N is reached through application of ray tracing for a chosen start direction, e 0 . the area that can be reached by rays from the sources. A ray aiming is regarded as successful if the norm of the difference in Eq. (6) is smaller than 1 nm. Local minima that prevent the algorithm from converging are detected and such rays are masked. Using the Jacobian matrix, the ray aiming usually converges in about four steps. Without the Jacobian matrix this takes approximately 12 to 15 steps, depending on the particular situation.

Application example: tilted-wave interferometer
The tilted-wave interferometer (TWI) is an interferometric measurement device for form measurements of asphere (Garbusi et al., 2008) and freeform surfaces. Figure 4 illustrates the simulation of interferograms for a measurement of an asphere. For SimOptDevice, the TWI is a particular application example since virtual experiments are not only used for designing the instrument and measurement setup, but are also an integrated part of data analysis. More precisely, we use SimOptDevice to obtain a numerical model of the experiment in dependence on the unknown form of the asphere or freeform surface under test. The latter is retrieved as the solution of a nonlinear inverse problem where in our implementation SimOptDevice constitutes the model. The inverse problem is solved iteratively, which is facilitated through the parallel computing capabilities of SimOptDevice. The nonlinear problem is linearized with the help of Jacobian matrices (Fortmeier et al., 2014). The solution of the inverse problem is found with the MATLAB ® routine fsolve. It is possible to use different solver algorithms with fsolve. We typically choose trust-region-dogleg or trust-region-reflective. Trust region algorithms locally approximate the cost function by a quadratic function and progress towards the minimum of this approximation.
In order to accomplish an accurate recovery of the specimen surface, it is necessary to obtain an accurate model of the TWI (Garbusi et al., 2008). Differences between the numerical model of the instrument and the actual instrument have to be characterized very accurately prior to a measurement. For this calibration process, well-known test specimens (typically spheres) are measured with the TWI (Baer et al., 2014). The numerical model is then updated by adjusting the Zernike polynomial parametrization of two reference surfaces in order to account for the remaining discrepancies.
SimOptDevice has also been used for uncertainty evaluation (Fortmeier et al., 2017) of the TWI and for the exploration of new measurement concepts with the TWI .

Conclusions
SimOptDevice is a versatile library for conducting virtual opto-mechanical experiments that has been applied successfully in several projects and studies. SimOptDevice can model a large number of optical elements and sensors which can be combined flexibly to cover a wide range of experimental setups.
We have explained the mathematical concepts within the library. A detailed description of our ray tracing and ray aiming procedures and of the determination of Jacobian matrices, needed for efficiently solving the nonlinear inverse problems, was given. This will be useful for readers interested in implementing virtual optical experiments. We also conclude that SimOptDevice can be used to simulate very complex opto-mechanical systems.
Data availability. No data sets were used in this article.
Competing interests. The authors declare that they have no conflict of interest.

Special issue statement.
This article is part of the special issue "Sensors and Measurement Systems 2018". It is a result of the "Sensoren und Messsysteme 2018, 19. ITG-/GMA-Fachtagung", Nürnberg, Germany, from 26 June 2018 to 27 June 2018.