## Abstract

Modeling techniques for light-shaping systems with freeform surface are presented from a physical-optics point of view. We apply the modeling techniques to different light-shaping systems with freeform surfaces designed by “ray mapping method”. The simulation results show that the design is not always valid. Diffraction effects occur, especially in paraxial situations. We discuss the accuracy of the design via physical-optics simulation, and find an explanation in the geometric-optics assumption of the design algorithm being sufficient only if the optical system results in homeomorphic behavior for the electric field between the input and target.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

For applications in illumination, e.g. indoor and outdoor lighting, freeform surfaces provide more freedom than traditional lenses for the design of the geometry of the optical structure, especially for non-symmetric illumination problems [1–4]. However, the optical design of freeform surfaces with the aim of shaping a light source into a specified energy distribution still poses important challenges [5]. One of the most often used design methods is the so-called “ray mapping” algorithm [6–11], which mainly contains two steps: (1) a one-to-one mapping is searched between the energy distribution of the input and that of the target; (2) a freeform surface is designed to realize that mapping. Another possible design method is the “Monge-Ampère equation method” [2,4,12,13], whereby a partial differential equation for the function representing the surface is set up, and numerical algorithms are applied to discretizing and solving said equation.

For both types of methods mentioned above, the geometric-optics assumption is always presupposed in the design of the freeform surface, which establishes a one-to-one mapping for the energy distribution between the input and the target, or in literature, also called the intensity law of geometrical optics [14]. On the other hand, for analyzing the design result, a ray-tracing simulation is performed, which gives a result corresponding exactly to the target, and is usually given as proof of the validity of the design. However, using a simulation technique which is also under the geometric-optics assumption is not enough for a full analysis of the accuracy of the design, since any eventual physical-optics effects would go unnoticed, both in the design and in the subsequent simulation. Therefore, to answer the question of whether the geometric-optics assumption is accurate enough for the design algorithm, the designed result has to be investigated via physical-optics simulation. A related question, the resolution limit of the usage of a freeform surface, was discussed by Zwick et al [15]. In our paper we focus on the limitation of the validity of geometrical optics in freeform modeling and design, and how this limitation can be mathematically expressed and investigated in the framework of physical optics. We also introduce a fast physical-optics modeling approach which enables 2D calculations, in contrast to the 1D simulations in Zwick et al [15].

In this work, the freeform surfaces for different tasks are designed, with the geometric-optics assumption applied. Then the field tracing techniques are introduced for analyzing the designed results. Together with the examples, the accuracy of the geometric-assumption in the design is discussed based on the modeling techniques, especially the application of different Fourier transform algorithms.

## 2. Task description and design

For illumination applications, the design task is usually to shape the light in its far field. Here, two simple light-shaping systems are shown in Fig. 1 for the purposes of discussion. A freeform lens with a square aperture is sought in order to shape the input field to a given irradiance distribution with two different sizes. In both cases, the freeform lens is configured with a fixed, predefined, planar surface, and the freeform surface to be designed. The N-BK7 glass is used for the lens, with the refractive index of $1.52$ at wavelength of $532 \, \mathrm{nm}$. The size of the target irradiance pattern for the first case is $900\times 900\, \mathrm{mm}$ (shown in Fig. 1(a)), while for the other one is smaller at $150\times 150\, \mathrm{mm}$ (shown in Fig. 1(b)). The source is taken to be a plane-wave field, with wavelength $532 \, \mathrm{nm}$, and field size $1\times 1\, \mathrm{mm}$. For a more straightforward discussion, the target plane is set perpendicular to the optical axis and located $1\, \mathrm{m}$ away from the freeform lens. We have selected the parameters for our examples to demonstrate the theoretical premise, which is the main object of this manuscript, without paying much attention to the manufacturability or otherwise of the device, since such considerations lie beyond the scope of this paper. Consequently, manufacturability would be an added consideration in the process of designing such an optical element, mostly uncoupled from the points analyzed in this paper.

