5.2. Tracing Procedure

5.2.1. Tracing Process

The tracing process primarily involves three key steps: detecting surface hits, calculating refraction indices and directions, and managing rays that do not make contact with a surface. For ideal lenses, filters, and apertures, only one instance of surface hit detection is required, whereas a standard lens requires two.

Utilizing multiple threads for the majority of tracing tasks significantly enhances computational speed. Each thread is assigned a specific set of rays, ensuring that all processing, including ray creation and system-wide propagation, occurs exclusively within that thread.

Although not explicitly detailed here, the refraction calculation function also computes new ray polarization and weights.

../_images/TracePAP.svg

Fig. 5.6 Tracing process flowchart.

5.2.2. Refraction

The law of refraction is typically expressed in terms of input and output angles. However, for the tracing process, a vector-based formulation is more convenient as it avoids the need to calculate any angles directly.

The following figure illustrates the refraction of a ray on a curved surface.

../_images/refraction_interface.svg

Fig. 5.7 Refraction on a curved interface.

The refractive indices of the media are denoted as \(n_1\) and \(n_2\), while \(s\) and \(s'\) represent the input and output propagation vectors, respectively. It is essential to normalize these vectors, along with the normal vector \(n\), to ensure accurate calculations. Additionally, \(s\) and \(n\) must point in the same half-space direction, which is mathematically ensured by the condition \(s \cdot n \geq 0\).

An equation for this vector-based formulation of the refraction law can be found in [1] and [2]:

(5.23)\[s^{\prime}=\frac{n_1}{n_2} s-n\left\{\frac{n_1}{n_2}(n s) -\sqrt{1-\left(\frac{n_1}{n_2}\right)^{2}\left[1-(n s)^{2}\right]}\right\}\]

In the event of total internal reflection (TIR), the root argument becomes negative. In optrace, TIR rays are absorbed since reflections are not modeled within the system. When this occurs, the user receives a message to the standard output.

5.2.3. Polarization

The following calculations are similar to those described in [3].

The polarization vector \(E\) can be decomposed into two components: \(E_\text{p}\), which lies in the plane formed by the surface normal and the incidence vector, and \(E_\text{s}\), which is perpendicular to this plane. When refraction occurs at an interface, the \(E_\text{s}\) component remains unchanged for both the incident and refracted ray vectors, \(s\) and \(s'\). However, the \(E_\text{p}\) component is rotated around \(E_\text{s}\) towards \(s'\), resulting in a new component \(E\text{p}'\).

It’s important to note that, in our calculations, all vectors are unity vectors. The magnitudes of the polarization components are encoded in the scaling factors \(A_\text{tp}\) and \(A_\text{ts}\).

../_images/refraction_interface_polarization.svg

Fig. 5.8 Ray polarization components before and after refraction.

Case 1:

For \(s \parallel s'\) the new polarization vector is equal to the initial polarization.

Case 2

For \(s \nparallel s'\) the new polarization vector differs. In optics, the polarization vector and its components are orthogonal to the propagation direction. Moreover, the two polarization components are perpendicular to each other. Given that all vectors mentioned are unity vectors, we can perform the following calculations:

(5.24)\[\begin{split}\begin{align} E_\text{s} &= \frac{s' \times s}{|| s' \times s ||}\\ E_\text{p} &= E_\text{s} \times s\\ E_\text{p}' &= E_\text{s} \times s'\\ \end{align}\end{split}\]

Since \(||E_\text{p}|| = ||E_\text{s}|| = ||E|| = 1\), the amplitude components are then:

(5.25)\[\begin{split}\begin{align} A_\text{tp} &= E_\text{p} \cdot E\\ A_\text{ts} &= E_\text{s} \cdot E\\ \end{align}\end{split}\]

For the new polarization unity vector, which also composed of two components, we finally get:

(5.26)\[E' = A_\text{ps} E_\text{s} + A_\text{tp} E_\text{p}'\]

5.2.4. Transmission

The new ray powers are determined by calculating the transmission, which can be derived from the properties already established during the refraction calculations. According to the Fresnel equations, the transmission of light depends on its polarization direction. The following equations illustrate this behavior, as detailed in [4].

(5.27)\[t_{\mathrm{s}}=\frac{2\, n_{1} \cos \varepsilon}{n_{1} \cos \varepsilon+n_{2} \cos \varepsilon'}\]
(5.28)\[t_{\mathrm{p}}=\frac{2\, n_{1} \cos \varepsilon}{n_{2} \cos \varepsilon+n_{1} \cos \varepsilon'}\]
(5.29)\[T=\frac{n_{2} \cos \varepsilon'}{n_{1} \cos \varepsilon} \left( (A_\text{ts} t_\text{s})^2 + (A_\text{tp} t_\text{p})^2 \right)\]

The polarization components \(A_\text{ts}\) and \(A_\text{tp}\) are obtained from equations (5.25). The cosine terms involved are calculated using the direction and normal vectors as follows: \(\cos \varepsilon = n \cdot s\) and \(\cos \varepsilon' = n \cdot s'\).

For light that hits the surface perpendicularly, this results in an expression that is independent of polarization, as noted in [5]:

(5.30)\[T_{\varepsilon=0} = \frac{4 n_1 n_2 }{(n_1 + n_2)^2}\]

5.2.5. Refraction at an Ideal Lens

A ray with an unnormalized direction vector \(s_0\) intersects the lens at the point \(P = (x_0, y_0, 0)\) on a lens with a focal length \(f\). The corresponding point on the focal plane is \(P_f = (x_f, y_f, f)\).

According to optics, ideal parallel rays converge at the same position on the focal plane. Therefore, a ray with the same direction, but intersecting the lens at the optical axis, can be used to determine the position \(P_f\). This is illustrated in Fig. 5.9.

../_images/ideal_refraction.svg

Fig. 5.9 Geometry for refraction on an ideal lens.

Cartesian Representation

Calculating positions \(x_f\) and \(y_f\) is straightforward by evaluating the linear ray equations \(x(z)\) and \(y(z)\) at \(z=f\).

For \(x_f\), we have:

(5.31)\[x_f = \frac{s_{0x}}{s_{0z}} f\]

Similarly, for \(y_f\):

(5.32)\[y_f = \frac{s_{0y}}{s_{0z}} f\]

It is important to note that \(s_{0z} = 0\) is prohibited, as all rays are required to have a positive z-direction component.

With the point \(P_f\) known, the outgoing propagation vector \(s_0'\) can be calculated as follows:

(5.33)\[\begin{split}s_0' = P_f - P = \begin{pmatrix} \frac{s_{0x}}{s_{0z}}f - x_0 \\ \frac{s_{0y}}{s_{0z}}f - y_0 \\ f \end{pmatrix}\end{split}\]

Normalizing this vector gives us:

(5.34)\[s' = \frac{s_0'}{||s_0'||}\]

Angular Representation

Consider the x-component of the propagation vector:

(5.35)\[s_{0x}' = \frac{s_{0x}}{s_{0z}}f - x_0\]

Dividing it by \(f\) yields:

(5.36)\[\frac{s_{0x}'}{f} = \frac{s_{0x}}{s_{0z}} - \frac{x_0}{f}\]

From Fig. 5.9, it follows that \(\tan \varepsilon_x' = \frac{s_{0x}'}{f}\) and \(\tan \varepsilon_x = \frac{s_{0x}}{s_{0z}}\). Therefore, we have:

(5.37)\[\tan \varepsilon_x' = \tan \varepsilon_x - \frac{x_0}{f}\]

Analogously, in the y-direction, we obtain:

(5.38)\[\tan \varepsilon_y' = \tan \varepsilon_y - \frac{y_0}{f}\]

This angular representation is a formulation also found in [6].

5.2.6. Filtering

When passing through a filter, a ray with power \(P_i\) and wavelength \(\lambda_i\) is attenuated according to the filter’s transmission function \(T_\text{F}(\lambda)\):

(5.39)\[\begin{split}P_{i+1} = \begin{cases} P_{i}~ T_\text{F}(\lambda_i) & \text{for}~~ T_\text{F}(\lambda_i) > T_\text{th}\\ 0 & \text{else}\\ \end{cases}\end{split}\]

Additionally, ray powers are set to zero if the transmission falls below a specific threshold \(T_\text{th}\). This approach avoids ghost rays that continue to be propagated during raytracing but carry very little power. As their contribution to image formation is negligible, absorbing them as soon as possible helps speed up the tracing process.

As a side note, apertures are also implemented as filters, but with \(T_\text{F}(\lambda) = 0\) for all wavelengths.

5.2.7. Geometry Checks

Geometry checks before tracing include:

  • All tracing-relevant elements must be within the outline.

  • There must be no object collisions.

  • Elements must follow a defined, sequential order.

  • Ray sources must be available.

  • All ray sources must precede all other types of elements.

Collision checks are performed by first sorting the elements and then comparing positions on adjacent surfaces. This is implemented by randomly sampling many lateral surface positions and checking the correct order of the axial coordinates. While this method doesn’t guarantee the absence of collisions, the sequentiality is verified for each ray during raytracing, and warnings are emitted if issues are detected.

5.2.8. Intersection Calculation

5.2.8.1. Surface Extension

To simplify the handling of non-intersecting rays and ensure each ray has the same number of sections, the surface is extended so that all rays intersect a surface. Gaps on the surface are filled, and the surface edge is extended radially towards infinity.

An intersection is calculated for the extended surface. The ray is marked as hitting or non-hitting based on a surface mask afterwards.

../_images/surface_extension.svg

Fig. 5.10 Surface Extension

5.2.8.2. Intersection of a Ray with a Plane

The intersection of all planar surface types (CircularSurface, RectangularSurface, RingSurface) as well as the TiltedSurface can be computed analytically. The following definitions hold:

(5.40)\[\begin{split}\text{surface normal vector:}~~~~ \vec{n} &= (n_x, n_y, n_z)\\ \text{surface center vector:}~~~~ \vec{q} &= (x_0, y_0, z_0)\\ \text{point on ray or surface:}~~~~ \vec{p} &= (x, y, z)\\ \text{ray support vector:}~~~~ \vec{p_0} &= (x_{0p}, y_{0p}, z_{0p})\\\end{split}\]

The surface point normal equation is:

(5.41)\[(\vec{p} - \vec{q})\cdot \vec{n} = 0\]

The ray equation in dependence of ray parameter \(t\) is as follows:

(5.42)\[\vec{p} = \vec{p_0} + \vec{s} \cdot t\]

Inserting these equations into each other results in:

(5.43)\[(\vec{p_0} + \vec{s}\cdot t_\text{h} - \vec{q}) \cdot \vec{n} = 0\]

Rearranging provides the ray parameter for the hit position \(t_\text{h}\) for the case where \(\vec{s} \cdot \vec{n} \neq 0\):

(5.44)\[t_\text{h} = \frac{(\vec{q} - \vec{p_0})\cdot \vec{n}}{\vec{s} \cdot \vec{n}}\]

This result can be inserted into the ray equation to determine the intersection position. If \(\vec{s} \cdot \vec{n} = 0\), the surface normal and ray direction vector are perpendicular. If the ray lies in the plane, there are infinite solutions. If not, there is none.

5.2.8.3. Intersection of a Ray with a Conic Surface

Intersections with a ConicSurface can also be calculated analytically.

(5.45)\[\begin{split}\text{Ray support vector:}~~~~ \vec{p} &= (p_x, p_y, p_z)\\ \text{Ray direction vector:}~~~~ \vec{s} &= (s_x, s_y, s_z)\\ \text{Center of surface:}~~~~ \vec{q} &= (x_0, y_0, z_0)\end{split}\]

Equating the ray equation and the conic section leads to:

(5.46)\[p_z + s_z t = z_0 + \frac{\rho r^2}{1 + \sqrt{1-(k+1)\rho^2 r^2}}\]

The squared radius is:

(5.47)\[r^2 = (p_x + s_x t - x_0)^2 + (p_y+s_y t - y_0)^2\]

Some work in rearranging leads to the following equation, valid for \(\rho \neq 0\) (conic section is not a plane):

(5.48)\[A t^2 + 2 B t + C = 0\]

With the new parameters:

(5.49)\[\begin{split}A &= 1 + k s_z^2\\ B &= o_x s_x + o_y s_y - \frac{s_z}{\rho} + (k+1) o_z s_z\\ C &= o_x^2 + o_y^2 - 2\frac{o_z}{\rho} + (k+1) o_z^2\\ \vec{o} &= \vec{p} - \vec{q} = (o_x, o_y, o_z)\end{split}\]

The solutions for \(t\) are (\(\rho \neq 0\)):

(5.50)\[\begin{split}t = \begin{cases} \frac{-B \pm \sqrt{B^2 -CA}}{A} & \text{for}~~ A \neq 0, ~~ B^2 - CA \geq 0 \\ -\frac{C}{2B} & \text{for}~~ A = 0, ~~B \neq 0\\ \{\mathbb{R}\} & \text{for}~~ A = 0, ~~B = 0, ~~C = 0\\ \emptyset & \text{else} \end{cases}\end{split}\]

The first case produces two intersections, while in the second, the ray touches the surface tangentially. In the fourth case, there is no intersection. In cases one and two, it must be verified whether the intersection lies within the valid radial range of the surface, as the conic equation describes the surface for all possible \(r\). If there are still two valid intersections in the first case, the one with the lower axial position is used. The third case is only relevant for rays lying along a linear section of a surface, such as a plane or a cone surface. This is further explained below.

Equation (5.50) yields the same results as described in [7], although in his derivation, the constants \(A, B, C\) include an additional curvature scaling and the factor of 2 is incorporated inside \(B\).

In the case of \(\rho = 0\) (plane) equation (5.46) simplifies to \(p_z + s_z t = z_0\), resulting in:

(5.51)\[\begin{split}t = \begin{cases} - \frac{p_z - z_0}{s_z} & \text{for}~~ s_z \neq 0 \\ \{\mathbb{R}\} & \text{for}~~ s_z = 0, p_z = z_0\\ \emptyset & \text{for}~~ s_z = 0,~p_z\neq z_0\\ \end{cases}\end{split}\]

The first case produces exactly one intersection and corresponds to the special case of equation (5.44) for \(\vec{n} = (0, 0, 1)\). In the second case, the ray lies within the plane, resulting in infinite intersection positions. In the third case, the ray is parallel but offset from the plane, producing no intersections. We enforce \(s_z > 0\) in our simulator, making only the first case relevant.

For \(\rho \rightarrow \infty\) and \(r \in \mathbb{R}\), the conic section becomes a cone for \(k < -1\), a line from \(z = z_0\) to \(z \rightarrow \infty\) for \(k = -1\), or a point at \((x_0, y_0, z_0)\) for \(k > -1\). In the first two cases, the conic sections consist of long linear segments. Depending on the cases of equation (5.50), a ray can intersect once, twice, not hit the surface, or lie within it. For a point (\(k > -1\)), there can be no intersection or a single one. In the latter case, \(B^2 - CA = 0\) holds, producing a single intersection in the first case of equation (5.50).

5.2.9. Normal Calculation

5.2.9.1. General

Surface normals are required to calculate the refraction, polarization, and transmittance of a ray. The equation for an analytical, unnormalized normal vector is presented below, as detailed in [10].

(5.55)\[\begin{split}\vec{n} = \begin{pmatrix} -\frac{\partial z}{\partial x}\\ -\frac{\partial z}{\partial y}\\ 1\\ \end{pmatrix}\end{split}\]

To use this vector for calculations, it must be normalized:

(5.56)\[\vec{n}_0 = \frac{\vec{n}}{|| \vec{n} ||}\]

5.2.9.2. Planar Surfaces

For the surface types Rectangle, Circle, and Ring, the normal vector \(n\) is simply \(n = (0, 0, 1)\). For the surface type TiltedSurface, a user-specified value for \(n\) is provided, which is normalized internally.

5.2.9.3. Asphere

The radial derivative of an asphere equation is given by:

(5.57)\[m_r = \frac{\partial z}{\partial r} = \frac{\rho r}{\sqrt{1 - (k+1)\rho^2 r^2}} + \sum_{i=1}^{m} 2i \cdot a_i \cdot r^{2i - 1} = \tan \alpha\]

This radial component must be rotated around the vector \((0, 0, 1)\) by the angle \(\phi\). According to the general normal vector calculation, the x- and y-components of the unnormalized vector are computed as follows:

(5.58)\[\begin{split}\frac{\partial z}{\partial x} &= m_r \cos \phi\\ \frac{\partial z}{\partial y} &= m_r \sin \phi\\\end{split}\]

These components can be applied to the equation described above.

5.2.9.4. Conic Surface

The derivative of the conic function represents a simplified case of the asphere scenario, where polynomial terms are absent:

(5.59)\[m_r = \frac{\rho r}{\sqrt{1 - (k+1)\rho^2 r^2}}\]

Here, \(m_r\) is equivalent to \(\tan \alpha\), the tangent of the angle between the surface normal and the vector \((0, 0, 1)\). Using geometrical relations, one can derive \(\sin \alpha\), which leads to the radial component \(n_r\) of the normalized normal vector:

(5.60)\[n_r = -\sin{\alpha} = -\frac{m_r}{\sqrt{m_r^2+1}} = -\frac{\rho r}{\sqrt{1- k\rho^2 r^2}}\]

The normalized vector \(n\) for the ConicSurface is expressed as follows:

(5.61)\[\begin{split}\vec{n} = \begin{pmatrix} n_r \cos \phi\\ n_r \sin \phi\\ \sqrt{1- n_r^2} \end{pmatrix}\end{split}\]

For a SphericalSurface, this expression simplifies further to:

(5.62)\[\begin{split}\vec{n} = \begin{pmatrix} -\rho r \cos \phi \\ -\rho {}r \sin \phi\\ \sqrt{ 1 - \rho^2 r^2}\\ \end{pmatrix}\end{split}\]

5.2.9.5. FunctionSurface

For the surface types FunctionSurface1D and FunctionSurface2D, the derivatives required for constructing the normal vector are typically calculated numerically. An exception to this arises when a derivative function is directly provided to the object.

The central first derivative has the form:

(5.63)\[f'(x) = \frac{f(x+\varepsilon) - f(x-\varepsilon)}{2 \varepsilon}\]

For partial surface derivatives, the expressions are approximately given by:

(5.64)\[\begin{split}\frac{\partial z}{\partial x} &\approx \frac{z(x + \varepsilon, ~y) - z(x - \varepsilon, ~y)}{2\varepsilon}\\ \frac{\partial z}{\partial y} &\approx \frac{z(x, ~y + \varepsilon) - z(x, ~y - \varepsilon)}{2\varepsilon}\\\end{split}\]

For a FunctionSurface1D, both derivatives are the same due to symmetry, so only a single value needs to be computed. The step width \(\varepsilon\) is selected as the maximum of the optimal step width \(\varepsilon_\text{o}\) and the minimal positional step width \(\varepsilon_\text{p}\):

(5.65)\[\varepsilon = \max (\varepsilon_\text{o}, ~\varepsilon_\text{n})\]

The optimal step width, as derived from [11], is given by:

(5.66)\[\varepsilon_\text{o} = \sqrt[3]{3 \varepsilon_\text{f} \left| \frac{f(x)}{f^{(3)}(x)} \right|}\]

Here, \(\varepsilon_\text{f}\) represents the machine precision for the floating-point type being utilized. Under the assumption of predominantly spherical surfaces, the primary function component is \(x^2\). Higher polynomial orders are less significant, a reasonable assumption is \(\left| \frac{f(x)}{f^{(3)}(x)} \right| = 50\).

Although this may vary for different functions, due to the influence of the cube root, a quotient 1000 times larger modifies \(\varepsilon_\text{o}\) by only about a factor of \(10\). With a 64-bit floating-point number, where \(\varepsilon_\text{f} \approx 2.22 \times 10^{-16}\), the result is \(\varepsilon_\text{o} \approx 3.22 \times 10^{-5}\). Since optrace units are measured in millimeters, this corresponds to \(32.2\,\) nm.

In order to ensure that numerical differences in \(f\) are not only representable, but also that \(x + \varepsilon\) is distinct from \(x\), it is essential to satisfy the condition \(x + \varepsilon > x (1 + \varepsilon_\text{f})\) for every coordinate \(x\) on the surface.

Given \(r_\text{max}\) as the maximum absolute distance on the surface, the minimal bound can be expressed as:

(5.67)\[\varepsilon_\text{p} = r_\text{max} ~\varepsilon_\text{f}\]

It is advisable to center the surface at \(r = 0\) prior to differentiation to minimize \(r_\text{max}\).

5.2.9.6. DataSurface

A DataSurface object utilizes the scipy.interpolate.InterpolatedUnivariateSpline (FunctionSurface1D) or the scipy.interpolate.RectBivariateSpline (FunctionSurface2D) to interpolate surface values. To ensure smooth curvature transitions, an interpolation order of k=4 is applied in all lateral dimensions.

By invoking the spline objects with scipy.interpolate.InterpolatedUnivariateSpline.__call__() or scipy.interpolate.RectBivariateSpline.__call__(), it is possible to specify a derivative parameter, enabling the output of surface derivatives. These derivatives can then be used to compute surface normals, as described in the methods discussed earlier.

It is important to note that due to the nature of interpolation, the minimum and maximum surface values may exceed the specified data range. A warning informs the user about the occurence of large deviations.

References