The freeform surface is designed using the algorithm presented in this paper’s authors’ previous work [16], which can also be understood as a “ray mapping” method developed under the geometric-optics assumption. The algorithm assumes a homeomorphism, i.e. a bijective map between the coordinate of the irradiance function on the input plane and the target, with the input plane defined right in front of the lens. Therefore, the energy conservation law can be written locally in the form of

The search for a smooth, bijective mapping $\boldsymbol{\rho}^{\prime }(\boldsymbol{\rho}$) in Eq. (1) can be mathematically derived as the $L^{2}$ Monge-Kantorovich problem or optimal mass transport (OMT) problem. Numerical methods have been proposed to solve the problem [17–20]. Here, we consider the irradiance defined on 2D supports in $\mathbb {R}^{2}$, and implement the algorithm from Prins [20] with some modifications, since said algorithm is flexible enough to solve the mapping between two supports with different boundaries. In both light-shaping tasks, we use $301\times 301$ sampling points for calculating the mapping. The sampling points are also used later for the surface design.

Once the mapping is solved, an iterative algorithm is performed to design the surface structure [16]. An initial guess for the surface is necessary to start the algorithm. Then, in each iteration, with the mapping, the input wavefront and the current surface profile, the required normal vector of the surface is found by Snell’s law. Then a B-spline method is used to integrate the surface height function from its gradient, which can be concluded directly from the normal vector. The difference of the surface height value between two iterations is calculated, and the $L^{2}$ norm of the relative difference is used as a merit function. The algorithm is applied iteratively till the surface converges. Usually the surface converges fast within 5 iterations.

The height profiles of the designed freeform surface for both cases are shown in Fig. 2. The freeform surface is represented by a B-spline function. From the result, one can see that the profiles are similar in both cases, but have different gradients.

## 3. Modeling techniques

In order to analyze the validity of the design result, simulations are performed with a physical-optics modeling technique, field tracing, which combines various modeling methods for different subdomains of an optical system [21,22]. The resulting fields of the individual subdomains are then interconnected, sequentially or non-sequentially, providing a total solution for the whole optical system. In order to perform the simulation in an efficient way, modeling methods are applied either in the spatial domain or the spatial-frequency domain, depending on the mathematical characteristics of the method in question. Different algorithms for the calculation of the Fourier transform are available for transforming the field between the domains [23,24].

A field-tracing diagram, as shown in Fig. 3, illustrates the modeling sequence used in the light path of the light-shaping system. In Fig. 3, all the $\boldsymbol{E}(\boldsymbol{\rho}$) are the vectorial electric fields, i.e. $\boldsymbol{E}(\boldsymbol{\rho})=(E_x(\boldsymbol{\rho}), E_y(\boldsymbol{\rho}), E_z(\boldsymbol{\rho})$), defined in the $\boldsymbol{\rho}$ domain. $\boldsymbol{E}^{\mathrm {ln}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,f}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \right \vert_{\mathrm {P}}{^\mathrm {out,f}}$ are fields defined on the plane $\mathrm {P}{^\mathrm {in,p}}$ which is in front of the planar surface, and the plane $\mathrm {P}{^\mathrm {out,p}}$ behind the planar surface; $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,f}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \right \vert_{\mathrm {P}}{^\mathrm {out,f}}$ are fields defined on the plane $\mathrm {P}{^\mathrm {in,f}}$ in front of the freeform surface, and the plane $\mathrm {P}{^\mathrm {out,f}}$ behind it; $\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}') \left. \right \vert_{\mathrm {P}}{^\mathrm {tar}}$ is the field on the target plane. $\boldsymbol{\tilde{E}}^{\mathrm{out}}(\boldsymbol{\kappa}) \left \vert_{\mathrm {P}}{^\mathrm {out,f}}, \boldsymbol{\tilde{E}}^{\mathrm{out}}(\boldsymbol{\kappa}) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$ are the fields defined in the $\boldsymbol{\kappa}$ domain, with $\boldsymbol{\kappa}=(k_x, k_y)$, the $x\, \mathrm{and} \, y$ component of the wave vector.

$\boldsymbol {\mathcal{B}}$ denotes the LPIA operator for both surfaces of the freeform lens, which mathematically is in a matrix form [25]. LPIA is applied in the system for modeling the propagation of the field through the freeform lens,

Each instance of application of the LPIA operator $\boldsymbol {\mathcal{B}}$ contains mainly three sequential steps [25]: 1. a coordinate transformation for the field on the surface from the global coordinate system to the (position-dependent) surface coordinate system; 2. applying the also position-dependent Fresnel coefficients to the local field components; 3. a coordinate transformation for the field on the surface from the local coordinate system back to the global one. LPIA is a pointwise operator for the electric field, which gives homeomorphism behavior between $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,p}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \right \vert_{\mathrm {P}}{^\mathrm {out,p}}, \boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,f}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \right \vert_{\mathrm {P}}{^\mathrm {out,f}}$. For the example shown in Fig. (3), Eq. (2a) leads to an aperture effect only. However, in the case of a curved first surface and/or a non-planar input field, Eq. (2a) enables the propagation of the input field through the first surface including all relevant effects.

We would like to mention that LPIA is also used by Zwick et al in their derivation of Eq. (4) in [15]. However, the authors assume small angles and neglect boundary effects on amplitude and polarization. Thus, they use the input field on the surface itself in the propagation integral formula in space domain (see Eq. (3) in [15]). We use a full LPIA operator, including amplitude and polarization effects, which enables modeling situations with large angles. Moreover, we replace the convolution-type integral formula by its fast implementation in the $\boldsymbol{\kappa}$ domain, by applying different algorithms for the calculation of the Fourier transforms [23,24,26]. This is discussed next.

$\boldsymbol{\tilde{\mathcal{P}}}$ in Fig. 3 is the free space propagation operator. Regarding the free-space propagation inside the freeform lens medium and in air behind the freeform lens, the field tracing is a combination of Fourier transforms and a propagation operator in the $\boldsymbol{\kappa}$ domain. In the following, the free-space propagation behind the freeform lens is taken as an example to illustrate the algorithm. As shown in Fig. 3, the field at the target plane is obtained by

${\boldsymbol {\mathcal {F}}}\, \mathrm{and} \, {\boldsymbol {\mathcal {F}}}^{-1}$ in the field tracing diagram are the direct and inverse Fourier operators for transforming the fields between two domains. Numerically, there are several algorithms available to compute the Fourier transform in Eq. (3) [23,24], e.g. Fast Fourier Transform (FFT), semi-analytical Fourier Transform (SFT), or Homeomorphic Fourier Transform (HFT).

HFT and inverse HFT are approximated pointwise Fourier transforms, which produce a one-to-one map between the coordinates of the fields expressed in the $\boldsymbol{\rho}$ domain and in the $\boldsymbol{\kappa}$ domain [23]. Specifically, we denote the electric field in the space domain by

With the above in mind, the HFT operation is formulated as follows:

HFT (or inverse HFT) is purely a mathematical tool which performs the Fourier transform via a pointwise calculation, whose accuracy is mathematically determined by the validity of the stationary-phase approximation and the bijective characteristic of the wavefront phase terms $\psi (\boldsymbol{\rho}$) (or ${\tilde {\psi }}(\boldsymbol{\kappa}$)). If the direct (or the inverse) Fourier transforms for the electric field in Fig. 3 can be approximated by the HFT (or the inverse HFT), we say that the electric field is in its geometric zone.

Therefore, if all the direct and the inverse Fourier transforms in the field tracing can be approximated by the HFT and the inverse HFT, since the $\boldsymbol {\mathcal{B}}\, \mathrm{and} \, \boldsymbol{\tilde{\mathcal{P}}}$ are also pointwise operators, the whole field tracing simulation of the light-shaping system gives a pointwise calculation, which forms the homeomorphism between $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,p}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}^{\prime }) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$, i.e. $\boldsymbol{\rho}^{\prime }(\boldsymbol{\rho}$) is a bijective mapping function. And that is the typical assumption usually made in freeform surface design algorithms.

## 4. Simulation results and discussion

The ray tracing simulation for the freeform surface of case 1 (Fig. 1(a)) is shown in Fig. 4(a) – source plane and (b) – target plane. The ray positions are connected to illustrate the deformation of the input mesh by the freeform surface. It also shows that the freeform surface slightly deforms the square shape of the input. This could be avoided by a further optimization step. For this paper it is of no concern, since we are interested in investigating the validity of the mapping assumption itself, which is to be studied next.

In order to verify the assumption in the design algorithm is valid for both examples in Section 2, it is to be analyzed whether $\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {out,p}}, \boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \right \vert_{\mathrm {P}}{^\mathrm {in,f}}$, $\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {out,f}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}^{\prime }) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$ are in their geometric zone, or in another words, if the HFT or the inverse HFT can be applied accurate to perform for all the direct or inverse Fourier transforms in the field tracing diagram. Therefore, the system with the designed freeform surface is simulated firstly with HFT and inverse HFT. And we investigate their accuracy by comparing the simulation result against the one obtained with the rigorous FFT and inverse FFT.

In case 1, ${\boldsymbol {\mathcal {F}}}\, \mathrm{and} \, {\boldsymbol {\mathcal {F}}}^{-1}$ in Fig. 3 are selected as HFT and inverse HFT to start with. Therefore, a pointwise calculation between $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,p}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}^{\prime }) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$ is performed in the whole system simulation. The irradiance at the target plane is measured in the detector, and shown in Fig. 4(c). The result in Fig. 4(c) coincides with the target pattern in Fig. 1, however, as discussed, that is not sufficient proof of the validity of the designed surface. Next, the Fourier transform for $\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {out,p}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}) \right \vert_{\mathrm {P}}{^\mathrm {out,f}}$ are switched to the FFT, and inverse Fourier transform for calculating $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,f}}\, \mathrm{and} \, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}^{\prime }) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$ is switched to an inverse FFT. The detected irradiance at the target plane is shown in Fig. 4(d). The resulting irradiance distribution also mostly coincides with the target one, with the exception of only slight errors which occur at the edges of the pattern, small enough to be neglected. Therefore, the validity of the designed freeform surface for the task is proved. Both techniques give coincident results, that shows the HFT and the inverse HFT are good approximate algorithms for the direct and inverse Fourier transform in this case. Consequently, the full sequence from $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left \vert_{\mathrm {P}}{^\mathrm {in,p}}\, \mathrm{to} \,\boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}^{\prime }) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$ exhibits a homeomorphic behavior, so that the geometric-optics assumption $\boldsymbol{\rho}^{\prime }(\boldsymbol{\rho}$) is quite a good approximation for the design in this example.

For case 2 (Fig. 1(b)), simulations with the homeomorphic assumption and the rigorous calculation are performed as in the previous example. The resulting irradiance distribution on the target plane is shown in Fig. 5(a) for the simulation with the HFT and the inverse HFT, and in Fig. 5(b) for the one with the FFT and the inverse FFT. The result in Fig. 5(a) again coincides with the target pattern as expected. However, the other one shows obvious diffraction effects, and also with low resolution in the irradiance distribution, which together produce a large deviation from the target of the design task. The blurring of the result in Fig. 5(b) is in nice agreement with the conclusions of Zwick et al [15], and provides a 2D demonstration. Figure 5(a) and Fig. 5(b) together show that the simulation with the HFT and the inverse HFT in this case is far from a good approximation, and that reveals that there is no homeomorphism between $\boldsymbol{E}^{\mathrm {in}}(\boldsymbol{\rho}) \left\vert_{\mathrm {P}}{^\mathrm {in,p}} \,\mathrm{and}\, \boldsymbol{E}^{\mathrm{out}}(\boldsymbol{\rho}^{\prime }) \right \vert_{\mathrm {P}}{^\mathrm {tar}}$. Therefore, the freeform surface which is designed via the geometric-optics assumption cannot fulfill the requirements of the light-shaping task.

With the two examples in hand, we would like to point out that for those design algorithms which include the homeomorphic assumption for the field (or irradiance distribution) between the input and target planes, such as “ray mapping” design algorithms, the accuracy of the assumption should be determined by physical-optics simulation with the designed result. If the analysis is done by the geometric-optics based “ray tracing” technique, the simulation result will evidently always coincide with the target pattern for the design, since both design and analysis follow the same assumption, and additional physical effects which might appear in reality, such as diffraction, will go unnoticed. As a rule of thumb, the more non-paraxial (both the input and output side) the system, the more accurate the geometric-optics assumption.

The simulations were done with the software VirtualLab Fusion [28]. One physical-optics simulation on a PC (2.39 GHz Intel Core Quad CPU), in the current implementation, takes a few seconds. It is worth mentioning that, in field tracing, the sampling of the field for the propagation typically dominates the calculation time and this sampling is not related to the selected resolution in the detector.

## 5. Conclusion

As a conclusion, algorithms working under the homeomorphic assumption provide a fast way to design freeform surfaces for light shaping. However, from a physical-optics point of view, the critical point for determining the validity of the design result is whether the fields through the light-shaping system are in their geometric zone. The accuracy of the HFT applied in the field tracing can help to analyze the field zone. If the fields are not in their geometric zone, the light-shaping task cannot be fulfilled by a freeform component designed with those algorithms. A further optimization process should be applied for the resulting structure.

## Acknowledgments

The authors would like to thank our colleague Ms. Olga Baladron Zorita for revising the manuscript. We are also grateful to the company Wyrowski Photonics GmbH for procuring a VirtualLab Fusion [28] license free of charge for the simulations included in this work.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **K. P. Thompson and J. P. Rolland, “Freeform optical surfaces: A revolution in imaging optical design,” Opt. Photonics News **23**(6), 30–35 (2012). [CrossRef]

**2. **H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A **19**(3), 590–595 (2002). [CrossRef]

**3. **Z. Feng, L. Huang, G. Jin, and M. Gong, “Designing double freeform optical surfaces for controlling both irradiance and wavefront,” Opt. Express **21**(23), 28693–28701 (2013). [CrossRef]

**4. **R. Wu, K. Li, P. Liu, Z. Zheng, H. Li, and X. Liu, “Conceptual design of dedicated road lighting for city park and housing estate,” Appl. Opt. **52**(21), 5272–5278 (2013). [CrossRef]

**5. **R. Wu, Z. Feng, Z. Zheng, R. Liang, P. Benítez, J. C. Miñano, and F. Duerr, “Design of freeform illumination optics,” Laser Photonics Rev. **12**(7), 1700310 (2018). [CrossRef]

**6. **A. Bruneton, A. Bäuerle, R. Wester, J. Stollenwerk, and P. Loosen, “High resolution irradiance tailoring using multiple freeform surfaces,” Opt. Express **21**(9), 10563–10571 (2013). [CrossRef]

**7. **J. Chen, Z. Huang, T. Liu, M. Tsai, and K. Huang, “Freeform lens design for light-emitting diode uniform illumination by using a method of source-target luminous intensity mapping,” Appl. Opt. **54**(28), E146–E152 (2015). [CrossRef]

**8. **Z. Feng, B. D. Froese, and R. Liang, “Freeform illumination optics construction following an optimal transport map,” Appl. Opt. **55**(16), 4301–4306 (2016). [CrossRef]

**9. **C. Bosel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express **24**(13), 14271–14282 (2016). [CrossRef]

**10. **C. Gannon and R. Liang, “Ray mapping with surface information for freeform illumination design,” Opt. Express **25**(8), 9426–9434 (2017). [CrossRef]

**11. **K. Desnijder, P. Hanselaer, and Y. Meuret, “Flexible design method for freeform lenses with an arbitrary lens contour,” Opt. Lett. **42**(24), 5238–5241 (2017). [CrossRef]

**12. **K. Brix, Y. Hafizogullari, and A. Platen, “Designing illumination lenses and mirrors by the numerical solution of monge ampare equations,” J. Opt. Soc. Am. A **32**(11), 2227–2236 (2015). [CrossRef]

**13. **S. Chang, R. Wu, L. An, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a monge-ampare type equation,” J. Opt. **18**(12), 125602 (2016). [CrossRef]

**14. **M. Born and E. Wolf, * Principles of Optics* (Cambridge University Press, 1999), 7th ed.

**15. **S. Zwick, R. Feßler, J. Jegorov, and G. Notni, “Resolution limitations for tailored picture-generating freeform surfaces,” Opt. Express **20**(4), 3642–3653 (2012). [CrossRef]

**16. **L. Yang, R. Knoth, C. Hellmann, and F. Wyrowski, “Non-paraxial diffractive and refractive laser beam shaping,” Proc. SPIE10518 (2018).

**17. **S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, “Optimal mass transport for registration and warping,” Int. J. Comput. Vis. **60**(3), 225–240 (2004). [CrossRef]

**18. **R. Wu, Y. Zhang, M. M. Sulman, Z. Zheng, P. Benitez, and J. C. Minano, “Initial design with l2 monge-kantorovich theory for the monge–ampare equation method in freeform surface illumination design,” Opt. Express **22**(13), 16161–16177 (2014). [CrossRef]

**19. **M. M. Sulman, J. Williams, and R. D. Russell, “An efficient approach for the numerical solution of the monge-ampere equation,” Appl. Numer. Math. **61**(3), 298–307 (2011). [CrossRef]

**20. **C. Prins, R. Beltman, J. ten Thije Boonkkamp, W. IJzerman, and T. Tukker, “A least-squares method for optimal transport using the monge–ampere equation,” SIAM J. Sci. Comput. **37**(6), B937–B961 (2015). [CrossRef]

**21. **M. Kuhn, F. Wyrowski, and C. Hellmann, “Non-sequential optical field tracing, in Advanced Finite Element Methods and Applications,” vol. 66 of * Lecture Notes in Applied and Computational Mechanics*T. Apel and O. Steinbach, eds. (Springer Berlin Heidelberg, 2013), pp. 257–273.

**22. **F. Wyrowski, “Unification of the geometric and diffractive theories of electromagnetic fields,” Proc. DGaO (2017).

**23. **O. Baladron-Zorita, Z. Wang, C. Hellmann, and F. Wyrowski, “Isolating the Gouy phase shift in a full physical-optics solution to the propagation problem,” J. Opt. Soc. Am. A **36**(9), 1551–1558 (2019). [CrossRef]

**24. **Z. Wang, S. Zhang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Application of the semi-analytical fourier transform to electromagnetic modeling,” Opt. Express **27**(11), 15335–15350 (2019). [CrossRef]

**25. **R. Shi, C. Hellmann, and F. Wyrowski, “Physical-optics propagation through curved surfaces,” J. Opt. Soc. Am. A **36**(7), 1252–1260 (2019). [CrossRef]

**26. **Z. Wang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Theory and algorithm of the homeomorphic fourier transform for optical simulations,” Opt. Express **28**(7), 10552–10571 (2020). [CrossRef]

**27. **O. Bryngdahl, “Optical map transformations,” Opt. Commun. **10**(2), 164–168 (1974). Test [CrossRef]

**28. **Physical optics simulation and design software “Wyrowski VirtualLab Fusion”, developed by Wyrowski Photonics GmbH, distributed by LightTrans International UG, Jena, Germany.